Next Article in Journal
Two Modified Chaotic Maps Based on Discrete Memristor Model
Next Article in Special Issue
Symmetric Mass Generation
Previous Article in Journal
DII-GCN: Dropedge Based Deep Graph Convolutional Networks
Previous Article in Special Issue
Symmetries of Thirring Models on 3D Lattices
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Large-S and Tensor-Network Methods for Strongly-Interacting Topological Insulators

1
Scuola Internazionale Superiore di Studi Avanzati, Via Bonomea 265, 34136 Trieste, Italy
2
Department of Mathematical Sciences, University of Liverpool, Liverpool L69 3BX, UK
3
Instituto de Física Teórica, Consejo Superior de Investigaciones Científicas, Universidad Autónoma de Madrid, Cantoblanco, 28049 Madrid, Spain
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(4), 799; https://doi.org/10.3390/sym14040799
Submission received: 18 March 2022 / Revised: 4 April 2022 / Accepted: 6 April 2022 / Published: 12 April 2022
(This article belongs to the Special Issue New Applications of Symmetry in Lattice Field Theory)

Abstract

:
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 discretisation 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 discretisations leads to a novel Heisenberg–Ising compass model with critical lines that flow with the value of r.

1. Introduction

1.1. Topological Matter and Relativistic Field Theories

Our most accurate description of nature is based on a four-dimensional quantum field theory (QFT) of fermionic matter coupled with gauge fields: the standard model of particle physics [1]. In this context, 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 nonperturbative 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, and 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 nonlinear constraint. These models serve to develop and test tools, such as bosonisation [6] and large-N expansions [7], the predictions of which can be benchmarked with efficient numerical methods for 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 do not solve specific real problems in high-energy physics that can be falsified experimentally. Instead, within high-energy physics, these toy QFTs act as theoretical laboratories.
Remarkably, during the last decades, we have witnessed a change in status for low-dimensional QFTs. Rather than looking at high energies and small-length scales, one may instead focus on nonrelativistic condensed-matter systems which, at long wavelengths and small energies [8], display certain universal behaviours determined by emergent relativistic QFTs. These effective descriptions [9] provide a more flexible framework in comparison to that provided by 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 have emerged include graphene [10], Weyl semimetals [11], and topological insulators and superconductors [12]. Regarding the emergence of scalar fields subjected to nonlinear constraints, quantum magnets play a prominent role [13].
A further step in this direction was provided by the so-called quantum simulators (QSs) [14]: well-isolated quantum many-body 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 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 lattice field theory approach [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 greater than that spacing, and the continuum description sets in. In recent years, we have taken very promising experimental steps in this direction for Dirac fermions [24,25], 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 the QFTs of self-interacting Dirac fermions. As discussed below, these QFTs were inspired by the Thirring and Gross–Neveu models, the origin of which can be traced back to the seminal contributions of E. Fermi [37,38] and Y. Nambu and G. Jona-Lasinio [39,40]. In particular, we explore specific lattice discretisations based on the so-called Wilson fermions [41], which make direct connections between these four-Fermi QFTs and the aforementioned topological insulators [12,42,43,44,45,46,47,48], allowing us to study the effect of electron–electron interactions. We note that recent advances in cold-atom 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, provide motivation for the 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 nonperturbatively [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 discussed in a series of recent works, the native Hubbard interactions [55,56] of cold-atom QSs of spin-orbit 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 can be regularised via a Wilson-type discretisation. Such a discretisation 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 unclear as one switches on the interactions. Here, as a consequence of spontaneous symmetry breaking, one can find phases of 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 nonlinear sigma model with discrete Z 2 symmetry, a constrained QFT amenable to a large-S expansion in the limit where the effective spin is 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 nonperturbative 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 tunnellings, each of which can be independently controlled by the lasers that form a so-called optical Raman lattice [29,30]. Accordingly, reaching the regime of r = 1 requires 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 interesting and useful to understand the effects of 0 < r < 1 . It is worth mentioning that, for asymptotically-free theories where the continuum limit is recovered from the lattice model in the regime of weak coupling strengths g 2 , one can generally expect that the value of Wilson parameter r will not change the properties of the continuum QFT. From the perspective of condensed matter and topological insulators, however, we are not only interested in this weak-coupling regime but in the whole phase diagram where the specific value of r can modify the layout and even allow for new strong-coupling phases. The goal of the present work is to explore this possibility.

1.2. 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
H 0 = 1 2 ( Π ( x ) · Π ( x ) j Φ ( x ) · j Φ ( x ) ) ,
where we use natural units = c = 1 and Einstein’s convention of repeated-index summation. Here, the fields and conjugate momenta Π ( t , x ) = t Φ ( t , x ) fulfil the canonical algebra [ Φ f 1 ( t , x 1 ) , Π f 2 ( t , x 2 ) ] = i δ f 1 , f 2 δ d ( x 1 x 2 ) , and j { 1 , , d } labels the spatial coordinates of a D = ( d + 1 ) -dimensional Minkowski spacetime with the metric η = diag ( 1 , 1 , , 1 ) . This QFT describes N uncoupled scalar bosons and is invariant under a continuous internal symmetry Φ ( x ) , Π ( x ) o Φ ( x ) , o Π ( x ) , where o O ( N ) is an arbitrary rotation. In order to couple the different flavours, a quartic self-interaction that respects this internal symmetry can be introduced:
H int = λ 0 4 ! Φ ( x ) · Φ ( x ) Φ 0 2 2
where λ 0 is the bare coupling strength. Here, we introduce Φ 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 σ [1] with σ ( x ) = Φ 1 ( x ) / Φ 0 , π ( x ) = ( Φ 2 ( x ) , Φ N ( x ) ) t / Φ 0 corresponding to the symmetry-breaking and Goldstone components, respectively. Instead of expanding around the SSB ground state, one may focus on the strong-coupling limit λ 0 , where the ground-state minimises the interaction energy (2) by imposing a nonlinear constraint on the fields
σ 2 ( x ) + π 2 ( x ) = 1 ,
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 ) nonlinear sigma model shares important features with non-Abelian gauge theories, such as asymptotic freedom for N > 2 [5], the existence of topologically non-trivial solutions called instantons [64], or large-N methods and dimensional transmutation [7,65].
In this work, we discuss how similar constraints can also appear 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
H 0 = Ψ ¯ ( x ) i ( I N γ j ) j Ψ ( x ) ,
where we introduce 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 discussion around Equation (2), we can now couple the flavours via a four-Fermi [37,38,39,40] contact interaction
H int = g 2 2 N Ψ ¯ ( x ) Ψ ( x ) 2 ,
where g 2 is the bare coupling strength. Although not directly apparent, as in the bosonic case, we show below that the strong-coupling limit g 2 leads to a constraint similar to the one of Equation (3), where the σ and π fields are 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 introduce a particular lattice discretisation, which plays an important role in determining the specific Z 2 SSB and its connection with topological insulators.

1.3. 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 Λ d , focusing on cases d = 1 , 2
Λ d = { j = 1 d n j a j e j : n j Z N j } ,
where { a j } are the lattice spacings along the { e j } unit vectors, and N j are the corresponding number of lattice sites along each axis (see Figure 1). Let us also note that for d = 1 , 2 spatial dimensions, one can use the following irreducible representations of the gamma matrices
d = 1 , γ 0 = σ z , γ 1 = i σ y , d = 2 , γ 0 = σ z , γ 1 = i σ y , γ 2 = i σ x ,
which are proportional to the Pauli matrices.
Discretising the spatial derivatives that appear in Equation (4) using central differences leads to the appearance of 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:
H = j = 1 d Ψ ¯ ( x ) i ( I N γ j ) 2 a j + r ( I N I s ) 2 a j Ψ ( x + a j e j ) + Ψ ¯ ( x ) m 2 d + r 2 a j Ψ ( x ) + H . c . g 2 2 N ( Ψ ¯ ( x ) Ψ ( x ) ) 2 ,
where we introduce the bare mass m and the aforementioned dimensionless Wilson parameter r. In Figure 1, we present a schematic diagram of this Wilsonian discretisation 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 = μ = 0 d a μ n μ e μ , such that recovering the time-continuum limit requires temporal anisotropies that permit the limit a 0 0 . In the case where there is no 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 being 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 of r = 1 brings the technical advantage that tunnelling terms are proportional to projection operators P j ± 1 2 ( 1 ± i γ j ) with P j ± 2 = 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. We 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 S O ( 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 nonunity Wilson parameters and anisotropic lattice constants, we start by focusing on the noninteracting limit g 2 = 0 . In this case, it is straightforward to compute the half-filled groundstates | ϵ ( k ) corresponding to the Dirac vacua [58,60,61], where we introduce the quasi-momentum within the first Brillouin zone k BZ = × j π / a j , π / a j . Associated with this band structure, one finds a Berry connection A ( k ) = ϵ ( k ) | i k | ϵ ( k ) [72,73], which characterises the principal fibre bundle associated with 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 the Wilson loop for a cycle in momentum space to be defined
W Z = e i φ Z , φ Z = d k A ( k ) = N π θ ( 2 r + m a 1 ) θ ( m a 1 ) ,
where θ ( x ) is the Heaviside step function. Therefore, one finds a trivial band insulator with W Z = + 1 for m a 1 > 0 or m a 1 < 2 r . Alternatively, a nontrivial topological insulator W Z = 1 arises when 2 r < m a 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 noninteracting phase diagram is completely equivalent if one simply rescales the dimensionless mass with the Wilson parameter m a 1 m a 1 / r .
For d = 2 , the Berry curvature B ( k ) = k A ( k ) [72] allows us to define the first Chern number
N Ch = d 2 k 2 π B ( k ) = N θ ( 2 r ξ 2 + m a 1 ) θ ( m a 1 ) + θ ( 2 r ( 1 + ξ 2 ) + m a 1 ) θ ( 2 r + m a 1 ) .
where ξ 2 = a 1 / a 2 is an anisotropy ratio, and we assume that ξ 2 1 . Here, some comments are in order. In the isotropic limit ξ 2 = 1 , when setting r = 1 , the ground state corresponds to quantum anomalous Hall ( QAH ) phase with N Ch = N for 2 < m a 1 < 0 and N Ch = + N for 4 < m a 1 < 2 , whereas it is a trivial band insulator with N Ch = 0 for m a 1 > 0 or m a 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 m a 1 axis comprising the QAH phases, surrounded by trivial band insulators at both sides. It is worth noting that both the QAH and 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 nonzero topological invariants 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], by 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 < m a 1 < 2 ξ 2 , which separates 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 nonunity Wilson parameter 0 < r < 1 is rather trivial in the noninteracting case—one can simply rescale m a 1 m a 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 that are different from the above trivial and topological insulators. In the following sections, we extend our previous studies presented in [58,60,61]. We study 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 is exploited to predict the shape of the phase diagram and possible phase transitions.

2. Strong Couplings and Effective Spin Models

2.1. 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 Equation (8) can be rewritten as
H int = g 2 ψ f , 1 ( x ) ψ f , 1 ( x ) ψ f , 2 ( x ) ψ f , 2 ( x ) ,
up to an irrelevant shift in 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 ground state will have a single fermion per site n , which has the freedom to select one of the two possible spinor configurations | n , | n . Within this subspace, the operators that fulfil the S U ( 2 ) algebra [ S ^ a ( t , x ) , S ^ b ( t , x ) ] = i δ n , n c ϵ abc S ^ c ( t , x ) at spatial points x = j n j a j e j , x = j n j a j e j become
S ^ ( x ) = S j a j ψ f ( x ) σ ψ f ( x ) S ( x ) = S σ n .
Here, S = 1 / 2 , and σ n is an operator acting on the projected Hilbert space H eff = n span { | n , | n } , that is defined by the tensor product of the identity I 2 on all sites except for n , where one can apply the vector of Pauli matrices σ = ( σ x , σ y , σ z ) . As discussed in [58,60,61] for d = 1 , 2 and for the 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 onsite terms and interact with each other via nearest-neighbour couplings, as depicted in Figure 2. The physical mechanism underlying these nearest-neighbour spin–spin interactions is the so-called super-exchange [83,84], and the most-general effective Hamiltonian can be written as follows
H eff = x Λ d a j = 1 d J j , a S a ( x ) S a ( x + a j e j ) + h a S a ( x ) .
Here, we use the label a { x , y , z } to distinguish the internal spin components from the spatial coordinates j { 1 , , d } , introduce a set of spin–spin couplings J j , a describing the strength of the interactions between neighbouring spins connected by a e j link and couple their internal spin components via a S a S a interaction (see Figure 2).
In previous works [58,60,61], where we used the standard choice of r = 1 , the nature of the spin–spin couplings was restricted to be of the Ising type. For d = 1 , where the coupling strength g 2 is dimensionless, the spin couplings found were
J 1 , a = 2 g 2 a 1 δ a , y , h a = 2 m + 1 a 1 δ a , z ,
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 S ( x ) ( S x ( x ) , S y ( x ) , S z ( x ) ) , x Λ 1 . This symmetry can also be combined with a reflection about the lattice centre, such that
S ( x ) ( S x ( x ) , S y ( x ) , S z ( x ) ) .
It is interesting to note that, for the representation of the gamma matrices in Equation (7), this Z 2 symmetry corresponds to
ψ f ( t , x ) γ 0 ψ f ( t , x ) ,
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 parity-breaking pseudo-scalar π condensate
S y ( x ) Π 5 = ψ ¯ f ( x ) i γ 5 ψ f ( x ) ,
where γ 5 = γ 0 γ 1 is the chiral gamma matrix. Note that, in the quantum Ising model [85], the ground state always displays nonzero magnetisation along the direction of the transverse field for any h z 0 . In the language of Dirac spinors, this corresponds to a nonzero value of the so-called scalar σ condensate
S z ( x ) Σ = ψ ¯ f ( x ) ψ f ( x ) ,
and the fact that it is generically nonzero 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 nonzero 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 this condensate nicely, even in the N = 1 limit and, moreover, provides an analytical expression for the critical line that has been proven to be very accurate when compared with matrix-product-state numerics [58]. In this manuscript, we 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
J 1 , a = 2 a 2 g 2 a 1 δ a , y , J 2 , a = 2 a 1 g 2 a 2 δ a , x , h a = 2 m + 1 a 1 + 1 a 2 δ a , z ,
which, again, have the correct units of energy, as shown in Equation (14). The corresponding spin model (13) is an instance of the so-called 90° compass model [88,89] with directional spin couplings S y S y ( S x S x ) along the e 1 ( e 2 ) spatial axis and a transverse magnetic field that is 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 ) | > 0 for | J 2 , x | > | J 1 , y | , achieved for a 1 > a 2 , and a FM y characterised by | S y ( 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 S O ( 1 , 2 ) . Rather than breaking parity, a nonzero value of the corresponding fermion π condensates
S x ( x ) Π 1 = ψ ¯ f ( x ) γ 1 ψ f ( x ) , S y ( x ) Π 2 = ψ ¯ f ( x ) γ 2 ψ f ( x ) ,
breaks the 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 to in [60,61] as Lorentz-breaking condensates, which contrast the parity-breaking pseudo-scalar condensate (17) of d = 1 .
In contrast, the regime with a nonvanishing transverse field h z 0 has not been studied in as 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 ) | > 0 ( | S 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 nonzero anisotropies ξ 2 > 0 that are finite ξ 2 1 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 introduce a parameter z e t a , which serves to locate the critical point and is equal to ζ = 2 in the regime of large spatial anisotropies. For other finite and nonzero anisotropies ξ 2 = a 1 / a 2 , the critical point will change, and ζ 2 will no longer be found exactly but must be estimated using numerical or analytical approximations [60,61].

2.2. Heisenberg–Ising Chains for d = 1

So far, we have reviewed known results that apply for the Wilson parameter r = 1 . Moving away from this limit modifies the super-exchange mechanism, leading to the additional spin–spin couplings depicted in Figure 2a. These still have the general form of Equation (13) but lead to different strengths with respect to those expressed in Equations (14)–(19). In particular, for d = 1 , we find that the expression of the external field is
h a = 2 m + r a 1 δ a , z
whereas the spin–spin couplings following the super-exchange mechanism of virtual double occupancies are now
J 1 , x = 1 r 2 g 2 a 1 , J 1 , y = 1 + r 2 g 2 a 1 , J 1 , z = 1 r 2 g 2 a 1 .
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 transformation that takes S y ( x ) ( 1 ) n 1 S y ( x ) . Under this transformation, the SSB ferromagnetic ground state is transformed into a classical Néel pattern of alternating spins. The discussion of the possible SSB orderings for the 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 and 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 signs, these can always be inverted by a spin rotation along the remaining axis with an alternating angle. The specific transformation in this case takes S ( x ) ( S x ( x ) , ( 1 ) n 1 S y ( x ) , ( 1 ) n 1 S z ( x ) ) , such that J 1 , a J = 1 / g 2 a > 0 , a = { x , y , z } , and the spin–spin interactions clearly display the S U ( 2 ) symmetry of the Heisenberg model. Additionally, if the external field is nonzero, this maps onto a staggered transverse field under the above transformation, such that
H eff H ^ eff = x Λ 1 J S ( x ) · S ( x + a 1 e 1 ) + h z e i k s · x S z ( x ) ,
where we introduce the wavevector k s = π a 1 e 1 . This transformation unveils a continuous U ( 1 ) symmetry with respect to rotations along the internal z axis that 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 staggered-fermion 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 the vanishing transverse field h z = 0 , the Heisenberg chain can be exactly solved via the Bethe ansatz [100,101] and does not support a long-range order, as shown via the inverse scattering method [102].
Given the clear difference between the r = 1 and r = 0 limits, one may expect different ground states with distinct magnetic orders, i.e., fermion condensates, for Wilson parameters 0 < r < 1 . In this regime, we recall that the noninteracting phase diagram comprises regions of nontrivial 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 U ( 1 ) symmetry by combining a pair of transformations. To do this, first, the spins are rotated about the internal z axis in an alternate fashion as S ( x ) ( ( 1 ) n 1 S x ( x ) , ( 1 ) n 1 S y ( x ) , S z ( x ) ) , which effectively changes the spin–spin couplings to J 1 , x = J 1 , z = J , and J 1 , y = J Δ , where
J = r 2 1 g 2 a 1 , Δ = 1 + r 2 1 r 2 .
Next, a rotation about the internal x axis S ( x ) ( S x ( x ) , S z ( x ) , S y ( x ) ) is applied; the spin chain maps onto the XXZ model [97,103,104], also known as a Heisenberg–Ising model, under an additional longitudinal field
H eff H ^ eff = J x Λ 1 ( S x ( x ) S x ( x + a 1 e 1 ) + S y ( x ) S y ( x + a 1 e 1 ) + Δ S z ( x ) S z ( x + a 1 e 1 ) + g y S y ( x ) ) .
As shown previously, for g y = h z / J = 0 , there is 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 S U ( 2 ) -symmetric point Δ = 1 [107], separating a critical phase at Δ < 1 from an Ising SSB 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 a 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 the FM y order, i.e., a 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 to those discussed for the quantum Ising model (14) may still appear, albeit at critical points that flow with the Wilson parameter.

2.3. Heisenberg–Ising Compass Models for d = 2

Before checking the validity of the above conjecture, let us discuss the effect of nonunit Wilson parameters on the effective spin model for the d = 2 case depicted in Figure 2b. Instead of the microscopic couplings presented in Equation (19), the super-exchange for 0 < r < 1 leads to an external field
h a = 2 m + r a 1 + r a 2 δ a , z ,
whereas the spin–spin couplings transform into
J 1 , x = a 2 ( 1 r 2 ) g 2 a 1 , J 1 , y = a 2 ( 1 + r 2 ) g 2 a 1 , J 1 , z = a 2 ( 1 r 2 ) g 2 a 1 , J 2 , x = a 1 ( 1 + r 2 ) g 2 a 2 , J 2 , y = a 1 ( 1 r 2 ) g 2 a 2 , J 2 , z = a 1 ( 1 r 2 ) g 2 a 2 .
In the limit r 1 , we recover the quantum compass model with the directional spin–spin couplings of Equation (19), which supports the idea of 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 a 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 transformed as S ( x ) ( S x ( x ) , ( 1 ) n 1 + 1 S y ( x ) , ( 1 ) n 1 + 1 S z ( x ) ) ,   x = ( n 1 a 1 , ( 2 n 2 1 ) a 2 ) , whereas for even rows, they transform as S ( x ) ( ( 1 ) n 1 S x ( x ) , S y ( x ) , ( 1 ) n 1 S z ( x ) ) , x = ( n 1 a 1 , ( 2 n 2 ) a 2 ) , where ( n 1 , n 2 ) Z N 1 × Z N 2 . The resulting Hamiltonian is
H eff H ^ eff = x Λ 1 j = 1 , 2 J S ( x ) · S ( x + a j e j ) + h z e i k s · x S z ( x ) ,
where the corresponding staggering wave-vector now reads k s = π a 1 e 1 + π a 2 e 2 . For a 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 the Bethe ansatz [100] and has been a subject of intense research in the past [13]. In contrast to d = 1 , all analytic and numerical evidence supports a ground state displaying a long-range, anti-ferromagnetic order in this case.
Let us now discuss the case of a nonzero Wilson parameter 0 < r < 1 , where Equation (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 being 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 for Equation (25), there is 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 a symmetric paramagnet (PM) when the transverse field is larger than the leading spin–spin coupling. 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° compass model | h z | < | J 2 , x | / ζ and | J 2 , x | > | J 1 , y | for the FM x –PM transition and the | 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, they can only be determined numerically or by using certain approximations that are discussed below. We start by deriving a path-integral representation of these effective spin models that connects to variants of the constrained QFTs discussed previously.

3. Z 2 Nonlinear Sigma Models

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 allows us to understand how constraints in QFTs, such as Equation (33), can appear in our fermionic model (8). This requires the use of spin coherent states [112], which can be defined by performing a S U ( 2 ) rotation on a fiducial state that fulfils S 2 ( x ) | S , + S x = S ( S + 1 ) | S , + S x and S z ( x ) | S , + S x = S | S , + S x . The coherent-state basis, depicted in Figure 3, can thus be defined by the action of the following operator on the tensor product of fiducial states for each lattice site
| { ω ( τ , x ) } = e x Λ d i θ ( τ , x ) sin ϕ ( τ , x ) S x ( x ) cos ϕ ( τ , x ) S y ( x ) x Λ d | S , + S x .
Here, we introduce 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 ) | 2 = 1 parametrises all possible spin directions
ω ( x ) = sin θ ( x ) cos ϕ ( x ) e x + sin θ ( x ) sin ϕ ( x ) e y + cos θ ( x ) e z .
Note that x = ( τ , x ) now represents the Wick-rotated spacetime points, where the imaginary time τ extent is related to inverse temperature via τ [ 0 , β ] [113]. One can readily check that S ( x ) = S ω ( 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 x + ω y ( x ) e y .
One can now rewrite the partition function as a path integral
Z = [ D Ω ] e S E , [ D Ω ] = x Λ d d 3 Ω ( 2 S + 1 ) 4 π δ ( Ω 2 ( τ , x ) 1 ) ,
where the vector field Ω ( τ , x ) can, in principle, take values in Ω ( τ , x ) = Ω ( τ , x ) ω ( τ , x ) R 3 but gets constrained to lie on the S 2 sphere through the path integral measure
Ω 2 ( τ , x ) = 1 , x Λ d .
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 ) on its respective sphere. These trajectories are closed due to the periodic boundary conditions along the τ direction Ω ( 0 , x ) = Ω ( β , x ) . Altogether, the action is expressed as
S E = x Λ d i Γ x d Ω · A ( Ω ( τ , x ) ) + 0 β d τ a j J j , a S 2 Ω a ( τ , x ) Ω a ( τ , x + a j e j ) + h a S Ω a ( τ , x ) ,
where the first contribution corresponds to the aforementioned Berry phase [73] 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 , which is where the fiducial state | S , + S x points to. Hence, the effective vector potential and magnetic field are those generated by the magnetic monopole
A ( Ω ( τ , x ) ) = S | Ω ( τ , x ) | e z × ω ( τ , x ) 1 + e z · ω ( τ , x ) , B ( Ω ( τ , x ) ) = Ω × A = S ω ( τ , x ) | Ω ( τ , x ) | 2 .
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 ) nonlinear 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 continuous S U ( 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 J j a j 1 2 a j ( Ω ( τ , x + a j e j ) Ω ( τ , x ) ) 2 J j a j 2 j Ω · j Ω , which clearly resembles the spatial derivative terms presented in Equation (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 Equation (1), appears 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 ( Ω ) · τ Ω , albeit still being different from the kinetic terms of the O ( 3 ) nonlinear sigma model. In a seminal work [96], F.D.M. Haldane showed that, for anti-ferromagnetic Heisenberg couplings J j = J > 0 , an expansion of Ω ( τ , x ) about the saddle point of the Euclidean action, corresponding to the alternating antiferromagnetic configuration of the classical limit, yields the kinetic term of the O ( 3 ) nonlinear sigma model exactly. Moreover, in d = 1 , the Berry phase also contributes a topological theta term θ = 2 π S which, depending on the half-integer or integer value of the spin S, makes the O ( 3 ) nonlinear sigma model massless or massive [114].
Let us now explore how this situation changes for our current spin models. Since the effective spin models do not have the continuous S U ( 2 ) symmetry of the Heisenberg limit for generic 0 < r < 1 , we first need to understand the nature of the saddle point of Equation (34) controlling the large-S limit, which will eventually lead to a different type of constrained nonlinear sigma model. By inspecting the action (34) and the magnetic monopole fields (35), one can readily see that the d Ω · A ( Ω ( τ , x ) ) term can only depend on the equatorial components of the scalar field via x Λ d ( Ω x τ Ω y Ω y τ Ω x ) . This term remains invariant under the Z 2 symmetry (15) which, in this context, reads
Ω ( τ , x ) ( Ω x ( τ , x ) , Ω y ( τ , x ) , Ω z ( τ , x ) ) .
Accordingly, the action of (34) with the constraint (32) (equivalent to the non-linear sigma constraint (33) via Equation (31)) describes a Z 2 nonlinear sigma model that arises naturally in the strong-coupling limit of the original model (8), even for a single fermion flavour N = 1 .

4. Large- S Limit and Saddle-Point Equations

Recall that the spin operators presented in Equation (12) correspond to S = 1 / 2 . In this section, we 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 will get suppressed in this limit due to the averaging of the rapid oscillations associated with the pure-imaginary Wess–Zumino term. Accordingly, the large-S limit will be controlled by static fields Ω ( τ , 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 [115] plays no role in determining the ground state in the large-N limit. Moreover, the Euclidean action can be rewritten as S E = 1 eff s E , where eff 1 / S plays the role of an effective Planck constant. In the absence of dynamic and kinetic terms, the Euclidean action per spin s E = β V eff can be expressed in terms of the following effective potential
V eff ( { Ω } ) = x Λ d a j J ˜ j , a Ω a ( x ) Ω a ( x + a j e j ) + h a Ω a ( x ) ,
where we note that, analogous to the large-N limit [7] of the four-Fermi term (5), the spin-spin couplings (22)–(27) must be rescaled to give finite contributions for S
J j , a = J ˜ j , a S ,
where J ˜ j , a are finite and have nonzero coupling strengths.
Accordingly, the large-S limit is controlled by the saddle point of the potential (40) given below, in which quantum fluctuations are suppressed 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 nonlinear constraint or work directly with the generalised ‘coordinates’ { Ω ( 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 nonlinear equations x Λ d , namely
V eff ( { θ ( x ) , ϕ ( x ) } ) θ ( x ) θ , ϕ = V eff ( { θ ( x ) , ϕ ( x ) } ) ϕ ( x ) θ , ϕ = 0 .

4.1. Large-S Ising Magnetism for d = 1

Let us start by discussing the solutions to these saddle-point equations for d = 1 . Following the discussion presented in Section 2.2, 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 nonlinear equations shown in Equation (39) by restricting them to translationally-invariant configurations within a 2-site unit cell { θ ( 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 Figure 3a). Under this simplification, the effective potential reads
V eff ( θ A , θ B , ϕ A , ϕ B ) N 1 = sin θ A sin θ B J ˜ 1 , x cos ϕ A cos ϕ B + J ˜ 1 , y sin ϕ A sin ϕ B + J ˜ 1 , z cos θ A cos θ B + 1 2 h z ( cos θ A + cos θ B ) ,
which leads to a system of four non-linear equations via (39). Let us note that this choice of angles allows for both ferromagnetic and antiferromagnetic configurations, the latter reducing the translational invariance of the model from a 1-site to a 2-site unit cell. If the angles in the A/B sites differ, one should enlarge the ansatz to rule out other possible inhomogeneous modulations. However, as discussed below in detail, we always find that the saddle-point configurations retain the full translational invariance.
To gain some knowledge about the ground state, one can numerically search for a global minimum of this potential using a coarse-grained discretisation of the angles and performing a grid search that restricts 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 nonlinear constrained problem using an interior-point algorithm of nonlinear programming [116], choosing the outcomes of the global coarser minimisation as initial points. In practice, we also initialise the minimisation for all different combinations of the cardinal states of the two S 2 spheres associated with the 2-site unit cell and check 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 starting points. This analysis yields saddle-point solutions { θ A , θ B , ϕ A , ϕ B } for different values of the Wilson parameter r { 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1 } .
In Figure 4a, we use these numerical large-S solutions to represent the corresponding SSB order parameter, which corresponds to the pseudo-scalar condensate S y ( x ) Π 5 in Equation (17). This figure clearly depicts a SSB region hosting an Ising ferromagnet FM y , corresponding to the parity-breaking Aoki phase with a nonzero 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 nonanalytically. To find the accurate location of this point, in Figure 4b, we show the corresponding susceptibilities χ y = S y ( 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, with this limit, can be found exactly
ϕ A = ϕ B π 2 , 3 π 2 , θ A = θ B = π arccos h z | J 1 , y | .
We thus find that, for | h z / J 1 , y | < | h z / J 1 , y | c = 1 , the spins align according to
S ( x ) = ± S 1 h z J 1 , y 2 e y S h z | J 1 , y | e z ,
where the two possible signs ± account for the two-fold degeneracy associated with the Z 2 parity SSB.
Now, as 0 < r < 1 is varied, 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 ground state expectation values (42), a mean-field decoupling of the additional terms of the Heisenberg–Ising chain (25) that arise when r 1
a = x , z J 1 , a S a ( x ) S a ( x + a 1 e 1 ) a = x , z 2 J 1 , a S a ( x ) S a ( x ) ,
would only contribute with terms along the internal z direction, effectively shifting the transverse field to h z h ˜ z = h z + 2 J 1 , z S z ( x ) . The saddle-point solution (42) now yields a self-consistent equation, which can be readily solved for the couplings in Equation (22):
h z J 1 , y c = 1 1 r 2 1 + r 2 .
In order to test the validity of this prediction, we present a contour plot of the Ising order parameter in Figure 5 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 Figure 4b. These are depicted with red stars, and their analytical prediction (44) is represented with a white dashed line. The level of agreement is very good, and the main source for the flow of the critical coupling with the Wilson parameter can be identified. 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 increases h ˜ 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 Equation (42). Coming back to the context of four-Fermi-Wilson lattice field theories (8), we cna 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 [117,118]. For strictly vanishing r = 0 , where the naive-fermion discretisation is recovered, the Aoki phase is no longer present. Right at the central branch, where h z = 0 , the strong-coupling ground state lies exactly at the critical line and would correspond to the gapless ground state 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 the main sources for the flow of the critical point qualitatively. In fact, for the quantum Ising model at r = 1 , where a large-S predicts a critical point | h z / J 1 , z | = 1 in Figure 4a,b, the exact solution gives instead instead | h z / J 1 , z | = 1 / 2 [85]. In Section 5, we present more accurate predictions of the critical points using the quasi-exact DMRG algorithm based on matrix product states.

4.2. Large-S Compass Magnetism for d = 2

Let us now discuss the d = 2 case, where the A and B sublattices correspond to two interpenetrating square lattices with the lattice spacing ( a 1 2 + a 2 2 ) 1 / 2 , rotated with respect to the original rectangular lattice by angles φ , π φ , where φ = arctan ( a 2 / a 1 ) (see Figure 3b). Analogous to d = 1 , we can restrict the configurations of spin coherent states to be translationally invariant within these sublattices { θ ( x ) , ϕ ( x ) } { θ s , ϕ s } , where s { A , B } , and we can still account for antiferromagnetic configurations when the respective A , B angles differ or ferromagnetic ones when they are equal. However, this choice may not suffice to capture the ground state 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 Equation (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 angles θ A * = θ B * = π / 2 , ϕ A * = ϕ B * { 0 , π } . Additionally, since the perturbative coupling between neighbouring columns fulfils J 1 , x > 0 , spins in adjacent columns might minimise the ground state energy by choosing opposing azimuthal angles ϕ A * = ϕ B * { π , 0 } . Since this is inconsistent with the A , B sublattice layout, allowing for this possible ordering in the parametrisation of the constrained effective potential requires the number of configurations to be augmented by considering a 4-site unit cell { θ ( x ) , ϕ ( x ) } { θ s , c , ϕ s , c } , where c { l , r } labels the left and right corners, as depicted in Figure 3b. By directly incorporating the nonlinear constraint, as we did for Equation (40), the effective potential reads
2 V eff N 1 N 2 = c ( j sin θ A , c sin θ B , c ˜ ( j ) J ˜ j , x cos ϕ A , c cos ϕ B , c ˜ ( j ) + J ˜ j , y sin ϕ A , c sin ϕ B , c ˜ ( j ) + J ˜ j , z cos θ A , c cos θ B , c ˜ ( j ) + h z 2 ( cos θ A , c + cos θ B , c ) ) ,
where the function c ˜ ( 2 ) = c is introduced when the spin–spin couplings occur along the e 2 spatial direction and c ˜ ( 1 ) = c ¯ along 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 nonlinear 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° compass model [88] in a transverse field, and one can find
ϕ s , c π 2 , 3 π 2 , θ s , c = π arccos h z | J 1 , y | , if a 1 < a 2 , ϕ s , c 0 , π , θ s , c = π arccos h z | J 2 , x | , if a 1 > a 2 .
The SSB order parameters for these solutions are
S ( x ) = ± S 1 h z J 1 , y 2 e y S h z | J 1 , y | e z , if a 1 < a 2 , S ( x ) = ± S 1 h z J 2 , x 2 e x S h z | J 2 , x | e z , if a 1 > a 2 ,
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 ) Π 1 corresponds to ferromagnetic ordering along the internal x axis FM x , which corresponds to the inversion-breaking π condensate of Equation (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 ) Π 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 nonlinear 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 nonanalytic behaviours, which can be used to infer the locations of the critical points and how these flow as one varies r. In Figure 6, we present a stack of two-dimensional contour plots that summarises the large-S phase diagram and shows a clear dependence on both the Wilson parameter and the anisotropy parameter
ξ 2 = a 1 a 2 .
In this contour plot, we represent the difference between the two possible SSB order parameters S y ( x ) S x ( x ) , such that negative (positive) values signal a FM x (FM y ) phase with a Π 1 ( Π 2 ) Lorentz-breaking condensate and are depicted on a 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, it can be observed 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 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 Equation (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 ). By setting r 1 , one can solve the corresponding self-consistent equations approximately, leading to
h z J 1 , y c 1 1 r 2 + ( 1 r 2 ) ξ 2 2 1 + r 2 1 ( 1 r 2 ) ξ 2 2 1 + r 2 , if a 1 < a 2 , h z J 2 , x c 1 1 r 2 + ( 1 r 2 ) ξ 2 2 ( 1 + r 2 ) 1 ( 1 r 2 ) ξ 2 2 1 + r 2 , if a 1 > a 2 .
The white dashed lines of Figure 6 represent the analytical predictions of the critical points. As can be observed, they match the numerical critical lines well for large anisotropies ξ 2 1 ( ξ 2 1 ), but the discrepancy increases as one approaches the isotropic case ξ 2 = 1 .

5. Tensor-Network Numerical Simulations

In this section, we test the validity of the large-S predictions of Section 2 by means of variational algorithms based on tensor network states (TNSs) [119,120,121].
The quantum state of a lattice model composed of N s d ˜ -level systems, i.e., spins, can be written on the basis of the 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 dimension d ˜
| ψ = i 1 , i 2 , , i N s C i 1 , i 2 , , i N s | i 1 , i 2 , , i N s .
Here, we introduce the indexes i n { 1 , , d ˜ } , such that the description of the state requires d ˜ 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 have a more concise description based on TNSs [122,123]. Obtained from a contraction of low-rank tensors on so-called virtual indices, TNSs economically approximate the states of a system with local interactions in thermal equilibrium. The number of required parameters scales only polynomially with system size [124], circumventing the previous exponential growth of the most generic description. In fact, these variational states are based on powerful insights related to the area law [125,126]. The area law places bounds on the 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 have benefited from the advent of White’s density matrix renormalisation group (DMRG) [127], famous for its extraordinary accuracy in solving one-dimensional quantum systems, which is intimately connected with a tensor decomposition known as the matrix product state (MPS) [128,129]. 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 through 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 = α = 1 D | α , α , often referred to as an entangled bond. By applying the map that plays the role of the aforementioned projector
A = A i , α β | i α β |
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 Figure 7a). In the d = 1 -dimensional models discussed in this work, we consider physical spins S = 1 / 2 , such that i n = s x n { , } , d ˜ = 2 , and n 1 , , N s , with N s = N 1 being the number of sites in the spatial chain. In general, the dimension of the entangled state | I can be site-dependent and A s x n [ x n ] is written for the D n × D n + 1 matrix corresponding to site x n ; the states then have the form
| ψ = { s x n } Tr A s x 1 [ s x 1 ] A s x 2 [ s x 2 ] A s x N s [ s x N s ] | s x 1 , s x 2 , , s x N s
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 Figure 7b. In this case, the number of parameters needed to describe a physical state in the MPS language scales as O ( N 1 d D 2 ) with d being the physical dimension of the spins.
A natural generalisation of MPS to two, or even more, spatial dimensions is represented by projected entangled pair states (PEPS [130]). 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 Figure 7c, 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 [130,131,132] as
| ψ = s x n F A s x 1 [ x 1 ] , A s x 2 [ x 2 ] , , A s x N 1 N 2 [ x N 1 N 2 ] | s x 1 , s x 2 , , s x N 1 N 2 .
This PEPS is represented by a network of N 1 N 2 tensors A s x n [ x m ] , 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 d D 4 complex coefficients. Therefore, the PEPS is characterized by O ( N s d D 4 ) parameters. The function F contracts all tensors A s x m [ x m ] according to this pattern and then performs a trace to obtain a scalar quantity, such that Equation (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 ground state 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. In contrast, for d = 2 , analogous 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 ground state is reached. This assumes that this ground state is unique and that the energy gap is nonzero, as done with the time-evolving block-decimation method (TEBD) for one-dimensional chains [133,134]. In the following text, we use this method with the thermodynamic limit for the infinite PEPS state (iPEPS) [135,136].

5.1. 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 the spin–spin couplings defined in Equation (22) following subjection to an additional transverse field in Equation (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 Section 4.1, giving more accurate predictions of the phase diagram and critical points presented in Figure 5.
As discussed in Section 4.1, the Heisenberg–Ising chain presents a critical line separating the ferromagnet FM y and the paramagnet PM. In Figure 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 the bond dimension D = 200 (red stars) or by the self-consistent mean-field method (yellow dashed lines), which exploits exact solutions of the transverse-field Ising model to derive a self-consistent equation for the transverse magnetisation that can be solved analytically in the limit | J 1 , x / J 1 , y | = | J 1 , z / J 1 , y | 1 . This yields
h z J 1 , y c = 1 2 4 π 1 r 2 1 + r 2 ,
which must be compared to the large-S estimate (54), depicted by a white dashed line in the figure.
As can be observed in Figure 8, for the 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. By 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 in the critical lines, e.g., for r = 1 h z / J 1 , y c = 1 within a 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 (pseudo-scalar condensate) will shrink 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 Figure 8, we present a contour plot of the SSB order parameter, which corresponds to the pseudo-scalar condensate S y ( x ) Π 5 defined in Equation (17). In order to avoid numerical problems due to the incomplete symmetry breaking of the magnetisation M y = S y ( x ) = 1 N 1 x Λ 1 S y ( x ) , we determine instead the corresponding structure factors
S y y ( k ) = 1 N 1 2 n 1 , n 1 e i k a 1 ( n 1 n 1 ) S y ( n 1 ) S y ( n 1 ) .
The zero-momentum component of these structure factors yields 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 on a yellow-green scale with a nonzero pseudo-scalar condensate, namely a Ising ferromagnet FM y , which is separated from the region by a blue scale, where the parity is preserved, namely through a paramagnet PM with zero magnetisation along the internal y axis.
Let us now give some more details on the methodology used to numerically extract the critical points shown in Figure 8. By calculating the ferromagnetic magnetisations and, particularly, the corresponding susceptibilities, we can identify the critical points occurring at a nonzero external field h z > 0 . In Figure 9a, we present the susceptibility χ y = M y / h z for different values of the Wilson parameter r, which clearly peaks at specific values of the ratio h z / J 1 , y , which moves to the left as r is decreased. In Figure 9b, we display the ferromagnetic susceptibility χ y for different numbers of sites N s , fixing r = 0.82 and vary 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 the transverse field h z c diverges with the size of the chain, and by fitting the maxima of h z c to h z c ( N 1 ) = h z c ( 1 + a N 1 1 + b N 1 2 ) , 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 magnetisation 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 perturbations to the Ising limit in the form of the additional spin–spin couplings (22) when r < 1 .

5.2. Tensor-Network Compass Magnetism for d = 2

In this section, we show the results obtained by using the above iPEPS algorithm for the Heisenberg–Ising compass model (13) with the spin–spin couplings defined in Equation (27) and subject to an additional transverse field in Equation (26) that works directly in the infinite-lattice limit. In particular, we compute 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 evaluate observable quantities on it, such as the ground state energy and the local order parameters related to the ferromagnetic phases.
In Section 4.2, 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 found analytically (49), corresponding to the white dashed lines of Figure 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 for a nonzero magnetic field h z > 0 .
We start by setting the spatial anisotropy to ξ 2 = 3.16 , which corresponds to a specific stacked plane in the large-S phase diagram of Figure 6, where the FM x order competes with the disordered PM. Our numerical iPEPS results for this competition are presented in Figure 10. The lines correspond to the critical points where the FM x and PM phase transitions occur, either obtained through an imaginary time evolution based on infinite iPEPS with the bond dimension D = 2 (red stars) or by large-S approximate prediction (49) (orange dashed lines). We observe a similar trend as for 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 nonvanishing 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 with d = 1 (see the red stars in Figure 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 with the d = 1 case, showcasing that mean-field predictions typically improve as the dimensionality increases.
Let us again discuss the details on how we extracted these critical points numerically. In the inset of Figure 11, we present the magnetisation M x = S x ( x ) as a function of a 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 magnetisation attains a nonzero 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 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 the bond dimension D > 2 , which would permit 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 a long-range order in contrast with the 1D version. We note that these questions could also be addressed with other methods, such as the Monte Carlo simulation.

6. Conclusions and Discussion

In this work, we explored the limit of strong Hubbard interactions in models of correlated topological insulators that arise for spin-orbit coupled fermions in lattices with 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 discretisation, demonstrating 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 discretisation to r = 1 , as a nonunity value has trivial consequences in the noninteracting limit. However, the role of 0 < r < 1 in the presence of interactions is not clear a priori, and this was the topic 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 the 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 derived an effective spin model for the limit of strong four-Fermi interactions and found a specific dependence of the couplings on the Wilson parameter r. For 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 coupling. We formulated a path-integral representation of the partition function, which connects the strongly-interacting limit with a constrained QFT: a nonlinear sigma model with a discrete Z 2 symmetry. This permitted exact solutions to be found for a large-S limit, which enabled 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 was tested against tensor-network numerical simulations, which showed that the large-S diagrams are qualitatively correct. On the other hand, the numerical results gave more accurate estimates of the flow of the critical lines and, in some cases, allowed us to infer the correct universality class of the lines that differs from the large-S mean-field-type scaling.
As an outlook, we believe that the present manuscript, together with other recent works [57,58,59,60,61], provides a rich cross-disciplinary toolbox that can be used 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 the abandonment of half-filling and exploration of correlated topological phases at nonzero fermion densities. This would connect to recent works on the existence of inhomogeneous kink crystals in the Gross–Neveu model with discrete chiral symmetry [137,138], extend it to ( 2 + 1 ) -dimensional spacetimes, and explore whether such inhomogeneities can also appear in the parity- and Lorentz-breaking condensates that are characteristic of Wilson-type discretisations.

Author Contributions

Conceptualization, S.H. and A.B.; methodology, E.T. and A.B.; software, E.T.; validation, E.T., A.B.; formal analysis, E.T.; investigation, E.T.; resources, E.T. and A.B.; data curation, E.T. and A.B.; writing—review and editing, E.T., S.H., A.B.; visualization, E.T. and A.B.; supervision, A.B.; project administration, A.B.; funding acquisition, A.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

A.B. acknowledges support from PGC2018-099169-B-I00 (MCIU/AEI/FEDER, UE) through a grant from the IFT Centro de Excelencia Severo Ochoa CEX2020-001007-S, funded by MCIN/AEI/ 10.13039/501100011033, a grant from QUITEMAD+S2013/ICE-2801 and a grant from the CSIC Research Platform on Quantum Technologies PTI-001. S.H. is supported by an STFC grant ST/T000813/1.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Peskin, M.E.; Schroeder, D.V. An Introduction to Quantum Field Theory; Addison-Wesley: Reading, PA, USA, 1995. [Google Scholar]
  2. Greensite, J. Introduction to the Confinement Problem; Springer Nature: Berlin/Heidelberg, Germany, 2020. [Google Scholar]
  3. Thirring, W.E. A soluble relativistic field theory. Ann. Phys. 1958, 3, 91–112. [Google Scholar] [CrossRef]
  4. Gross, D.J.; Neveu, A. Dynamical symmetry breaking in asymptotically free field theories. Phys. Rev. D 1974, 10, 3235–3253. [Google Scholar] [CrossRef]
  5. Polyakov, A. Interaction of goldstone particles in two dimensions. Applications to ferromagnets and massive Yang-Mills fields. Phys. Lett. B 1975, 59, 79–81. [Google Scholar] [CrossRef]
  6. Coleman, S. Quantum sine-Gordon equation as the massive Thirring model. Phys. Rev. D 1975, 11, 2088–2097. [Google Scholar] [CrossRef] [Green Version]
  7. Coleman, S. Aspects of Symmetry: Selected Erice Lectures; Cambridge University Press: Cambridge, UK, 1985. [Google Scholar] [CrossRef]
  8. Wilson, K.G.; Kogut, J. The renormalization group and the ϵ expansion. Phys. Rep. 1974, 12, 75–199. [Google Scholar] [CrossRef]
  9. Anderson, P.W. More Is Different. Science 1972, 177, 393–396. [Google Scholar] [CrossRef] [Green Version]
  10. Castro Neto, A.H.; Guinea, F.; Peres, N.M.R.; Novoselov, K.S.; Geim, A.K. The electronic properties of graphene. Rev. Mod. Phys. 2009, 81, 109–162. [Google Scholar] [CrossRef] [Green Version]
  11. Armitage, N.P.; Mele, E.J.; Vishwanath, A. Weyl and Dirac semimetals in three-dimensional solids. Rev. Mod. Phys. 2018, 90, 015001. [Google Scholar] [CrossRef] [Green Version]
  12. Qi, X.L.; Zhang, S.C. Topological insulators and superconductors. Rev. Mod. Phys. 2011, 83, 1057–1110. [Google Scholar] [CrossRef] [Green Version]
  13. Manousakis, E. The spin-½ Heisenberg antiferromagnet on a square lattice and its application to the cuprous oxides. Rev. Mod. Phys. 1991, 63, 1. [Google Scholar] [CrossRef]
  14. Feynman, R.P. Simulating physics with computers. Int. J. Theor. Phys. 1982, 21, 467. [Google Scholar] [CrossRef]
  15. Cirac, J.I.; Zoller, P. Goals and opportunities in quantum simulation. Nat. Phys. 2012, 8, 264–266. [Google Scholar] [CrossRef]
  16. Wiese, U.J. Ultracold quantum gases and lattice systems: Quantum simulation of lattice gauge theories. Ann. Phys. 2013, 525, 777–796. [Google Scholar] [CrossRef] [Green Version]
  17. Zohar, E.; Cirac, J.I.; Reznik, B. Quantum simulations of lattice gauge theories using ultracold atoms in optical lattices. Rep. Prog. Phys. 2015, 79, 014401. [Google Scholar] [CrossRef] [Green Version]
  18. Dalmonte, M.; Montangero, S. Lattice gauge theory simulations in the quantum information era. Contemp. Phys. 2016, 57, 388–412. [Google Scholar] [CrossRef]
  19. Bañuls, M.C.; Blatt, R.; Catani, J.; Celi, A.; Cirac, J.I.; Dalmonte, M.; Fallani, L.; Jansen, K.; Lewenstein, M.; Montangero, S.; et al. Simulating lattice gauge theories within quantum technologies. Eur. Phys. J. D 2020, 74, 165. [Google Scholar] [CrossRef]
  20. Aidelsburger, M.; Barbiero, L.; Bermudez, A.; Chanda, T.; Dauphin, A.; González-Cuadra, D.; Grzybowski, P.R.; Hands, S.; Jendrzejewski, F.; Jünemann, J.; et al. Cold atoms meet lattice gauge theory. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 2022, 380, 20210064. [Google Scholar] [CrossRef]
  21. Klco, N.; Roggero, A.; Savage, M.J. Standard Model Physics and the Digital Quantum Revolution: Thoughts about the Interface. arXiv 2021, arXiv:2107.04769. [Google Scholar] [CrossRef]
  22. Kogut, J.B. An introduction to lattice gauge theory and spin systems. Rev. Mod. Phys. 1979, 51, 659–713. [Google Scholar] [CrossRef]
  23. Kogut, J.; Susskind, L. Hamiltonian formulation of Wilson’s lattice gauge theories. Phys. Rev. D 1975, 11, 395–408. [Google Scholar] [CrossRef]
  24. Gerritsma, R.; Kirchmair, G.; Zähringer, F.; Solano, E.; Blatt, R.; Roos, C.F. Quantum simulation of the Dirac equation. Nature 2010, 463, 68–71. [Google Scholar] [CrossRef] [Green Version]
  25. Gerritsma, R.; Lanyon, B.P.; Kirchmair, G.; Zähringer, F.; Hempel, C.; Casanova, J.; García-Ripoll, J.J.; Solano, E.; Blatt, R.; Roos, C.F. Quantum Simulation of the Klein Paradox with Trapped Ions. Phys. Rev. Lett. 2011, 106, 060503. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Tarruell, L.; Greif, D.; Uehlinger, T.; Jotzu, G.; Esslinger, T. Creating, moving and merging Dirac points with a Fermi gas in a tunable honeycomb lattice. Nature 2012, 483, 302–305. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Duca, L.; Li, T.; Reitter, M.; Bloch, I.; Schleier-Smith, M.; Schneider, U. An Aharonov-Bohm interferometer for determining Bloch band topology. Science 2015, 347, 288–292. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Jotzu, G.; Messer, M.; Desbuquois, R.; Lebrat, M.; Uehlinger, T.; Greif, D.; Esslinger, T. Experimental realization of the topological Haldane model with ultracold fermions. Nature 2014, 515, 237–240. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Kang, J.H.; Han, J.H.; Shin, Y. Realization of a Cross-Linked Chiral Ladder with Neutral Fermions in a 1D Optical Lattice by Orbital-Momentum Coupling. Phys. Rev. Lett. 2018, 121, 150403. [Google Scholar] [CrossRef] [Green Version]
  30. Liang, M.C.; Wei, Y.D.; Zhang, L.; Wang, X.J.; Zhang, H.; Wang, W.W.; Qi, W.; Liu, X.J.; Zhang, X. Realization of Qi-Wu-Zhang model in spin-orbit-coupled ultracold fermions. arXiv 2021, arXiv:2109.08885. [Google Scholar]
  31. Martinez, E.A.; Muschik, C.A.; Schindler, P.; Nigg, D.; Erhard, A.; Heyl, M.; Hauke, P.; Dalmonte, M.; Monz, T.; Zoller, P.; et al. Real-time dynamics of lattice gauge theories with a few-qubit quantum computer. Nature 2016, 534, 516–519. [Google Scholar] [CrossRef] [Green Version]
  32. Schweizer, C.; Grusdt, F.; Berngruber, M.; Barbiero, L.; Demler, E.; Goldman, N.; Bloch, I.; Aidelsburger, M. Floquet approach to Z2 lattice gauge theories with ultracold atoms in optical lattices. Nat. Phys. 2019, 15, 1168–1173. [Google Scholar] [CrossRef] [Green Version]
  33. Kokail, C.; Maier, C.; van Bijnen, R.; Brydges, T.; Joshi, M.K.; Jurcevic, P.; Muschik, C.A.; Silvi, P.; Blatt, R.; Roos, C.F.; et al. Self-verifying variational quantum simulation of lattice models. Nature 2019, 569, 355–360. [Google Scholar] [CrossRef]
  34. Mil, A.; Zache, T.V.; Hegde, A.; Xia, A.; Bhatt, R.P.; Oberthaler, M.K.; Hauke, P.; Berges, J.; Jendrzejewski, F. A scalable realization of local U(1) gauge invariance in cold atomic mixtures. Science 2020, 367, 1128–1130. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Surace, F.M.; Mazza, P.P.; Giudici, G.; Lerose, A.; Gambassi, A.; Dalmonte, M. Lattice Gauge Theories and String Dynamics in Rydberg Atom Quantum Simulators. Phys. Rev. X 2020, 10, 021041. [Google Scholar] [CrossRef]
  36. Yang, B.; Sun, H.; Ott, R.; Wang, H.Y.; Zache, T.V.; Halimeh, J.C.; Yuan, Z.S.; Hauke, P.; Pan, J.W. Observation of gauge invariance in a 71-site Bose–Hubbard quantum simulator. Nature 2020, 587, 392–396. [Google Scholar] [CrossRef] [PubMed]
  37. Fermi, E. Versuch einer Theorie der β-Strahlen. I. Z. Phys. 1934, 88, 161–177. [Google Scholar] [CrossRef]
  38. Wilson, F.L. Fermi’s Theory of Beta Decay. Am. J. Phys. 1968, 36, 1150–1160. [Google Scholar] [CrossRef]
  39. Nambu, Y.; Jona-Lasinio, G. Dynamical Model of Elementary Particles Based on an Analogy with Superconductivity. I. Phys. Rev. 1961, 122, 345–358. [Google Scholar] [CrossRef] [Green Version]
  40. Nambu, Y.; Jona-Lasinio, G. Dynamical Model of Elementary Particles Based on an Analogy with Superconductivity. II. Phys. Rev. 1961, 124, 246–254. [Google Scholar] [CrossRef]
  41. Wilson, K.G. Quarks and Strings on a Lattice. In New Phenomena in Subnuclear Physics; Springer: Boston, MA, USA, 1977; pp. 69–142. [Google Scholar] [CrossRef]
  42. Qi, X.L.; Wu, Y.S.; Zhang, S.C. Topological quantization of the spin Hall effect in two-dimensional paramagnetic semiconductors. Phys. Rev. B 2006, 74, 085308. [Google Scholar] [CrossRef] [Green Version]
  43. Qi, X.L.; Hughes, T.L.; Zhang, S.C. Topological field theory of time-reversal invariant insulators. Phys. Rev. B 2008, 78, 195424. [Google Scholar] [CrossRef] [Green Version]
  44. Ryu, S.; Schnyder, A.P.; Furusaki, A.; Ludwig, A.W.W. Topological insulators and superconductors: Tenfold way and dimensional hierarchy. New J. Phys. 2010, 12, 065010. [Google Scholar] [CrossRef]
  45. Bermudez, A.; Mazza, L.; Rizzi, M.; Goldman, N.; Lewenstein, M.; Martin-Delgado, M.A. Wilson Fermions and Axion Electrodynamics in Optical Lattices. Phys. Rev. Lett. 2010, 105, 190404. [Google Scholar] [CrossRef]
  46. Mazza, L.; Bermudez, A.; Goldman, N.; Rizzi, M.; Martin-Delgado, M.A.; Lewenstein, M. An optical-lattice-based quantum simulator for relativistic field theories and topological insulators. New J. Phys. 2012, 14, 015007. [Google Scholar] [CrossRef] [Green Version]
  47. Kaplan, D.B.; Sun, S. Spacetime as a Topological Insulator: Mechanism for the Origin of the Fermion Generations. Phys. Rev. Lett. 2012, 108, 181807. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Zache, T.V.; Hebenstreit, F.; Jendrzejewski, F.; Oberthaler, M.K.; Berges, J.; Hauke, P. Quantum simulation of lattice gauge theories using Wilson fermions. Quantum Sci. Technol. 2018, 3, 034010. [Google Scholar] [CrossRef] [Green Version]
  49. Hands, S.; Kocic, A.; Kogut, J. Four-Fermi Theories in Fewer Than Four Dimensions. Ann. Phys. 1993, 224, 29–89. [Google Scholar] [CrossRef] [Green Version]
  50. Hands, S. Fixed Point Four-Fermi Theories. arXiv 1997, arXiv:9706018. [Google Scholar]
  51. Hohenadler, M.; Assaad, F.F. Correlation effects in two-dimensional topological insulators. J. Phys. Condens. Matter 2013, 25, 143201. [Google Scholar] [CrossRef]
  52. Rachel, S. Interacting topological insulators: A review. Rep. Prog. Phys. 2018, 81, 116501. [Google Scholar] [CrossRef] [Green Version]
  53. Neupert, T.; Chamon, C.; Iadecola, T.; Santos, L.H.; Mudry, C. Fractional (Chern and topological) insulators. Phys. Scr. 2015, 2015, 014005. [Google Scholar] [CrossRef] [Green Version]
  54. Bergholtz, E.J.; Liu, Z. Topological flat band models and fractional chern insulators. Int. J. Mod. Phys. B 2013, 27, 1330017. [Google Scholar] [CrossRef] [Green Version]
  55. Jaksch, D.; Bruder, C.; Cirac, J.I.; Gardiner, C.W.; Zoller, P. Cold Bosonic Atoms in Optical Lattices. Phys. Rev. Lett. 1998, 81, 3108–3111. [Google Scholar] [CrossRef]
  56. Hofstetter, W.; Cirac, J.I.; Zoller, P.; Demler, E.; Lukin, M.D. High-Temperature Superfluidity of Fermionic Atoms in Optical Lattices. Phys. Rev. Lett. 2002, 89, 220407. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Jünemann, J.; Piga, A.; Ran, S.J.; Lewenstein, M.; Rizzi, M.; Bermudez, A. Exploring Interacting Topological Insulators with Ultracold Atoms: The Synthetic Creutz-Hubbard Model. Phys. Rev. X 2017, 7, 031057. [Google Scholar] [CrossRef]
  58. Bermudez, A.; Tirrito, E.; Rizzi, M.; Lewenstein, M.; Hands, S. Gross-Neveu-Wilson model and correlated symmetry-protected topological phases. Ann. Phys. 2018, 399, 149–180. [Google Scholar] [CrossRef]
  59. Tirrito, E.; Lewenstein, M.; Bermudez, A. Topological chiral currents in the Gross-Neveu model extension. arXiv 2021, arXiv:2112.07654. [Google Scholar]
  60. Ziegler, L.; Tirrito, E.; Lewenstein, M.; Hands, S.; Bermudez, A. Correlated Chern insulators in two-dimensional Raman lattices: A cold-atom regularization of strongly-coupled four-Fermi field theories. arXiv 2020, arXiv:2011.08744. [Google Scholar]
  61. Ziegler, L.; Tirrito, E.; Lewenstein, M.; Hands, S.; Bermudez, A. Large-N Chern insulators: Lattice field theory and quantum simulation approaches to correlation effects in the quantum anomalous Hall effect. Ann. Phys. 2022, 439, 168763. [Google Scholar] [CrossRef]
  62. Gell-Mann, M.; Lévy, M. The axial vector current in beta decay. Il Nuovo Cimento (1955–1965) 1960, 16, 705–726. [Google Scholar] [CrossRef]
  63. Coleman, S. There are no Goldstone bosons in two dimensions. Commun. Math. Phys. 1973, 31, 259–264. [Google Scholar] [CrossRef]
  64. Polyakov, A.M.; Belavin, A. Metastable States of Two-Dimensional Isotropic Ferromagnets. Jetp Lett. 1975, 22, 503–506. [Google Scholar]
  65. Bardeen, W.A.; Lee, B.W.; Shrock, R.E. Phase transition in the nonlinear σ model in a (2+ϵ)-dimensional continuum. Phys. Rev. D 1976, 14, 985–1005. [Google Scholar] [CrossRef]
  66. Gattringer, C.; Lang, C.B. Quantum Chromodynamics on the Lattice: An Introductory Presentation; Springer: Berlin/Heidelberg, Germany, 2009. [Google Scholar]
  67. Nielsen, H.; Ninomiya, M. Absence of neutrinos on a lattice: (I). Proof by homotopy theory. Nucl. Phys. B 1981, 185, 20–40. [Google Scholar] [CrossRef]
  68. Nielsen, H.; Ninomiya, M. Absence of neutrinos on a lattice: (II). Intuitive topological proof. Nucl. Phys. B 1981, 193, 173–194. [Google Scholar] [CrossRef]
  69. Montvay, I.; Münster, G. Quantum Fields on a Lattice; Cambridge Monographs on Mathematical Physics; Cambridge University Press: Cambridge, UK, 1994. [Google Scholar]
  70. Amato, A.; Aarts, G.; Allton, C.; Giudice, P.; Hands, S.; Skullerud, J.I. Electrical Conductivity of the Quark-Gluon Plasma Across the Deconfinement Transition. Phys. Rev. Lett. 2013, 111, 172001. [Google Scholar] [CrossRef] [Green Version]
  71. China Lattice QCD Collaboration (CLQCD); Li, X.; Chen, Y.; Meng, G.Z.; Feng, X.; Gong, M.; He, S.; Li, G.; Liu, C.; Liu, Y.B.; et al. Hadron scattering in an asymmetric. J. High Energy Phys. 2007, 2007, 53. [Google Scholar] [CrossRef] [Green Version]
  72. Berry, M.V. Quantal phase factors accompanying adiabatic changes. Proc. R. Soc. Lond. A Math. Phys. Sci. 1984, 392, 45–57. [Google Scholar] [CrossRef]
  73. Xiao, D.; Chang, M.C.; Niu, Q. Berry phase effects on electronic properties. Rev. Mod. Phys. 2010, 82, 1959–2007. [Google Scholar] [CrossRef] [Green Version]
  74. Nakahara, M. Geometry, Topology and Physics; CRC Press: Boca Raton, FL, USA, 2017. [Google Scholar]
  75. Zak, J. Berry’s phase for energy bands in solids. Phys. Rev. Lett. 1989, 62, 2747–2750. [Google Scholar] [CrossRef]
  76. Schnyder, A.P.; Ryu, S.; Furusaki, A.; Ludwig, A.W.W. Classification of topological insulators and superconductors in three spatial dimensions. Phys. Rev. B 2008, 78, 195125. [Google Scholar] [CrossRef] [Green Version]
  77. Kitaev, A. Periodic table for topological insulators and superconductors. AIP Conf. Proc. 2009, 1134, 22–30. [Google Scholar] [CrossRef] [Green Version]
  78. Liu, C.X.; Zhang, S.C.; Qi, X.L. The Quantum Anomalous Hall Effect: Theory and Experiment. Annu. Rev. Condens. Matter Phys. 2016, 7, 301–321. [Google Scholar] [CrossRef]
  79. Kaplan, D.B. A method for simulating chiral fermions on the lattice. Phys. Lett. B 1992, 288, 342–347. [Google Scholar] [CrossRef] [Green Version]
  80. Kaplan, D.B.; Sen, S. Index theorems, generalized Hall currents and topology for gapless defect fermions. arXiv 2021, arXiv:2112.06954. [Google Scholar]
  81. Golterman, M.F.; Jansen, K.; Kaplan, D.B. Chern-Simons currents and chiral fermions on the lattice. Phys. Lett. B 1993, 301, 219–223. [Google Scholar] [CrossRef] [Green Version]
  82. Sen, S. Chern insulator transitions with Wilson fermions on a hyperrectangular lattice. Phys. Rev. D 2020, 102, 094520. [Google Scholar] [CrossRef]
  83. Anderson, P.W. Antiferromagnetism. Theory of Superexchange Interaction. Phys. Rev. 1950, 79, 350–356. [Google Scholar] [CrossRef]
  84. Anderson, P.W. Theory of Magnetic Exchange Interactions:Exchange in Insulators and Semiconductors; Solid State Physics; Academic Press: Cambridge, MA, USA, 1963; Volume 14, pp. 99–214. [Google Scholar] [CrossRef]
  85. Pfeuty, P. The one-dimensional Ising model with a transverse field. Ann. Phys. 1970, 57, 79–90. [Google Scholar] [CrossRef]
  86. Aoki, S. New phase structure for lattice QCD with Wilson fermions. Phys. Rev. D 1984, 30, 2653–2663. [Google Scholar] [CrossRef]
  87. Sharpe, S.; Singleton, R. Spontaneous flavor and parity breaking with Wilson fermions. Phys. Rev. D 1998, 58, 074501. [Google Scholar] [CrossRef] [Green Version]
  88. Nussinov, Z.; Fradkin, E. Discrete sliding symmetries, dualities, and self-dualities of quantum orbital compass models and p + ip superconducting arrays. Phys. Rev. B 2005, 71, 195120. [Google Scholar] [CrossRef] [Green Version]
  89. Nussinov, Z.; van den Brink, J. Compass models: Theory and physical motivations. Rev. Mod. Phys. 2015, 87, 1. [Google Scholar] [CrossRef] [Green Version]
  90. Dorier, J.; Becca, F.; Mila, F. Quantum compass model on the square lattice. Phys. Rev. B 2005, 72, 024448. [Google Scholar] [CrossRef] [Green Version]
  91. Chen, H.D.; Fang, C.; Hu, J.; Yao, H. Quantum phase transition in the quantum compass model. Phys. Rev. B 2007, 75, 144401. [Google Scholar] [CrossRef] [Green Version]
  92. Orús, R.; Doherty, A.C.; Vidal, G. First order phase transition in the anisotropic quantum orbital compass model. Phys. Rev. Lett. 2009, 102, 077203. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  93. Heisenberg, W. Zur Theorie des Ferromagnetismus. Z. Phys. 1928, 49, 619–636. [Google Scholar] [CrossRef]
  94. Haldane, F. Continuum dynamics of the 1-D Heisenberg antiferromagnet: Identification with the O(3) nonlinear sigma model. Phys. Lett. A 1983, 93, 464–468. [Google Scholar] [CrossRef]
  95. Haldane, F.D.M. Nonlinear Field Theory of Large-Spin Heisenberg Antiferromagnets: Semiclassically Quantized Solitons of the One-Dimensional Easy-Axis Néel State. Phys. Rev. Lett. 1983, 50, 1153–1156. [Google Scholar] [CrossRef]
  96. Haldane, F.D.M. O(3) Nonlinear σ Model and the Topological Distinction between Integer- and Half-Integer-Spin Antiferromagnets in Two Dimensions. Phys. Rev. Lett. 1988, 61, 1029–1032. [Google Scholar] [CrossRef]
  97. Affleck, I. Field Theory Methods and Strongly Correlated Electrons. In Physics, Geometry and Topology; Lee, H.C., Ed.; Springer: Boston, MA, USA, 1990; pp. 1–13. [Google Scholar] [CrossRef]
  98. Susskind, L. Lattice fermions. Phys. Rev. D 1977, 16, 3031–3039. [Google Scholar] [CrossRef]
  99. Hands, S. Critical flavor number in the 2+1D Thirring model. Phys. Rev. D 2019, 99, 034504. [Google Scholar] [CrossRef] [Green Version]
  100. Bethe, H. Zur Theorie der Metalle. Z. Phys. 1931, 71, 205–226. [Google Scholar] [CrossRef]
  101. Hulthén, L. Über das Austauschproblem Eines Kristalles. Ph.D. Thesis, Almqvist & Wiksell, Stockholm, Sweden, 1938. [Google Scholar]
  102. Bogoliubov, N.; Izergin, A.; Korepin, V. Critical exponents for integrable models. Nucl. Phys. B 1986, 275, 687–705. [Google Scholar] [CrossRef]
  103. Kasteleijn, P. The lowest energy state of a linear antiferromagnetic chain. Physica 1952, 18, 104–113. [Google Scholar] [CrossRef]
  104. Luther, A.; Peschel, I. Calculation of critical exponents in two dimensions from quantum field theory in one dimension. Phys. Rev. B 1975, 12, 3908–3917. [Google Scholar] [CrossRef]
  105. Berezinsky, V.L. Destruction of Long-range Order in One-dimensional and Two-dimensional Systems Possessing a Continuous Symmetry Group. II. Quantum Systems. Sov. Phys. JETP 1972, 34, 610. [Google Scholar]
  106. Kosterlitz, J.M.; Thouless, D.J. Ordering, metastability and phase transitions in two-dimensional systems. J. Phys. C Solid State Phys. 1973, 6, 1181–1203. [Google Scholar] [CrossRef]
  107. Affleck, I. The quantum Hall effects, σ-models at Θ = π and quantum spin chains. Nucl. Phys. B 1985, 257, 397–406. [Google Scholar] [CrossRef]
  108. Orbach, R. Linear Antiferromagnetic Chain with Anisotropic Coupling. Phys. Rev. 1958, 112, 309–316. [Google Scholar] [CrossRef]
  109. Walker, L.R. Antiferromagnetic Linear Chain. Phys. Rev. 1959, 116, 1089–1090. [Google Scholar] [CrossRef]
  110. Baxter, R.J. One-Dimensional Anisotropic Heisenberg Chain. Phys. Rev. Lett. 1971, 26, 834. [Google Scholar] [CrossRef]
  111. Baxter, R.J. One-dimensional anisotropic Heisenberg chain. Ann. Phys. 1972, 70, 323–337. [Google Scholar] [CrossRef]
  112. Radcliffe, J.M. Some properties of coherent spin states. J. Phys. A Gen. Phys. 1971, 4, 313–323. [Google Scholar] [CrossRef]
  113. Fradkin, E. Field Theories of Condensed Matter Physics, 2nd ed.; Cambridge University Press: Cambridge, UK, 2013. [Google Scholar] [CrossRef]
  114. White, S.R.; Huse, D.A. Numerical renormalization-group study of low-lying eigenstates of the antiferromagnetic S=1 Heisenberg chain. Phys. Rev. B 1993, 48, 3844–3852. [Google Scholar] [CrossRef]
  115. Gomes, M.; Mendes, R.S.; Ribeiro, R.F.; da Silva, A.J. Gauge structure, anomalies and mass generation in a three-dimensional Thirring model. Phys. Rev. D 1991, 43, 3516–3523. [Google Scholar] [CrossRef] [PubMed]
  116. Nocedal, J.; Wright, S.J. Numerical Optimization, 2nd ed.; Springer: New York, NY, USA, 2006. [Google Scholar]
  117. Misumi, T.; Tanizaki, Y. Lattice gauge theory for the Haldane conjecture and central-branch Wilson fermion. Prog. Theor. Exp. Phys. 2020, 2020, 033B03. [Google Scholar] [CrossRef] [Green Version]
  118. Misumi, T.; Yumoto, J. Varieties and properties of central-branch Wilson fermions. Phys. Rev. D 2020, 102, 034516. [Google Scholar] [CrossRef]
  119. Verstraete, F.; Murg, V.; Cirac, J.I. Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems. Adv. Phys. 2008, 57, 143–224. [Google Scholar] [CrossRef] [Green Version]
  120. Orús, R. A practical introduction to tensor networks: Matrix product states and projected entangled pair states. Ann. Phys. 2014, 349, 117–158. [Google Scholar] [CrossRef] [Green Version]
  121. Ran, S.J.; Tirrito, E.; Peng, C.; Chen, X.; Tagliacozzo, L.; Su, G.; Lewenstein, M. Tensor Network Contractions: Methods and Applications to Quantum Many-Body Systems; Springer Nature: Berlin/Heidelberg, Germany, 2020. [Google Scholar]
  122. Cirac, J.I.; Verstraete, F. Renormalization and tensor product states in spin chains and lattices. J. Phys. A Math. Theor. 2009, 42, 504004. [Google Scholar] [CrossRef] [Green Version]
  123. Evenbly, G.; Vidal, G. Algorithms for entanglement renormalization: Boundaries, impurities and interfaces. J. Stat. Phys. 2014, 157, 931–978. [Google Scholar] [CrossRef] [Green Version]
  124. Molnar, A.; Schuch, N.; Verstraete, F.; Cirac, J.I. Approximating Gibbs states of local Hamiltonians efficiently with projected entangled pair states. Phys. Rev. B 2015, 91, 045138. [Google Scholar] [CrossRef] [Green Version]
  125. Eisert, J.; Cramer, M.; Plenio, M.B. Area laws for the entanglement entropy—A review. arXiv 2008, arXiv:0808.3773. [Google Scholar]
  126. Eisert, J. Entanglement and tensor network states. arXiv 2013, arXiv:1308.3318. [Google Scholar]
  127. White, S.R. Density matrix formulation for quantum renormalization groups. Phys. Rev. Lett. 1992, 69, 2863. [Google Scholar] [CrossRef] [PubMed]
  128. Rommer, S.; Östlund, S. Class of ansatz wave functions for one-dimensional spin systems and their relation to the density matrix renormalization group. Phys. Rev. B 1997, 55, 2164. [Google Scholar] [CrossRef] [Green Version]
  129. Dukelsky, J.; Martín-Delgado, M.A.; Nishino, T.; Sierra, G. Equivalence of the variational matrix product method and the density matrix renormalization group applied to spin chains. Europhys. Lett. (EPL) 1998, 43, 457–462. [Google Scholar] [CrossRef] [Green Version]
  130. Murg, V.; Verstraete, F.; Cirac, J.I. Variational study of hard-core bosons in a two-dimensional optical lattice using projected entangled pair states. Phys. Rev. A 2007, 75, 033605. [Google Scholar] [CrossRef] [Green Version]
  131. Verstraete, F.; Cirac, J.I. Renormalization algorithms for quantum-many body systems in two and higher dimensions. arXiv 2004, arXiv:0407066. [Google Scholar]
  132. Cirac, J.I.; Perez-Garcia, D.; Schuch, N.; Verstraete, F. Matrix product states and projected entangled pair states: Concepts, symmetries, theorems. Rev. Mod. Phys. 2021, 93, 045003. [Google Scholar] [CrossRef]
  133. Vidal, G. Classical simulation of infinite-size quantum lattice systems in one spatial dimension. Phys. Rev. Lett. 2007, 98, 070201. [Google Scholar] [CrossRef] [Green Version]
  134. Orus, R.; Vidal, G. Infinite time-evolving block decimation algorithm beyond unitary evolution. Phys. Rev. B 2008, 78, 155117. [Google Scholar] [CrossRef] [Green Version]
  135. Orús, R.; Vidal, G. Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction. Phys. Rev. B 2009, 80, 094403. [Google Scholar] [CrossRef] [Green Version]
  136. Corboz, P.; White, S.R.; Vidal, G.; Troyer, M. Stripes in the two-dimensional t-J model with infinite projected entangled-pair states. Phys. Rev. B 2011, 84, 041108. [Google Scholar] [CrossRef] [Green Version]
  137. Thies, M. Analytical solution of the Gross-Neveu model at finite density. Phys. Rev. D 2004, 69, 067703. [Google Scholar] [CrossRef] [Green Version]
  138. Lenz, J.; Pannullo, L.; Wagner, M.; Wellegehausen, B.; Wipf, A. Inhomogeneous phases in the Gross-Neveu model in 1 + 1 dimensions at finite number of flavors. Phys. Rev. D 2020, 101, 094512. [Google Scholar] [CrossRef]
Figure 1. Four-Fermi-Wilson field theories: (a) For d = 1 , the Dirac fermions indexes corresponding to the lattice and spinor degrees of freedom can be depicted as a synthetic two-leg ladder. The Wilsonian regularisation of the four-Fermi field theory can be depicted by horizontal (intraleg) and cross-link (interleg) tunnellings as well as interleg density–density interactions and an onsite energy imbalance (not shown). (b) For d = 2 , the fermions occupy a synthetic bilayer with intra- and interlayer tunnellings that depend on the direction { e j } j = 1 , 2 as well as an interlayer density–density interaction and an onsite energy imbalance (not shown). In both (a,b), if the legs or layers are understood as different spin states of the fermions, the tunnellings can be understood in terms of spin-orbit coupling in Hubbard lattice models.
Figure 1. Four-Fermi-Wilson field theories: (a) For d = 1 , the Dirac fermions indexes corresponding to the lattice and spinor degrees of freedom can be depicted as a synthetic two-leg ladder. The Wilsonian regularisation of the four-Fermi field theory can be depicted by horizontal (intraleg) and cross-link (interleg) tunnellings as well as interleg density–density interactions and an onsite energy imbalance (not shown). (b) For d = 2 , the fermions occupy a synthetic bilayer with intra- and interlayer tunnellings that depend on the direction { e j } j = 1 , 2 as well as an interlayer density–density interaction and an onsite energy imbalance (not shown). In both (a,b), if the legs or layers are understood as different spin states of the fermions, the tunnellings can be understood in terms of spin-orbit coupling in Hubbard lattice models.
Symmetry 14 00799 g001
Figure 2. 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 ) S a ( x + a 1 e 1 ) of strength J 1 , a for a { x , y , z } , here represented by solid lines on a red scale that determine 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 ) S a ( x + a j 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 shown in (a), the spins now interact vertically with strengths J 2 , a , here represented by solid lines on a blue scale that determine 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.
Figure 2. 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 ) S a ( x + a 1 e 1 ) of strength J 1 , a for a { x , y , z } , here represented by solid lines on a red scale that determine 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 ) S a ( x + a j 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 shown in (a), the spins now interact vertically with strengths J 2 , a , here represented by solid lines on a blue scale that determine 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.
Symmetry 14 00799 g002
Figure 3. Spin coherent states for bipartite lattices: (a) The one-dimensional chain Λ 1 of lattice spacing a 1 contains a 2-site unit cell with s = A (odd) and s = B (even) sites represented by blue and red circles. The low-energy properties of the half-filled chain in the strong-coupling limit can be spanned by the | s , | s states, which correspond to the north and south poles of the respective S 2 spheres. A generic spin coherent state can be described by the unit vector ω s with angles θ s , ϕ s . (b) For a rectangular lattice Λ 2 with spacings a 1 , a 2 , the two sub-lattices A , B are represented by blue and red symbols. The spin-coherent state basis is now composed of a 4-site unit cell where, in addition to the sublattice label s = A , B , we consider the left- and right- corners c = l , r , leading to θ s , c , ϕ s , c .
Figure 3. Spin coherent states for bipartite lattices: (a) The one-dimensional chain Λ 1 of lattice spacing a 1 contains a 2-site unit cell with s = A (odd) and s = B (even) sites represented by blue and red circles. The low-energy properties of the half-filled chain in the strong-coupling limit can be spanned by the | s , | s states, which correspond to the north and south poles of the respective S 2 spheres. A generic spin coherent state can be described by the unit vector ω s with angles θ s , ϕ s . (b) For a rectangular lattice Λ 2 with spacings a 1 , a 2 , the two sub-lattices A , B are represented by blue and red symbols. The spin-coherent state basis is now composed of a 4-site unit cell where, in addition to the sublattice label s = A , B , we consider the left- and right- corners c = l , r , leading to θ s , c , ϕ s , c .
Symmetry 14 00799 g003
Figure 4. Large- S Ising magnetism in the Heisenberg–Ising chain for d = 1 : (a) Ferromagnetic order parameter S y ( x ) as a function of the relative transverse field, considering various values of the Wilson parameter. For each value of r, the region with nonzero magnetisation corresponds to the long-range ordered FM y . (b) Chiral magnetic susceptibility χ y = S y ( x ) / h z , which peaks at the critical points of SSB.
Figure 4. Large- S Ising magnetism in the Heisenberg–Ising chain for d = 1 : (a) Ferromagnetic order parameter S y ( x ) as a function of the relative transverse field, considering various values of the Wilson parameter. For each value of r, the region with nonzero magnetisation corresponds to the long-range ordered FM y . (b) Chiral magnetic susceptibility χ y = S y ( x ) / h z , which peaks at the critical points of SSB.
Symmetry 14 00799 g004
Figure 5. Large- S phase diagram for d = 1 : Magnetisation contour plot, including the red stars that stand for the critical points obtained numerically in Figure 4b as well as the white dashed line for the predictions in Equation (44), which correspond to a straight line with negative unit slope h z / J 1 , y c = 1 | J 1 , x / J 1 , y | .
Figure 5. Large- S phase diagram for d = 1 : Magnetisation contour plot, including the red stars that stand for the critical points obtained numerically in Figure 4b as well as the white dashed line for the predictions in Equation (44), which correspond to a straight line with negative unit slope h z / J 1 , y c = 1 | J 1 , x / J 1 , y | .
Symmetry 14 00799 g005
Figure 6. Large- S compass magnetism for d = 2 : The large-S method predicts that the type of ferromagnetic SSB order (inversion-breaking) condensate will correspond to FM x or FM y for ξ 2 > 1 or ξ 2 < 1 , respectively. To visualise the phase diagram in a single figure, we present stacked contour plots of the difference in ferromagnetic order parameters S y ( x ) S x ( x ) as a function of the relative transverse field and the Wilson parameter via | J 1 , x / J 1 , y | = | J 2 , y / J 2 , x | = ( 1 r 2 ) / ( 1 + r 2 ) . The contour plot shows FM x in red scale, FM y on a yellow scale and PM on a blue scale. We also include dashed lines for the large-S analytical estimates (49).
Figure 6. Large- S compass magnetism for d = 2 : The large-S method predicts that the type of ferromagnetic SSB order (inversion-breaking) condensate will correspond to FM x or FM y for ξ 2 > 1 or ξ 2 < 1 , respectively. To visualise the phase diagram in a single figure, we present stacked contour plots of the difference in ferromagnetic order parameters S y ( x ) S x ( x ) as a function of the relative transverse field and the Wilson parameter via | J 1 , x / J 1 , y | = | J 2 , y / J 2 , x | = ( 1 r 2 ) / ( 1 + r 2 ) . The contour plot shows FM x in red scale, FM y on a yellow scale and PM on a blue scale. We also include dashed lines for the large-S analytical estimates (49).
Symmetry 14 00799 g006
Figure 7. Tensor network representations: (a) Matrix product states as a network of maximally entangled states | I shared between physical sites of the one-dimensional lattice to which local operations A j are applied on combined virtual space at each site. (b) Diagrammatic representation of MPS characterised by a three leg tensor A defined at every site throughout the tensor network. (c) Diagrammatic representation of PEPS corresponding to a square lattice.
Figure 7. Tensor network representations: (a) Matrix product states as a network of maximally entangled states | I shared between physical sites of the one-dimensional lattice to which local operations A j are applied on combined virtual space at each site. (b) Diagrammatic representation of MPS characterised by a three leg tensor A defined at every site throughout the tensor network. (c) Diagrammatic representation of PEPS corresponding to a square lattice.
Symmetry 14 00799 g007
Figure 8. Phase diagram of the Heisenberg–Ising chain: The phase diagram displays two regions hosting a long-range-ordered ferromagnetic phase (FM y ) and a paramagnetic phase (PM). The horizontal axis represents the magnetic field h z , whereas the vertical axis corresponds to the ratio of the tunnelling strengths J 1 , z / J 1 , y . The red stars (yellow dashed lines) show the critical points found from DMRG (self-consistent mean-field) numerics. These points are plotted on top of the contour plot of the magnetisation S y ( x ) obtained using DMRG.
Figure 8. Phase diagram of the Heisenberg–Ising chain: The phase diagram displays two regions hosting a long-range-ordered ferromagnetic phase (FM y ) and a paramagnetic phase (PM). The horizontal axis represents the magnetic field h z , whereas the vertical axis corresponds to the ratio of the tunnelling strengths J 1 , z / J 1 , y . The red stars (yellow dashed lines) show the critical points found from DMRG (self-consistent mean-field) numerics. These points are plotted on top of the contour plot of the magnetisation S y ( x ) obtained using DMRG.
Symmetry 14 00799 g008
Figure 9. Ferromagnetic and paramagnetic susceptibilities: (a) The ferromagnetic susceptibility χ M y for a fixed coupling strength J y = 1 and for different couplings J x . As the magnetic field h z is varied, it develops peaks at the critical points. In the inset, we show ferromagnetic magnetisation along the y direction. The system develops a nonzero expectation value for transverse fields below a critical value h z < h z c . (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.
Figure 9. Ferromagnetic and paramagnetic susceptibilities: (a) The ferromagnetic susceptibility χ M y for a fixed coupling strength J y = 1 and for different couplings J x . As the magnetic field h z is varied, it develops peaks at the critical points. In the inset, we show ferromagnetic magnetisation along the y direction. The system develops a nonzero expectation value for transverse fields below a critical value h z < h z c . (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.
Symmetry 14 00799 g009
Figure 10. Phase diagram of the Compass Heisenberg–Ising model: The phase diagram displays two regions hosting a long-range-ordered ferromagnetic phase (FM x ) and a paramagnetic phase (PM). The horizontal axis represents the magnetic field h z , whereas the vertical axis corresponds to the ratio of the tunnelling strengths J 2 , y / J 2 , x . The red stars (yellow dashed lines) show the critical points found from the iPEPS algorithm (large-S predictions) numerics. These points are plotted on top of the contour plot of the magnetisation S x ( x ) , using the iTEBD algorithm for the iPEPS ansatz.
Figure 10. Phase diagram of the Compass Heisenberg–Ising model: The phase diagram displays two regions hosting a long-range-ordered ferromagnetic phase (FM x ) and a paramagnetic phase (PM). The horizontal axis represents the magnetic field h z , whereas the vertical axis corresponds to the ratio of the tunnelling strengths J 2 , y / J 2 , x . The red stars (yellow dashed lines) show the critical points found from the iPEPS algorithm (large-S predictions) numerics. These points are plotted on top of the contour plot of the magnetisation S x ( x ) , using the iTEBD algorithm for the iPEPS ansatz.
Symmetry 14 00799 g010
Figure 11. Ferromagnetic susceptibility: The ferromagnetic susceptibility χ M x for J 2 , y = 1.0 and J 1 , x = 0.2 shows peaks for different Wilson parameters r and allows us to locate the critical points. In the upper inset, we show the magnetisation M x versus h z .
Figure 11. Ferromagnetic susceptibility: The ferromagnetic susceptibility χ M x for J 2 , y = 1.0 and J 1 , x = 0.2 shows peaks for different Wilson parameters r and allows us to locate the critical points. In the upper inset, we show the magnetisation M x versus h z .
Symmetry 14 00799 g011
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Tirrito, E.; Hands, S.; Bermudez, A. Large-S and Tensor-Network Methods for Strongly-Interacting Topological Insulators. Symmetry 2022, 14, 799. https://doi.org/10.3390/sym14040799

AMA Style

Tirrito E, Hands S, Bermudez A. Large-S and Tensor-Network Methods for Strongly-Interacting Topological Insulators. Symmetry. 2022; 14(4):799. https://doi.org/10.3390/sym14040799

Chicago/Turabian Style

Tirrito, Emanuele, Simon Hands, and Alejandro Bermudez. 2022. "Large-S and Tensor-Network Methods for Strongly-Interacting Topological Insulators" Symmetry 14, no. 4: 799. https://doi.org/10.3390/sym14040799

APA Style

Tirrito, E., Hands, S., & Bermudez, A. (2022). Large-S and Tensor-Network Methods for Strongly-Interacting Topological Insulators. Symmetry, 14(4), 799. https://doi.org/10.3390/sym14040799

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop