Gross-Neveu-Wilson model and correlated symmetry-protected topological phases

We show that a Wilson-type discretization of the Gross-Neveu model, a fermionic N-flavor quantum field theory displaying asymptotic freedom and chiral symmetry breaking, can serve as a playground to explore correlated symmetry-protected phases of matter using techniques borrowed from high-energy physics. A large- N study, both in the Hamiltonian and Euclidean formalisms, yields a phase diagram with trivial, topological, and symmetry-broken phases separated by critical lines that meet at a tri-critical point. We benchmark these predictions using tools from condensed matter and quantum information science, which show that the large-N method captures the essence of the phase diagram even at N = 1. Moreover, we describe a cold-atom scheme for the quantum simulation of this lattice model, which would allow to explore the single-flavor phase diagram.


I. INTRODUCTION
The understanding and classification of all possible phases of matter is one of the most important challenges of contemporary condensed-matter physics [1], and high-energy physics [2], finding also important implications in quantum information science [3]. Such a complex quest can benefit enormously from the complementary perspectives and tools developed by these different communities, calling for a crossdisciplinary dialogue that can lead to a very interesting collaborative approach. The theory of spontaneous symmetry breaking [4] and critical phenomena [5] are representative examples, where such an open dialogue has provided fundamental insight to unveil generic and universal properties in the classification of various phases of matter, and the transitions between them. However, these examples do not exhaust all possible phenomena [6], encouraging further efforts to provide a general classification encompassing other exotic orders.
Some of these studies were initially stimulated by the community working on quantum chaos [7], which looked for a complete classification of various random matrix ensembles depending on the symmetries, leading to the so-called tenfold way [8]. The ten-fold way turned out to be a fundamental tool for the classification of non-interacting phases of matter [9] which, in contrast to symmetry-broken phases, can exist within the same symmetry class [10][11][12]. In this case, transitions between different phases of matter can only occur via gap-closing continuous phase transitions, but there is neither symmetry-breaking, nor any underlying local order parameter. In contrast, these new phases are characterized by a topological invariant, the value of which changes abruptly across the symmetry-preserving critical point. This leads to the notion of symmetry-protected topological (SPT) phases [13], which includes the fermionic topological insulators and superconductors, but also other SPT phases of bosons and spins.
From a quantum-information perspective, the recent progress in the so-called tensor networks [14] has triggered the interest within this community in the general question of classifying topological phases of matter for generic interacting systems [15], including static and dynamical situations [16]. Note that, despite the considerable progress, a complete classification has so far been accomplished only for (1+1)-dimensional systems [17]. At such reduced dimensionalities, there is essentially a single gapped phase, which is trivial (i.e. it can be transformed into an uncorrelated product state by local unitaries) unless additional discrete symmetries are taken into account. Such symmetries may protect arXiv:1807.03202v2 [cond-mat.quant-gas] 10 Sep 2018 the phases, such that the states belonging to different symmetry sectors cannot be transformed into one another using local symmetry-preserving operations. A detailed understanding and characterization of the properties of these SPT phases, in the presence of interactions and strong correlations, is an open question of current interest. As argued in this work, these phases are not only relevant in condensed-matter systems, but also arise in the context of high-energy physics for certain lattice formulations of quantum field theories (QFTs).
In this paper, we focus on strongly-correlated SPT phases of a paradigmatic model of high energy physics: the Gross-Neveu model [18]. This QFT describes Dirac fermions with N flavors interacting via quartic interactions in 1 spatial and 1 time dimension, and was originally introduced as a toy model that shares several fundamental features with quantum chromodynamics. We consider a Wilson-type discretization of the QFT [19], and term the lattice version as the Gross-Neveu-Wilson model. Despite extensive studies of the GN model, a detailed characterization of strongly-correlated SPT phases has not been discussed in detail, to the best of our knowledge, neither in the large, nor in the finite N limit. The present work has the ambition of filling this gap using methods of contemporary theoretical physics and numerical simulations. Moreover, we present a scheme for the experimental realization of this discretized QFT using cold-atom quantum simulators. In this way, we hope that the Gross-Neveu model will get upgraded from a toy model used to understand some essential features of more realistic high-energy QFTs, into a cornerstone in the classification of correlated topological phases of interest in condensed matter and quantum information, which can also be explored in a realistic experiment of atomic, molecular and optical physics.
We now summarize our main results, and how they are organized in this paper: In Sec. II, we discuss generalities of the Gross-Neveu-Wilson model viewed from the complementary perspectives of high-energy, condensed-matter, and coldatom physics. This section is intended to bridge the specific knowledge gaps between these different communities, in our effort to provide a self-contained cross-disciplinary study. In Sec. III, we study the occurrence of correlated SPT phases in the model using tools common to high-energy physics. We discuss the phase diagram from the large-N expansion, including both a continuous time approach (i.e. Hamiltonian field theory on the lattice), and a discretized time approach (leading to Euclidean field theory on the lattice). This detailed study has allowed us to identify important details of the Euclidean approach, which must be carefully considered in order to understand the phase diagram of the model. In particular, we provide a neat picture where trivial gapped phases and correlated SPT phases, well-known in the condensedmatter and quantum information communities, and paritybreaking Aoki phases, well-known to the lattice field theory community, coexist in a rich phase diagram. Moreover, the large-N approach is exploited to indicate the existence of tricritical points where these three different phases are joined. We benchmark the large-N predictions using tools common to the condensed-matter and quantum-information communities, i.e. tensor-network techniques based on matrix-product-state variational ansatzs. As discussed in the text, these quasiexact numerics for a N = 1 realization of the Gross-Neveu-Wilson model confirm the large-N prediction of the phase diagram, and provide additional information that complements the large-N high-energy-inspired understanding of the model. Finally, we present a proposal for a potential experimental realization of the Gross-Neveu-Wilson model with ultra-cold atoms. In this way, relativistic models of high-energy physics could be explored with table-top non-relativistic experiments at ultra-low temperatures by focusing on low-energies and long wavelengths. The Gross-Neveu model is a relativistic QFT describing N species (flavors) of a massless Dirac field, which live in a (1+1)-dimensional spacetime and interact via four-fermion terms [18]. This model originates from its higher-dimensional counterparts, the so-called Nambu-Jona-Lasinio models [20,21], which were introduced as alternatives to non-Abelian gauge theories [22]. Pre-dating quantum chromodynamics (QCD) [23], these models offer a simplified framework to study essential features of the strong interaction, such as dynamical mass generation by chiral symmetry breaking. In addition to these features, the lower-dimensional Gross-Neveu model was introduced post-QCD as a tractable QFT displaying asymptotic freedom in a renormalizable framework. In contrast to some of its higher-dimensional cousins, this feature permits to derive rigorous results concerning the renormalization group and the convergence of perturbation theory [24].
In the continuum, this model is described by the following normal-ordered Hamiltonian H = dx : H : with Here, ψ n (x), ψ n (x) = ψ † n (x)γ 0 are two-component spinor field operators for the n-th fermionic species, and γ 0 = σ z , γ 1 = iσ y are the gamma matrices, which can be expressed in terms of Pauli matrices for a (1+1)-dimensional Minkowski spacetime, leading to the chiral matrix γ 5 = γ 0 γ 1 = σ x [25]. Therefore, the Gross-Neveu model describes a collection of N copies of massless Dirac fields coupled via the quartic interactions.
The first term in Eq. (1) corresponds to the kinetic energy of the massless Dirac fermions, where we use natural units h = c = 1, whereas the second term describes two-body interactions between pairs of fermions that scatter off each other with a strength g 2 /N. This model has a global, discrete chiral symmetry ψ n (x) → γ 5 ψ n (x), ∀x, as follows directly from the anti-commutation relations of the Dirac matrices. Additionally, a global U(N) internal symmetry becomes apparent by introducing Ψ(x) = (ψ 1 (x), · · · , ψ N (x)) t , after rewriting the Gross-Neveu Hamiltonian density as which is invariant under the transformation Ψ(x) → u ⊗ I 2 Ψ(x), ∀x, with the unitary matrix u ∈ U(N). We note that the fields have classical mass dimension d ψ = 1/2, while the interaction couplings are dimensionless d g = 0.
In the limit where the number of flavors N is very large, D. J. Gross and A. Neveu showed that this model yields a renormalizable QFT displaying asymptotic freedom, i.e. the interaction strength g 2 is a relevant perturbation in the infrared (IR), but becomes weaker at high energies in the ultraviolet (UV) limit [18]. Moreover, even if the discrete chiral symmetry prevents the fermions from acquiring a mass to all orders in perturbation theory, they showed that a mass can be dynamically generated through the spontaneous breaking of this chiral symmetry, which can be captured by large-N methods. In contrast to the Higgs mechanism, where masses can be generated by introducing additional scalar fields that undergo spontaneous symmetry breaking themselves, here a physical mass (i.e. gap) is generated dynamically as a non-perturbative consequence of the four-fermion interactions. These results are exact in the N → ∞ limit, and it is possible to calculate the leading corrections for a finite, but still large, N.
A different strategy to explore such non-perturbative effects is the so-called lattice field theory (LFT), which discretizes the fermion fields on a uniform lattice Λ s = aZ N s = {x : x/a = n ∈ Z N s }, where N s is the number of lattice sites, and a is the lattice spacing [28]. A naive discretization of the derivative of the Dirac operator yields the Hamiltonian H N = a ∑ x∈Λ s :H N :, which describes a system of interacting fermions hopping between neighboring sites of a one-dimensional lattice Here, the lattice fields fulfill the desired anti-commutation algebra in the continuum limit Unfortunately, this naive discretization also leads to spurious fermion doublers which, for g 2 = 0, correspond to massless Dirac fields appearing as long-wavelength excitations around the corners of the Brillouin zone [29]. In the present case, the Brillouin zone is BZ s = {k = 2πn/N s } = (−π/a, π/a] such that, in addition to the target massless Dirac field around k = 0, a single doubler arises around the corner k = π/a [30]. Note that, as soon as the interactions are switched on g 2 > 0, there will be scattering processes where the doubler affects the properties of the massless Dirac field, such that the continuum limit may differ from the desired QFT (1). Among several possible strategies to cope with the presence of such fermion doublers [28], K. Wilson considered introducing a momentum-dependent mass term, the socalled Wilson mass, that sends all the doublers to the cutoff of the lattice field theory [19]. In this way, one expects that these heavy fermions will not influence the universal longwavelength properties of the continuum limit.
For the Hamiltonian QFT of interest (1), this can be accomplished by introducing an additional Wilson term in the naive discretization (3) leading to H W = a∑ x∈Λ s : H W :, where which will be referred to as the Gross-Neveu-Wilson (GNW) model in this work. Here, r ∈ [0, 1] is the so-called Wilson parameter. In the continuum limit, and for g 2 = 0, the mass of the doubler around k = π/a becomes m π = 2r/a, while the Dirac fermion around k = 0 remains massless m 0 = 0. We will set r = 1 henceforth, such that the doubler mass coincides with the UV energy cutoff of the QFT. On the other hand, the Dirac field around k = 0 remains massless, and one expects that the IR limit will be governed by the desired chiral-invariant QFT. This situation gets more involved when the interactions are switched on g 2 > 0, as the additional Wilson terms (4) break explicitly the discrete chiral symmetry (i.e. rΨ(x)Ψ(x) → −rΨ(x)Ψ(x) under the discrete chiral transformation since γ 5 γ 0 γ 5 = −γ 0 ). Accordingly, the vanishing mass m 0 = 0 of the Dirac fermion around k = 0 is no longer protected by the discrete chiral symmetry, and it can become finite even for perturbative interactions in contrast to the continuum model. Since one is interested in recovering the QFT (1) for massless Dirac fermions, it is thus necessary to approach the continuum limit using a different strategy. The idea is to introduce an additional mass term in the lattice Hamiltonian (4) leading toH W = a∑ x∈Λ s :H W :, where we have introduced and m is a bare mass parameter. By tuning this mass as a function of the interaction strength m(g 2 ), one must search for a critical line m = m c (g 2 ) where the renormalized mass of the Dirac fermion around k = 0 vanishesm 0 = 0, such that the correlation length fulfills ξ a (i.e. a second-order quantum phase transition). In this case, the physical quantities of interest become independent of the underlying lattice, and one expects to recover the desired continuum QFT. The key question is to analyze if such continuum scale-invariant limit corresponds to the chiral-invariant Gross-Neveu model (1), or if a QFT of a different nature emerges in the IR limit. The answer to this question may depend on the possible phases of the lattice field theory (4) and the different critical lines in between them. Therefore, addressing this question requires a detailed non-perturbative approach using for instance large-N methods on the lattice, or Monte-Carlo methods from lattice field theory. In this work, we will present a detailed large-N analysis of the lattice GNW model, applying it to the prediction of its phase diagram, and benchmarking it with numerical simulations for the N = 1 single-flavour case.

B. Symmetry-protected topological phases for interacting fermions
A wide variety of phases transitions can be understood according to Landau's theory of spontaneous symmetry break-ing [4], which exploits the notion of symmetry and local order parameters to classify various phases of matter. Nowadays, we understand that Landau's theory does not exhaust all possibilities, as one can find different phases of matter within the same symmetry class that can only be connected via phase transitions where the symmetry is not broken. These so-called symmetry-protected topological (SPT) phases cannot be described by local order parameters, but require instead the use of certain topological invariants to characterize their groundstate (e.g. topological insulators and superconductors [31,32]). These topological invariants are in turn related to observables displaying quantized values that are robust with respect to perturbations that respect these symmetries (i.e. the topological numbers can only change via a gapclosing phase transition). Accordingly, these new phases of matter can be organized within different symmetry classes, as occurs for the fermionic topological insulators [9,33]. Despite having a gapped bulk, these insulators display a quantized conductance (e.g. integer quantum Hall effect [34]) related to a topological invariant (e.g. first Chern number [10]). A bulk-edge correspondence allows to understand this topological robustness by the appearance of current-carrying edge excitations through a band-inversion process, corresponding to mid-gap states that are exponentially localized within the boundaries of the material (e.g. one-dimensional edge modes where fermions cannot back-scatter due to disorder [35]).
The connection between SPT phases and LFTs is very natural for three-dimensional time-reversal-invariant topological insulators [31]. Here, the band-inversion process yielding the topological phase leads to an odd number of massless Dirac fermions localized within the boundaries of the material. As emphasized in [36,37], this band inversion can be understood in terms of lower-dimensional versions of domainwall fermions [38], whereby an odd number of Wilson doubler masses change their sign, contributing each with a twodimensional massless Dirac fermion localized at the boundary. In fact, we note that the Wilson-like terms in Eq. (4) arise very naturally in the low-energy description of topological insulating materials in various dimensionalities [32].
Let us now discuss how these topological effects also appear in the non-interacting limits of the GNW model (3)- (5). Here, the band inversion would occur when tuning the bare mass to lie within m ∈ (m π , m 0 ), where we recall that m 0 = 0 and m π = 2/a correspond to the masses of the Wilson fermions. To understand the SPT phase in this LFT, we consider periodic boundary conditions, such that the Hamiltonian in momentum space isH W = ∑ N n=1 ∑ k∈BZ s ψ † n (k)h k (m)ψ n (k), where we have introduced the flavor-independent singleparticle Hamiltonian By a straightforward diagonalization, one finds where ψ † n,η (k), ψ n,η (k) are the creation-annihilation operators of a fermionic excitation with flavor n in the energy band ε ± (k) = ± 1 a ma + 1 − cos ka 2 + sin 2 ka.
This band structure has a non-zero gap for m > 0, yielding an insulating phase. In order to show that this insulator is topological, and an instance of a SPT phase, we note that this band structure has an associated topological invariant that can be defined through the Berry connection A n (k) = i ε n,− (k)|∂ k |ε n,− (k) , where we have introduced the singleparticle negative-energy states |ε n,− (k) = ψ † n,− (k)|0 . In our case (6), the Berry connection can be expressed as which allows to construct a topological invariant, the so-called Zak's phase [39], as the integral of the Berry connection over the Brillouin zone. From Eq. (9), the total Zak's phase ϕ Zak = ∑ n BZ s dkA n (k) can be expressed as where θ (x) is Heaviside's step function. We note that, as occurs with the Chern number and the transverse conductivity in the quantum Hall effect [10], the topological Zak's phase can be related to an observable: the electric polarization [40,41].
Since the groundstate is constructed by filling all negativeenergy states |gs = ∏ k∈ BZ s |ε − (k) , the above integral over the whole Brillouin zone (10) characterizes the topological features of the LFT groundstate. Accordingly, this LFT hosts a SPT phase in the parameter regime m π < m < m 0 for N odd, which coincides with the band-inversion regime introduced above. This regime can be interpreted as the result of a massinversion process, whereby the mass of some of the Wilson fermions gets inverted. This becomes apparent after rewriting the Zak's phase in terms of the N Wilson masses m 0 = m,m π = m + 2/a.
Indeed, a non-trivial topological invariant (i.e. ϕ Zak /2π ∈ Z) can only be achieved when an odd number of fermion doubler pairs display a different mass sign We note that this SPT phase can be identified with a one-dimensional topological insulator in the so-called chiralorthogonal BDI class [9,32,33], which would display zeroenergy modes localized at the edges of the chain for open boundary conditions. Note that this chiral symmetry class is not related to the standard notion of chirality in QCD, which is indeed broken by the GNW model. Instead, it is related to the ten-fold Cartan's classification of symmetric spaces, and its connection to single-particle Hamiltonian via the time-evolution operator [9]. For the non-interacting GNW single-particle Hamiltonian (6), we find that time-reversal T yields T † h −k (m) * T = h k (m) where T = −iσ x σ y = γ 0 , and charge-conjugation C leads to C † h −k (m) * C = −h k (m) where C = iσ z σ y = γ 5 [25]. The combination of these two antiunitary symmetries is called chiral, or sub-lattice, symmetry S = TC, and yields S † h k (m)S = −h k (m) with S = γ 0 γ 5 = γ 1 . To avoid confusion with the chiral symmetry of high-energy physics, which is a fundamental ingredient in low-energy effective descriptions of QCD, and pivotal in our previous discussion of the GNW model, we will refer to the S as the sublattice symmetry. Since T 2 = C 2 = S 2 = +1, the corresponding GNW topological insulator with an odd number N of fermion flavors (12) is in the BDI class.
We note that the ten symmetric spaces that classify the topological insulators/superconductors also correspond to the target spaces of effective non-linear-sigma model describing the long-wavelength properties of the edge/boundary. When such a non-linear-sigma model includes a topological term, the edge modes are robust and evade Anderson localization in the presence of symmetry-preserving disorder [9]. This perspective allows us to understand the difference of N even/odd in the GNW model. For N even, there can be a symmetrypreserving disorder that couples the different flavors of the edge states, leading to scattering/localization and destroying the BDI topological protection. On the other hand, for N odd, at least one of the edge modes will remain robust against interflavor scattering, and thus evade Anderson localization. We note that similar parity effects can occur also in models with more than one fermion doubler in the regime where an even number of Wilson masses gets inverted, as occurs for higherdimensional time-reversal topological insulators [36].
In contrast to the LFT perspective described in Sec. II A, where one is mainly interested in searching for the secondorder quantum phase transitions to recover a continuum limit described by the QFT of interest (1); the study of symmetry-protected topological phases focuses on the topological gapped phases away from criticality. Interestingly, even in the non-interacting regime, the emerging QFTs governing their response to external fields turns out to be very different from the original discretized QFT, and can be described in terms of topological quantum field theories (e.g. Chern-Simons or axion QFTs) [42]. A generic question of current interest in the study of SPT phases is to explore the interplay of topological features and strong-correlation effects as interactions between the fermions are switched on [43].
For the GNW lattice model (3)-(5), the interactions do not modify the symmetry class as ΨΨ → ΨΨ, and Ψγ 5 Ψ → ∓Ψγ 5 Ψ under time-reversal and charge-conjugation transformations, respectively [25]. Accordingly, the quartic terms in Eq. (2), or its chiral extension introduced in Eq. (16) below, do not modify the aforementioned BDI symmetry class. A question of potential interest for both the SPT and LFT communities is the precise determination of the critical lines of the lattice model m c (g 2 ) for non-perturbative interactions. From this knowledge, the LFT community can explore the nature of the continuum QFT in the vicinity of the critical line, while the SPT community may study how the topological phase is modified in presence of interactions. As argued above, a possible tool to study non-perturbative effects could be large-N methods, or Monte-Carlo methods in Euclidean lattice field theory. We remark, however, that the standard Euclidean ap-proach where time is also discretized [28] can lead to qualitative differences of the phase diagram in the (m, g 2 ) plane. Discretizing time introduces additional fermion doublers, which may lead to additional critical lines that are not present in the Hamiltonian approach (3)-(5) with continuous time [44]. Although this is not relevant when one is only interested in the nature of the continuum QFT, it will be of relevance for topological insulators where one is interested in the finite region of phase space with the topological gapped phase. In this work, we will show that special care in the Euclidean lattice formulation is required in order to recover the relevant phase diagram.
For one-dimensional models, the study of lattice field theories in the Hamiltonian approach can be efficiently accomplished using variational methods based on Matrix Product States (MPS) [45]. In this work, we shall confront predictions of the large-N approximation with results from MPS numerical methods for the study of topological insulating phases in the GNW model with Wilson fermions.

C. Cold-atom quantum simulators of high-energy physics
As an alternative to Monte-Carlo numerical methods in lattice field theory, one may follow R. P. Feynman's insight [46], and develop schemes to control a quantum-mechanical device such that its dynamics reproduces faithfully that of the model of interest (i.e. quantum simulation). From this perspective, a very appealing application of the future fault-tolerant quantum computers will be their ability to function as universal quantum simulators [47] that can address complicated quantum many-body problems relevant for different disciplines of physics and chemistry. Prior to the development of quantum error correction and large-scale fault-tolerant quantum computers, one may consider building special-purpose quantum simulators that are designed to tackle a particular family of models. This is the case of cold-atom quantum simulators of lattice models [48,49], where neutral atoms are laser-cooled to very low temperatures in deep optical lattices [50].
In the continuum, neutral-atom systems are typically described by a Hamiltonian QFT, albeit a non-relativistic one [50], with H = d 3 x : (H 0 + V I ) : containing are field operators that createannihilate an atom of the n-th species in the internal state σ . This Hamiltonian contains (i) the kinetic energy for a multispecies gas of Alkali atoms of mass m n , where n ∈ {1, · · · N sp } labels the atomic species/isotope; (ii) the internal energy ε n σ of the atomic groundstate manifold, which typically consists of various hyperfine levels characterized by the quantum numbers associated to the total angular momentum σ ∈ {F, M F }; and (iii) the single-particle terms V σ ,σ n (x x x), which contain the trapping potential that confines the n-th atomic species and, possibly, additional radiation-induced terms that drive transitions between the different atomic levels σ → σ . In particular, we shall be interested in periodic trapping potentials due to the ac-Stark-shift of pairs of retro-reflected laser beams, which will depend on the atomic species, but not on the particular hyperfine level (i.e. state-independent optical lattices). We also consider laser-induced Raman transitions via highly off-resonant excited states. Altogether, this leads to where V n 0,ν is the ac-Stark shift for the n-th atomic species stemming from the retro-reflected beams with wave-vector k L,ν along the ν-axis. Additionally, ω n,ν is the frequency of a residual harmonic trapping due to the intensity profile of the lasers. Finally, Ω n,l σ ,σ is the two-photon Rabi frequency for the Raman transition induced by the l-th pair of laser beams with wave-vector (frequency) difference ∆ ∆ ∆k k k l (∆ω l ), and phase ϕ l .
In addition, at sufficiently low temperatures, the neutral atoms also interact by contact scattering processes leading to where the interaction strengths U n,n σ ,σ depends on the s-wave scattering lengths a σ σ of the corresponding channels, some of which can be controlled by Feshbach resonances [50]. We also note that fully-symmetric interactions between all species can be achieved by using alkali-earth atoms [51], which could be an interesting property for the experimental realization of higher number of flavors N in the Gross-Neveu-Wilson model.
As announced above, in the regime of deep optical lattices V n 0,ν E n R = k k k 2 L /2m n , one can introduce the basis of so-called Wannier functions, which are localized to the minima of the potential, and show that this non-relativistic QFT yields a family of Hubbard-type models with tunable parameters [52,53]. Therefore, by doing controlled table-top experiments with cold atoms, it becomes possible to explore the physics of strongly-correlated electrons in solids, which has opened a fruitful avenue of research in quantum simulations of condensed-matter models [48,49]. More recently, several works have explored the possibility of extending this cold-atom Hubbard toolbox [54] to the quantum simulation of high-energy physics, including relativistic QFTs [36,[55][56][57][58][59], gauge field theories [60][61][62][63][64], theories for coupled Higgs and gauge fields [65,66], and also theories of relativistic fermions interacting with Abelian/non-Abelian gauge fields [67][68][69].
In this work, we shall be concerned with a cold-atom realization of the Gross-Neveu model using a Wilson-fermion discretization (3)-(5). We note that there are cold-atom proposals to implement this QFT (1) with optical superlattices lattices by a different discretization [56], via the so-called staggered fermions [30,44]. Since we are interested in the connection of this model with correlated SPT phases, we will instead focus on the Wilson-fermion approach of Eqs. (3)-(5). Building on previous proposals for the quantum simulation of Wilson fermions [36,[71][72][73], we present in this work a simplified scheme to realize the GNW model using a two-component single-species Fermi gas confined in a one-dimensional optical lattice with laser-assisted tunneling.

III. CORRELATED SYMMETRY-PROTECTED TOPOLOGICAL PHASES IN THE GROSS-NEVEU-WILSON MODEL
A. Phase diagram from the large-N expansion As advanced in the previous sections, our goal is to determine the critical lines of the GNW model (3)-(5) as a function of the coupling strength m(g 2 ) for non-perturbative interactions. We start by developing a large-N expansion for the partition function Z = Tr exp −βH W , where β = 1/T is the inverse temperature for k B = 1. In the continuum, large-N methods were first employed by Gross and Neveu to prove that the groundstate of their eponymous model (1) displays a nonzero vacuum expectation value σ 0 = Ψ(x)Ψ(x) = 0 ∀x, as soon as a non-vanishing interaction g 2 > 0 is switched on [18]. In this way, the discrete chiral symmetry This non-perturbative result can be obtained using functional techniques to calculate an effective action for an auxiliary bosonic σ (x) field, which condenses due to the formation of particle-anti-particle pairs, and acquires a non-zero expectation value σ (x) = σ 0 = 0 in the chirally-broken phase. On the lattice (3)-(4), similar results are recovered in the continuum limit [74,75], provided that the additional bare mass (5) is adjusted to recover the discrete chiral symmetry.
Let us now comment on a generalization of the GNW model, where the above discrete chiral symmetry is upgraded to a continuous one Ψ(x) → I N ⊗e iθ γ 5 Ψ(x), ∀θ ∈ [0, 2π) [18]. This requires a modified four-fermion term In this case, in addition to the σ (x) field, it is natural to introduce an additional bosonic field Π(x), obtaining an effective action for both fields in the large-N limit. In Ref. [76], S. Aoki showed that the large-N results with lattice Wilson fermions lead to a richer phase diagram displaying new regions where a discrete parity symmetry In this case, the particleanti-particle pairs lead to the so-called pseudoscalar condensate Π(x) = Π 0 = 0, which necessarily breaks the parity transformation of the corresponding fermion bilinear due to the vacuum expectation value Interestingly, these results on the chiral GNW model were used to conjecture that these socalled Aoki phases would also appear in the phase diagram of lattice quantum chromodynamics [76]. However, in this context, these Aoki phases are considered as unphysical lattice artifacts not present in the continuum QFT.
In this section, we discuss the role of such Aoki phases in the GN model (16) with a Wilson-type discretization (3)-(5), and their interplay with the topological insulating phases discussed in the previous sections. In the context of symmetryprotected topological phases, such Aoki phases are not artifacts, but become instead physical phases of matter that shall delimit the region of the phase diagram that hosts a correlated SPT phase. Moreover, from the perspective of a cold-atom implementation, these phases might also be observed in future table-top experiments. We also note that the appearance of Aoki phases is not restricted to the GNW model, but also occurs in strong-coupling calculations of U(1) Wilson-type lattice gauge theories [37,77], which can be used to model the strongly-interacting limit of higher-dimensional topological insulators with long-range Coulomb interactions.
We remark that, in the limit of a single fermion flavor N = 1, which is the relevant case for the cold-atom implementation, the four-fermion interactions of Eq. (3) can be rewritten as (17) which follows from a so-called Fierz identity in the language of relativistic QFTs. Accordingly, besides a change in the coupling constant g 2 → g 2 /2, there is no further distinction between the N = 1 GNW model with discrete or continuous symmetry, such that the previous Aoki phases could in principle also occur in this limiting case. However, since their prediction is based on the N → ∞ results, we will have to benchmark large-N methods with other non-perturbative approaches valid for N = 1 (e.g. MPS numerical simulations or a potential cold-atom quantum simulation). Regarding the first approach, and give a detailed comparison of the large-N predictions with the MPS results of the phase diagram.

Continuous time: Hamiltonian field theory on the lattice
Let us first discuss the large-N phase diagram of the GNW model using a functional-integral representation of the partition function with a continuum Euclidean (i.e. imaginary) time τ. Introducing fermionic coherent states by means of mutually anti-commuting Grassmann variables Ψ k (τ), Ψ k (τ), which are defined at each point of the Brilllouin zone k ∈ BZ s and for each imaginary time τ ∈ (0, β ) [78], one can readily express the finite-temperature partition function as where the Euclidean action is is defined in terms of the singleparticle Hamiltonian in Eq. (6). Moreover, V g (Ψ , Ψ) results from substituting the fermion field operators by the Grassmann variables in the normal-ordered interaction (17), which leads to quartic interactions Let us note that the propagator associated to the free part of the action displays two poles at k ∈ {0, π/a} when −ma ∈ {0, 2}, which correspond to the aforementioned Dirac fermions around the corners of the Brillouin zone.
The first step in the large-N approximation is to introduce two auxiliary real scalar fields σ (x), Π(x) with classical mass dimension d σ = d Π = 1, such that the partition function can be expressed as a new functional integral over both Grassmann and real auxiliary fields Ψ] up to an irrelevant constant, such that the thermodynamic properties of the system are not modified by the introduction of the auxiliary fields. The idea is to chose a particular action where the four-fermion terms can be understood as effective interactions carried by the auxiliary bosonic fields. Moreover, assuming that these fields are homogeneous, the new action becomes whereH k = I N ⊗h k , and the fermionic single-particle Hamiltonian now depends on the auxiliary bosonic fields Essentially, the σ field modifies the mass term of the Dirac fermion, and a vacuum expectation value of the former would thus renormalize the fermion mass, resembling the dynamical mass generation of the continuum model. The second step in the large-N approximation is to integrate over the fermionic Grassmann fields, obtaining an effective action for the auxiliary bosons Z = [dσ dΠ]e −NS eff [σ ,Π] . This step can be readily performed since the Grassmann integral is Gaussian, which leads to where L = N s a is the length of the chain, and we have introduced an abbreviation for the integral over momentum and 2π , assuming already the zero-temperature limit which is the regime of interest of this work. Here, the energies of the new fermionic single-particle Hamiltonian ε k (m+σ , Π) have been expressed in terms of the function When the number of fermion flavors is very large N → ∞, the partition function Z = [dσ dΠ]e −NS eff [σ ,Π] with the effective action (22) yields a groundstate obtained from the saddle point equations ∂ σ S eff | (σ 0 ,Π 0 ) = ∂ Π S eff | (σ 0 ,Π 0 ) = 0. Nonvanishing values of σ 0 , Π 0 are related to the breaking of the discrete chiral or parity symmetries discussed above. For instance, the boundary of the aforementioned Aoki phases can be obtained from the self-consistent solution of these saddlepoint equations imposing Π 0 = 0. Using contour techniques for the frequency integrals, and substituting ka → k in the momentum integrals, we can express the pair of saddle-point equations as follows Here, we have used the complete elliptic integrals of the first and second kind as well as the following parameters In general, the solution of the pair of gap equations (24) must be performed numerically, and leads to the critical lines that delimit the Aoki phase (i.e. solid green lines in Fig. 1). These lines can be interpreted as different flows of the bare mass m c (g 2 ) that determine the second-order phase transitions where a scale-invariant QFT should emerge. Note that this figure displays a clear reflection symmetry with respect to the axis −ma = 1. In fact, using the expression of the elliptic integrals in terms of hypergeometric functions [79], it follows that 1)), which can be exploited to show that the gap equations (24) can be rewritten as These gap equations can now be related to the original ones in Eq. (24) under the following transformation which corresponds to the aforementioned reflection symmetry about −ma = 1, and leads to η 0 → −η 0 and θ 0 → θ 0 /(θ 0 −1). Accordingly, there should only be three distinct phases in the regime −ma ∈ [0, 2], with the Aoki phase being completely absent for ma > 0 and ma < −2.
To make a connection to the continuum results [18], and interpret this phase diagram in light of the symmetry-protected topological phases of Sec. II B, we note that a solution to the gap equations (24) can be found analytically in the regime of small interactions and masses g 2 , |ma| 1. In this case, one can assume that η 0 = 1 + δ η 0 with |δ η 0 | 1, and perform a Taylor expansion of Eq. (24) to find that the σ -field acquires the following non-zero vacuum expectation valuẽ The first contribution stems for the perturbative renormalization of the bare mass ma ≈ −g 2 /2π, while the 1/g 2 behavior of the second contribution highlights that the large-N expansion captures non-perturbative effects, recalling the chiral symmetry breaking by dynamical mass generation of the continuum case [18]. We also note that, as the UV cutoff is removed a → 0, the interaction strength must decrease g 2 → 0 to maintain a finite scalar condensate (29), which shows that the continuum GNW model is an asymptotically free QFT. As announced above, such a vacuum expectation value (29) then leads to a small renormalization of the Wilson masses (11),m k →m k (g 2 ), valid in the regime g 2 , |m| 1. We can thus ascertain that the large mass of the fermion doubler will only be perturbed slightly, remaining thus at the cutoff, and maintaining sgn(m π (g 2 )) = +1. Conversely, the sign of the light-fermion mass sgn(m 0 (g 2 )) may indeed change as the interactions g 2 are increased. According to Eq. (12), we can write the topological invariant in this regime as ϕ Zak = 1 2 Nπ 1 − sgn m 0 (g 2 ) , such that the region hosting a correlated BDI topological-insulating groundstate corresponds to the parameter region withm 0 (g 2 ) < 0.
In order to locate this region, we substitute the saddle-point solution (29) into Eq. (20), and perform a long-wavelength approximation |k| < Λ 1/a, yielding the effective freefermion actioñ up to an irrelevant constant. Here, we have introduced , where the single-particle Hamiltonian for a massive Dirac fermion is which allows us to identify the renormalized Wilson mass The leftmost red dashed line of Fig. 1 corresponds to the points where this renormalized mass vanishesm 0 (g 2 ) = 0. We note that this analytical solution matches the lower critical line obtained by the numerical solution of the gap equations (24) remarkably well, even considerably beyond the perturbative regime g 2 , |ma| 1. Following Eq. (32), the area below this line fulfills −m >σ 0 , such that the interacting Dirac fermion has a negative renormalized massm 0 (g 2 ) < 0, leading to ϕ Zak = Nπ and to an SPT phase for N odd.
An analogous behavior can be found in the regime g 2 , |ma+ 2| 1, where the light fermion is around k = π/a, while the heavy one corresponds to k = 0 (i.e. the Wilson fermions interchange their roles). Using the previous symmetry (28) to locate the critical linem π (g 2 ) = 0 in this parameter regime, we can readily predict the value of this renormalized mass m π = m + 2/a →m π (g 2 ) = m + 2/a −σ 0 . The vanishing of this mass leads to the rightmost red dashed line of Fig. 1, which again agrees very well with the numerical solution of the gap equations. Since the heavy fermion around k = 0 has a large negative mass, the topological invariant becomes ϕ Zak = 1 2 Nπ sgn m π (g 2 ) + 1 , and one can identify the symmetry-protected phase displaying ϕ Zak = Nπ for N odd, with the parameter region fulfillingm π (g 2 ) > 0, and thus −m < 2 −σ 0 (i.e. shaded yellow area below the dashed line).
At larger couplings and intermediate masses, we must resort to the numerical solution of the gap equations, and search for a region of phase diagram that can be adiabatically connected to these two areas hosting a topological phase. This is precisely the shaded yellow lobe of Fig. 1, which is separated from other phases by a gap-closing line. The area above these lines, given bym 0 (g 2 ) > 0, andm π (g 2 ) < 0, determines a regime where both renormalized Wilson masses have the same sign, such that the gapped phase has no topological features, corresponding either to a trivial band insulator (grey area in Fig. 1), or to the aforementioned Aoki phase where the Z 2 parity symmetry Ψ(x) → I N ⊗ γ 0 Ψ(x) is spontaneously broken (green area in Fig. 1).

Discretized time: Euclidean field theory on the lattice
We now move on to the discussion of the large-N phase diagram of the GNW lattice model using a discretized Euclidean time x 0 = τ. This is the most common formalism in lattice field theory computations [28], and can become the starting point to apply other methods such as Monte-Carlo numerical techniques. As emphasized below, it will be important to understand the connection between the lattice and Hamiltonian approaches, requiring a careful treatment of the continuumtime limit to understand lattice artifacts that can change qualitatively the shape of the phase diagram.
In Euclidean LFT, both space-and time-like coordinates {x ν } ν=0,1 are discretized into an Euclidean lattice Λ E = {x x x : x 0 /a 0 = n τ ∈ Z N τ , x 1 /a 1 = n s ∈ Z N s }, where N τ (N s ) is the number of lattice sites in the time (space) -like direction, and a 0 (a 1 ) is the corresponding lattice spacing. Therefore, a similar discussion to the one around Eqs. (3)-(5) must also be applied to the Euclidean time derivative appearing in the action (18), such that nearest-neighbor hoppings along the timelike direction also appear. Introducing fermionic coherent states on the Euclidean lattice, and their corresponding Grassmann variables Ψ x x x , Ψ x x x , the finite-temperature partition function can be expressed as Here, the action is divided into: (i) the free quadratic term which is expressed in terms of the Euclidean gamma matriceŝ γ 0 = γ 0 ,γ 1 = iγ 1 , and the unit vectors {e ν } of a rectangular lattice; and (ii) the interacting quartic term which is expressed in terms of the chiral matrixγ 5 = γ 5 .
(i) Lattice approach with dimensionless fields: Let us note that, in the lattice Wilson approach [80], it is customary to work with dimensionless fields ψ x x x = √ a 0 + a 1 Ψ x x x , and rewrite the action as follows is expressed in term of dimensionless tunnelings κ ν , and the dimensionless massm. Similarly, the interacting term is obtained from Eq. (35) by substituting the fields Ψ → ψ and the coupling constant g →g by the dimensionless ones.
Since the Grassmann variables must fulfill periodic (antiperiodic) boundary conditions along the space (time) -like directions, one can move into momentum space ψ k k k , ψ k k k , where the dimensionless quasi-momenta belong to the Euclidean Brillouin zone BZ E = {k k k : k 0 = 2π(n τ + 1/2)/N τ , k 1 = 2πn s /N s } = (0, 2π] 2 . Then, one can rewrite the action as where we have introduced S k k k (m) = I N ⊗ s k k k (m), together with the single-flavor action Let us note that, in contrast to the continuum-time free action (18), this Euclidean action leads to a propagator with four poles at k k k ∈ {(0, 0), (0, π), (π, 0), (π, π)} when the bare mass equals −m ∈ {0, 4κ 0 , 4κ 1 , 4(κ 0 + κ 1 )}, each of which corresponds to a long-wavelength Dirac fermion. Accordingly, there is an additional doubling due to the discretization of the the Euclidean time direction (i.e. the extra fermions with k 0 = π shall be referred to as time doublers). At this point, the discussion parallels that of the Hamiltonian formalism of Sec. III A 1 via the corresponding steps for the large-N approximation. First, the auxiliary dimensionless lattice fieldsσ x x x ,Π x x x are introduced, such that the action can be rewritten as Here, we have assumed again that the auxiliary fields are homogeneous, introducingS k k k (σ ,Π) = I N ⊗s k k k (σ ,Π), such that the new single-flavor action can be obtained from Eq. (38) usings k k k (m +σ , Π) = s k k k (m +σ ) + iγ 5Π . The second and third steps are the same, since the action is quadratic in Grassmann fields, and the saddle-point solutions control the large-N limit. In this case, the gap equations can be expressed as which are equivalent to those derived in [74] upon a different definition of the microscopic couplings.
We have solved this system of non-linear equations for different Euclidean lattices with N τ = ξ N s , setting N s = 512 sites in the space-like direction, and using ξ = a 1 a 0 ∈ {1, 2, 4, · · · , 128} to approach the time-continuum limit ξ → ∞ (see Fig. 2). Let us note that the dimensionless tunnelings can be expressed in terms of the anisotropy parameter as κ 0 = ξ /2(1 + ξ ), and κ 1 = 1/2(1 + ξ ). At this point, it is worth mentioning that the number of lattice sites in the timelike direction N τ is also modified in the LFT community to explore non-zero temperatures. In that case, however, the κ ν parameters remain constant as N τ is varied (i.e. the Euclidean lattice is rectangular, but the unit vectors remain the same).
In Fig. 2(a), we represent the solution of the gap equations for the isotropic lattice ξ = 1, such that κ 0 = κ 1 = 1 4 . We note that the characteristic trident-shaped phase diagram is in qualitative agreement with the results of S. Aoki [76]. In order to interpret this phase diagram in terms of the symmetryprotected topological phases, let us recall the distribution of the poles described below Eq. (38). Atg 2 = 0, we observe that the critical points separating the different phases correspond to −m ∈ {0, 1, 2}, which lie exactly at the aforementioned poles signaling the massless Dirac fermions. For −m ∈ (0, 1), the only Dirac fermion with a negative mass is that around we see that ϕ Zak = Nπ for −m ∈ (0, 1), corresponding to the BDI topological insulator for N odd. For −m ∈ (1, 2), the Wilson fermions around k k k = (0, π) and k k k = (π, 0) also invert their masses, leading to ϕ Zak = −Nπ, and yielding again an BDI topological insulator for N odd. These two areas, extend on to the neighboring lobes of Fig. 2(a) using a similar reasoning as the one presented around Eq. (30). Therefore, the whole region below the trident that delimits the parity-broken Aoki phase corresponds to the BDI topological insulator. We note, however, that the black dashed lines in this figure, and subsequent ones, do not follow from the solution of the large-N gap equations, but are included as a useful guide to the eye to delimit the SPT phases. In Sec. III B below, we will show that they indeed correspond to a critical line delimiting the SPT phase of a carefully-defined time-continuum limit. Let us start exploring how this phase diagram changes as the time-continuum limit is approached, and compare the results to those of Fig. 1. In Fig. 2(b), we represent the phase boundaries for an increasing number of lattice sites N τ = ξ N s with anisotropies ξ ∈ {1, 2, 4, 8}. Here, one can observe how the central prong of the Aoki phase separating the topologicalinsulating lobes is split into two peaks, each of which goes in a different direction as ξ is increased. We note that this behavior differs markedly from the finite-temperature studies, which show that the lobe structure disappears completely as N τ is varied [76]. Therefore, the anisotropy in the lattice constants gives rise to a different playground, which must be understood in terms of the symmetry-protected topological phases. Since κ 0 → 1/2, while κ 1 → 0, as the anisotropy ξ → ∞, one can identify the left-moving prong with the pole at k k k = (0, π) with mass −ma → 4κ 1 → 0, and thus approaching the lower left corner. Similarly, the right-moving one can be identified with the pole at k k k = (π, 0) with mass −ma = 4κ 0 → 2 approaching the lower right corner. As a result of this movement, and considering the signs of the corresponding Wilson masses, one finds that the region between these two poles correspond to a situation where both space-(time-) like doublers have a negative (positive) mass, such that the topological invariant vanishes ϕ Zak = 0, and one gets a trivial band insulator. Unfortunately, as the anisotropy increases, the two BDI topological lobes get smaller and smaller, such that the symmetryprotected topological phases vanish as we approach the timecontinuum limit, and the central lobe corresponds to a trivial band insulator (see Fig. 2(c)).
This result seems to be in contradiction with our findings for the Hamiltonian formalism in Fig. 1, which predict that the central lobe should correspond to the correlated SPT phase with ϕ Zak = Nπ. Moreover, since each of the two prongs now contain a pair of massless Dirac fermions, the continuum QFT that should emerge in the long-wavelength limit is no-longer that of the Gross-Neveu model for N flavors, but rather that of the Gross-Neveu model for 2N flavors, which would indeed modify the universal features of the phase transition, and not only the non-universal shape of the critical line. As mentioned at the beginning of this section, the Euclidean approach can lead to lattice artifacts that can modify qualitatively the phase diagram, and a detailed and careful account of the timecontinuum limit is required to understand them. We address precisely this issue in the two following subsections.
(ii) Large-N phase diagrams with rescaled couplings: We have found that one of the problems leading to the apparent contradiction between the phase diagrams is the standard use of dimensionless quantities in the Euclidean lattice approach (36). A detailed derivation of this action, which starts from the original action (33) rescaling the fields, shows that the dimensionless parameters are related to the original ones by the following expressioñ Although apparently innocuous, this rescaling changes qualitatively the shape of the phase diagram (see Fig. 3). In order to understand the main features of this phase diagram, the location of the non-interacting poles will be very useful again. For instance, at g 2 = 0, we note that the pole at −m = 4κ 0 gets mapped into −ma 1 = 4(1 + ξ )κ 1 . Therefore, as the time-continuum limit is approached, this pole tends to −ma 1 → 2 as ξ → ∞, and no longer to the origin. Likewise, both time-like doublers at −m ∈ {4κ 0 , 4(κ 0 +κ 1 )} are mapped into −ma 1 → ∞ in the time-continuum limit. Accordingly, in the region of interest displayed in Fig. 3, these time doublers have a very large positive mass. Inspecting the sign of the corresponding Wilson masses, we can conclude that the region −ma 1 ∈ (0, 2) will host an BDI topological insulator, while a trivial insulator will set in for −ma 1 > 2.
Following a similar reasoning as in previous subsections, we know that these critical points surrounding the topological phase will flow as the interactions are switched on and the σ field acquires an non-zero vacuum expectation value. Accordingly, we identify the lobe of Fig. 3 as the BDI topological insulator that also appeared in the continuum-time Hamiltonian formalism of Fig. 1. Moreover, the universal features are now in agreement as the critical lines are controlled by a single pole, and the long-wavelength limit should now be controlled Let us remark that, although the rescaled solution looks somewhat closer to the Hamiltonian results, there are still qualitative differences in the lattice approach that deserve a deeper understanding. For instance, the phase diagram does no longer display the mirror symmetry about −ma 1 = 1 (28).
(iii) Continuum limit and connection to the Hamiltonian approach: In order to understand these differences, and the connection to the gap equations continuum limit (24), let us consider the original action with dimensional fields (33). Following the same steps as before, one can integrate the fermion fields, Z = [dσ dΠ]e −NS eff [σ ,Π] , finding the following effective action which is the Euclidean lattice version of Eq. (22). Here, we have introduced the corresponding lengths L τ = N τ a 0 , L s = N s a 1 , together with the following function If we now take the limit of N → ∞, the saddle point conditions ∂ σ S eff | (σ 0 ,Π 0 ) = ∂ Π S eff | (σ 0 ,Π 0 ) = 0 lead to the following pair of gap equations, which are equivalent to Eqs. (40)- (41) but using dimensional couplings and dimensional fields, In order to make a connection to the gap equations obtained with the Hamiltonian formalism (24), we should take the continuum limit in the imaginary time direction N τ → ∞, and a 0 → 0, such that L τ = L s remains constant imposing ξ = a 1 /a 0 → ∞. To deal with the additional time doublers mentioned above, let us introduce a UV cutoff Λ τ 1/a 0 , and make a long-wavelength approximation around k 0 ∈ {0, π/a 0 }. We find that the gap equation (40) becomes where we have used the single-particle energies of Eq. (23) and the spatial Brillouin zone, after identifying a = a 1 . We note that the first line of this expression comes from the contribution around k 0 = 0, while the second line stems from the time doublers around k 0 = π/a 0 . We observe that the effective Wilson mass of these doublers becomes very large in the continuum limit m + 2/a 0 → ∞ if one keeps the bare mass m non-zero. Hence, these doublers become very massive, and their contribution to above gap equation should become vanishingly small as described below Eq. (43). To prove that, let us get rid of the cutoff Λ τ → ∞, and use ∑ |k 0 |<Λ τ → L τ ∞ −∞ dk 0 /2π. After performing the integral using contour techniques, we directly obtain 1 g 2 = π −π dk 1 4π 1 (ma + σ a + (1 − cos k 1 a)) 2 + sin 2 k 1 a + π −π dk 1 4π 1 (ma + σ a + (1 + 2ξ − cos k 1 a)) 2 + sin 2 k 1 a , where we have also taken the continuum limit in the spacelike direction. Using the definition of the complete elliptic integrals (25), this equation can be expressed Here, we have used the parameters of Eq. (26), together with which determine the contribution of the time doublers to the gap equation (i.e. second term of Eq. (48)). In the continuum limit, we take ξ → ∞, such thatη 0 → ∞, andθ 0 → 0. This makes K(θ 0 ) → π/2, such that the time-doubler contribution vanishes, and we recover exactly the gap equation of the Hamiltonian approach (24).
The continuum limit of the remaining gap equation (41) follows the same lines: we perform a long-wavelength approximation around the time doublers, let the cutoff Λ τ → ∞, and use contour integration to find Using the definition of the complete elliptic integrals (25), this equation becomes where the contribution of the time doublers is expressed in the second line. In this case, taking the time-continuum limit ξ → ∞, such that E(θ 0 ) → π/2, leads to which contains an additional −g 2 /2 term with respect to the gap equation of the Hamiltonian formalism (24). We thus find that, in contrast to the first gap equation (49), the contribution of the time doublers is no longer vanishing in this case, but can instead be understood as a finite renormalization of the bare mass ma → m r a = ma + g 2 /2.
It is precisely this renormalization which is responsible for the lack of the mirror symmetry (28) in Fig. 3, and its qualitative difference with respect to the Hamiltonian prediction of Fig. 1. These results can thus help us to identify the corresponding mirror symmetry, which is no longer about the vertical line −ma = 1, but instead about −m r a = 1, which corresponds to the red dashed line −ma = 1 + g 2 /2 of Fig. 3.
To study in more detail the onset of this symmetry in the continuum limit, and the quantitative agreement with the Hamiltonian prediction, we plot the phase diagram with the corresponding renormalized mass in Fig. 4, and superimpose the continuum-time prediction of Fig. 1. This figure shows the clear agreement between both approaches, and highlights the importance of performing a careful analysis of the continuum limit in order to avoid Euclidean lattice artifacts that can lead to qualitatively different predictions, even questioning the universal aspects of the emerging QFTs (see Fig. 2). It also highlights the fact that the time doublers, despite becoming infinitely heavy in the continuum limit, can leave an imprint in the non-universal properties of the low-energy phase diagram, such as the particular value of the critical points (see the tilted phase diagram of Fig. 4). From the perspective of the renormalization group, this effect does not come as a surprise, since the time doublers lie at the cutoff of the continuum-time limit of the lattice field theory, and their integration can thus renormalize the parameters of the long-wavelength light-fermion modes. In this case, a careful analysis of the gap equations has allowed us to extract an additive renormalization δ m = g 2 /2a which, as usual in discretized QFTs, depends on the remaining UV cutoff and shows that the bare mass must be fine tuned to a cutoff-dependent value to yield the physical mass of the low-energy excitations.
In order to address this point, we apply large-N techniques away from half filling via the introduction of a chemical potentialμ in the GNW model. Following the orthodox prescription for Euclidean LFTs [83], the hopping term κ ν (1 − sγ ν ) of the Euclidean action S E W (36) is modified to e sμδ ν0 κ ν (1 − sγ ν ), such that time-like hopping is promoted in the forwards direction by a factor eμ , and suppressed by e −μ when hopping backwards. As a consequence, one can study the phase diagram of the GNW model at finite densities by solving the gap equations (40)-(41) with the sum over over the time-like momenta now given by k 0 = 2π(n τ + 1/2)/N τ − iμ.
Moreover, using the Euclidean partition function, one finds that the conserved fermion charge density n q = − ∂ ln Z ∂μ is Settingμ → 0, this quantity becomes proportional to the expectation value of the time-like component of the vector , which is the discretized version of the continuum vector cur- : for Wilson fermions. Therefore, the time-like component is simply related to the fermion density in the continuum limit, and we can readily explore situations away from half-filling n q = 0. Interestingly, while the gap equations (40)- (41) remain symmetrical under the transformation (28) using the renormalized mass (54), the charge density n q (m,g 2 ,μ) has only an approximate symmetry.
We now solve the gap equations (40)-(41) with a dimensional chemical potential µ ≡ ξμ = 0, which yield the phase diagram of Fig. 5(a), where the axes have been rescaled to match those of Fig. 4. We see that, as a consequence of the non-zero-chemical potential, the leftmost and rightmost cusps of the half-filled phase diagram of Fig. 4 split into a couple of cusps each, such that the region hosting the Aoki phase becomes smaller. The decrease of the Aoki phase can be qualitatively understood as follows. For µ > 0, one expects that the charge density n q will eventually rise from the value n q = 0 characterizing the half-filled regime. As a consequence, a Fermi surface will be formed in a certain parameter regime, which consists of two disconnected Fermi points in 1+1d. This has the effect of disfavoring the particle-anti-particle pairing required for the pseudoscalar condensate Π(x) = 0, since the excitation energy for such a zero-momentum pair would be on the order of ∆ε ∼ 2µ. Accordingly, the Aoki phase should shrink as the chemical potential is increased. This is corroborated by the curves for n q (m r ) in Fig. 5(b), which rise above zero only around the borders of the dropletshaped region between the two newly-formed cusps for µ > 0 (see the grey regions in Fig. 5(a)). These are precisely the regions where the half-filled Aoki phase has been expelled from. We also note that Fig. 5(b) shows an approximate symmetry about −m r a = 1 which we expect to become exact for ξ → ∞.
Let us now describe how these results can be used to determine, in a controlled way, the extent of the Aoki phase at half-filling. The empty circles of Fig. 5(a) mark the so-called onset, beyond which the ground state has a non-zero charge density (i.e. n q > 0 for µ > µ o (g 2 )). By numerically obtain- ing such onsets for a variety of parameters and time discretizations, we obtain Fig. 6. We note that the variation for finite ξ is probably due to non-universal effects since in the sums over the Brillouin zone of Eqs. (40)- (40), as the chemical potential enters as e.g. sinh(ξ µa). One observes from this figure that all curves come closer in the limit µ → 0, and seem to approach a limit as ξ → ∞. This limiting value corresponds to the point where half-filled Aoki phase terminates, proving that these phase does not extend all the way down to the weak coupling limit, but only survives down to g 2 (µ o = 0) ≈ 0.8 according to the results of Fig. 6. Let us note that extracting this limiting value is numerically hard; for instance, the curvature appears slightly sensitive to temperature, as revealed by calculations on 512 × 1024ξ lattices. However, our approximate prediction g 2 (µ o = 0) ≈ 0.8 is consistent with the cusps of Fig. 4, where the Aoki phase terminates. Fig. 6 therefore strengthens our belief in the existence of a tricritical point at non-zero g 2 (µ o = 0); for couplings below this value there is a direct transition between trivial and topological insulating phases as m is varied, and no parity breaking Aoki phase is encountered in the middle.
So far, we have used the large-N results for a non-zero chemical potential to extract features of the half-filled phase diagram by taking the limit µ → 0 in a controlled manner. However, we note that another interesting question would be to study the fate of the symmetry-protected topological phases, and the appearance of other new phases of matter, in the GNW model away from half filling. In that respect, we note that our large-N results point towards the appearance of a new phase (i.e. droplet-shape region of Fig. 5). Since we have argued that the finite charge densities appear due to the formation of a Fermi surface, it is reasonable to expect that such densities will not drop abruptly to zero as we move away from the critical line. In that sense, the droplet-shaped region may either correspond to a metallic phase where the Fermi points occur at different momenta as the microscopic parameters are modified, or maybe to a kind of charge-density-wave where the fermionic density forms a regular periodic pattern. Un-derstanding the nature of this phase lies outside of the scope of the present work, and will be the subject of a future work. We advance at this point that the density-matrix renormalization group methods discussed in the following section could be adapted to study situations away from half filling, and are a potential tool to address the nature of this new phase. Moreover, we also note that the sign problem for µ = 0 can be safely avoided for any discretized Gross-Neveu or Nambu-Jonal-Lasinio models, such that Monte Carlo techniques [84] could also be applied to the present problem, and extensions thereof.

B. Large-N benchmark via matrix product states
In this section, we test the above large-N prediction for the single-favor GNW lattice model using numerical routines based on matrix product states (MPS) [14] (i.e. a variational version of real-space density-matrix renormalization group method [85]). On the one hand, this can be considered as the most stringent test of the validity of the large-N approach, as we are indeed very far from the large-N limit. On the other hand, the choice of N = 1 is also motivated by the fact that the single-flavor GNW model can be realized in cold-atom experiments following the scheme of Sec. II C described below. Note that the N = 1 flavor of the continuum Gross-Neveu QFT (2) with an additional mass term corresponds to the socalled massive Thirring model [86]. The discretization of this QFT using the Wilson approach allows us to discuss the occurrence of symmetry-protected topological phases in this LFT, and use it to benchmark the large-N predictions for the phase diagram of the GNW model with a finite number of flavors.

High-energy physics to condensed matter mapping
We consider the GNW lattice Hamiltonian (3)-(5) for a single fermion flavor N = 1. By performing a U(1) gauge transformation to the spinors Ψ(x) → e −i π 2a x Ψ(x), which can be understood as an instance of a Kawamoto-Smit phase rotation in LFTs [87], and using the algebraic properties of the gamma matrices, we can rewriteH W →H W = a∑ x∈Λ :H W :, wherẽ In this notation, the Hamiltonian looks similar to the Hubbard model [88], a paradigm of strongly-correlated electrons in condensed matter [89], with an additional spin-orbit coupling. Note that this formulation only differs from Eqs. (3)-(5) on the particular distribution of the complex tunnellings, which can be understood as a gauge transformation on a background magnetic field maintaining an overall π-flux. Indeed, the defining property of the above Kawamoto-Smit phases is that they yield a π-flux through an elementary plaquette. In order to understand the origin of this magnetic flux, let us introduce the following notation for the Dirac spinor Here, the dimensionless fermion operators c i, depend on a spinor index ∈ {u, d} that can be interpreted in terms of the upper ( = u) and lower ( ∈ d) legs of a synthetic ladder, and i ∈ {1, · · · , N s } labels the positions of the rungs of the ladder (see Fig. 7(a)). Considering our particular choice of gamma matrices γ 5 = σ x , γ 0 = σ z , the corresponding HamiltonianH W for the chosen Wilson parameter r = 1 can be rewritten as where we have introduced s = +1 (s = −1) for the upper = u (lower = d) leg of the ladder, and¯ = d, u for = u, d. As can be seen in Fig. 7(a), there is a net π-flux due to an Aharonov-Bohm phase that the fermion would pick when tunneling around an elementary plaquette.
In particular, Eq. (57) can be understood as a generalized Hubbard model on a ladder corresponding to the imbalanced Creutz-Hubbard model [71], which is an interacting version of the so-called Creutz ladder [90,91]. The first line in Eq. (57) describes the horizontal and diagonal tunneling of fermions with strengtht = 1/2a, which are subjected to an external magnetic π-flux threading the ladder (see Fig. 7(a)). One thus finds that the UV cutoff of the GNW model Λ = 1/a is provided by the maximum energy within the band structure Λ = 2t of the Creutz-Hubbard model. Likewise, one understands that the first term in the second line of Eq. (57) corresponds to an energy imbalance between both legs of the ladder ∆ε/2 = (m + 1/a), and yields a single-particle Hamiltonian in momentum space that is similar to Eq. (6), namely Finally, the last term of Eq. (57) amounts to a Hubbard-type density-density interaction V v n n,u n n,d between fermions residing on the same rung of the ladder, which repel themselves with a strength V v = g 2 /a. According to this discussion, the high-energy-physics GNW lattice model is gauge equivalent to the condensedmatter imbalanced Creutz-Hubbard model. Similarly to the high-energy physics convention of working with dimensionless parameters m/Λ = ma and g 2 , the condensed-matter community normalizes the couplings to the tunneling strengtht, such that the exact relation between the microscopic parameters of these two models is Let us also note that, in the condensed-matter context, the lattice constant d of the model (57) is fixed by the underlying Bravais lattice of the solid, which is typically set to d = 1 in the calculations (58) (i.e. lattice units). Note, however, that this does not preclude us from taking the continuum limit. In this case, the continuum limit corresponds to the low-energy limit, wheret = 1/2a (i.e. UV cutoff) is much larger than the energy scales of interest. By setting the model parameters in the vicinity of a second-order quantum phase transition, the relevant length scales fulfill ξ l d, and one recovers universal features that are independent of the microscopic lattice details, and can be described by a continuum QFT. In this section, we exploit the above mapping (59) to explore the phase diagram of the N = 1 lattice GNW model by importing some of the condensed-matter and quantuminformation techniques described in [71]. In particular, we will use the numerical matrix-product-state results to benchmark the large-N predictions. We remark that this mapping also becomes very useful in the reverse direction, as certain aspects of the Creutz-Hubbard model become clarified from the high-energy perspective of the GNW model.
In the parameter regime ∆ε/4t < 1, which corresponds to a bare mass −ma ∈ [0, 1], we found that the imbalanced Creutz-Hubbard model displays three distinct phases: an orbital paramagnet, an orbital ferromagnet, and an SPT phase [71]. The orbital paramagnet corresponds to a gapped phase of matter that is characterized by the absence of long-range order and any topological feature. Therefore, this phase should correspond to the trivial band insulator of the GNW model in Fig. 1.
The orbital ferromagnet, on the other hand, is a phase displaying an Ising-type long-range order due to the spontaneous breaking of a discrete orbital symmetry. Accordingly, it should correspond to the parity-broken Aoki phase of the Gross Neveu model in Fig. 1. To show this correspondence in more detail, let us comment on the orbital magnetization introduced for the Creutz-Hubbard ladder T 0 = T y i = 1 2 ic † i,u c i,d + c.c. = 0 ∀i, and show that it is related to an order parameter of the GNW model. The parity symmetry of the GNW model that is broken in the Aoki phase, namely Ψ(x) → ηI N ⊗ γ 0 Ψ(−x) with |η| 2 = 1, corresponds to c i,u → ηc N s −i,u , and c i,d → −ηc N s −i,d in the Creutz-Hubbard ladder. Hence, one finds that T 0 = T y i → − T y N s −i = −T 0 is spontaneously broken by the orbital ferromagnet. We thus see that, in the language of the synthetic Creutz-Hubbard ladder, the pseudoscalar condensate corresponds to an Ising-type ferromagnet with a non-zero orbital magnetization T 0 = −Π 0 a/2 = − Ψiγ 5 Ψ a/2. This connection also teaches us that one can perform a rigorous finite-size scaling of the pseudoscalar condensate to obtain accurate predictions of the critical lines enclosing the whole Aoki phase, instead of using the various mappings discussed in [71].
Finally, as shown explicitly in [71], the Creutz-Hubbard ladder also hosts a correlated SPT phase, which displays a double-degenerate entanglement spectrum [92] due to a couple of zero-energy edge modes. This phase should thus corresponds to the BDI symmetry-protected topological phase of the GNW model discussed throughout this work (see Fig. 1). Let us remark, however, that the topological insulator of the Creutz-Hubbard model lies in the symmetry class AIII, breaking explicitly the time-reversal T and charge-conjugation C symmetries, yet maintaining the sublattice symmetry. According to our discussion below Eq. (12), we see that the Creutz-Hubbard single-particle Hamiltonian breaks T: γ 0 (h CH −k ) * γ 0 = h CH k , and C: γ 5 (h CH −k ) * γ 5 = −h CH k , explicitly. On the other hand, the combination S = TC yields (γ 1 ) † h CH k γ 1 = −h CH k , such that the topological insulator of the Creutz-Hubbard ladder is in the AIII symmetry class. Therefore, the last element of our high-energy physics to condensed-matter dictionary is the mapping between the symmetry classes BDI ↔ AIII, which is a direct consequence of the above local gauge transformation/Kawamoto-Smit phase rotation. Although differences will arise regarding perturbations that explicitly break/preserve the corresponding symmetries (e.g. disorder), the phase diagram of the translationallyinvariant GNW model should coincide exactly with that of the Creutz-Hubbard ladder provided that one uses the relation between microscopic parameters in Eq. (59).
With this interesting dictionary for the correspondence of phases, and the microscopic parameter mapping in Eq. (59), we can use numerical matrix-product-state simulations, extending the parameter regime studied in [71] from ∆ε/4t < 1 to −1 < ∆ε/4t < 1. In this way, we can explore the full phase diagram diagram of the N = 1 GNW model, and compare it to our previous large-N predictions for −ma ∈ [0, 2]. Let us recall that the large-N approach fulfills (28), such that the obtained phase diagrams have a mirror symmetry about −ma = 1. However, it is not clear a priori if this symmetry is a property of the model, or if it is instead rooted in the approximations of the large-N prediction. We will be able to address this question with our new matrix-product-state simulations.
In Figs. 8 (a)-(b), we discuss a representative example of the finite-size scaling of the pseudoscalar condensate Π 0 = Ψiγ 5 Ψ for the transition between the trivial, or topological, band insulators and the Aoki phase. One clearly sees that the matrix-product-state numerical simulations for different lengths display a crossing that gives access to the critical point (main panel of Figs. 8 (a)-(b)), and that the data collapse of (inset of Figs. 8 (a)-(b)) corroborates that this critical point lies in the Ising universality class. Note that the pseudoscalar condensate gives no information about the phase transition between the trivial and topological insulators. In order to access this information, the mapping to the Creutz-Hubbard ladder becomes very useful, as it points to the possibility of using a generalized susceptibility associated to the variation of the scalar condensate σ 0 = ΨΨ with the bare mass χ σ 0 = ∂ σ 0 /∂ m. As shown in the main panel of Fig. 8 (c), this susceptibility diverges at the critical point of the thermodynamic limit, and can be used to perform a finite-size scaling.
Repeating this procedure for various critical points, we obtain the red empty circles displayed in Fig. 9, which are compared with the large-N results of Fig. 1 represented as solid lines. We can thus conclude that the large-N predictions are qualitatively correct, as they predict the same three possible phases, and the shape of the critical lines is qualitatively similar to the matrix-product-state prediction. Moreover, the agreement between the critical lines becomes quantitatively correct in the weak-coupling limit g 2 , |ma| 1, which is the regime where the asymptotically-free Gross-Neveu QFT (1) is expected to emerge from the lattice model. Since both the mass and the interaction strengths are relevant perturbations growing with the renormalization-group transformations, one expects that a continuum limit with physical parameters well below the UV cutoff can be recovered provided that g 2 , |ma| 1. Let us also remark that the matrix-product-state simulations are consistent with the mirror symmetry about −ma = 1 of the large-N gap equations. Therefore, it seems that this symmetry is an intrinsic property of the GNW model, which is easy to understand in the non-interacting limit, but not so obvious in the interacting case. On general grounds, Fig. 9 shows that the large-N prediction tends to overestimate the extent of the Aoki phase, predicting that the spontaneous breaking of the parity symmetry occurs for weaker interactions and smaller masses. This trend could be improved by considering next-to-leading-order (NLO) corrections to the saddle-point solution, and will be the subject of a future study. In this sense, our results suggest that large-N methods from a high-energy context can be a useful and systematic tool to study problems of correlated symmetry-protected topological phases in condensed matter.
We now comment on further interesting features that can be learned from this dictionary, and imported from condensed matter into the high-energy physics context. In Fig. 9, we have highlighted with a semi-transparent orange star the The red circles represent the critical points of the N = 1 Gross Neveu lattice model obtained by with matrix product states. The semitransparent green lines joining these points delimit the trivial band insulator, Aoki phase, and the BDI symmetry-protected topological phase. Note again that this SPT phase corresponds to the AIII topological insulator of the Creutz-Hubbard model .These lines are labelled by N = 1, and by the central charge c ∈ {1/2, 1} of the conformal field theory that controls the continuum QFT at criticality. These results serve to benchmark the large-N predictions for the critical lines, which are represented by dark solid lines, and labelled with "N = ∞". We also include the exact critical point at (−ma, g 2 ) = (1, 4), which is depicted by an orange star, and the strong-coupling critical lines that become exact in the limit of g 2 → ∞, which are depicted by dashed orange lines. The matrixproduct-states predictions match remarkably well these exact results.
critical point separating the topological and Aoki phases at (−ma, g 2 ) = (1, 4). This point corresponds to a Creutz-Hubbard model with vanishing imbalance ∆ε = 0, and strong repulsion V v = 8t. Interestingly, it is precisely at this point that an exact quantum phase transition is found by mapping the Creutz-Hubbard ladder onto an exactly-solvable quantum impurity Ising-type model via the so-called maximally-localized Wannier functions [71]. In this way, one learns that the lattice GNW model can be solved exactly for a particular limit with relatively strong couplings, and that the corresponding quantum phase transition must lie in the Ising class. From a high-energy perspective, the whole critical line separating the topological and Aoki phases should be controlled by the continuum QFT of a Majorana fermion, and not by the standard Dirac-fermion QFT expected at weak couplings (i.e. along the critical line separating the topological and trivial insulators). We have proved this rigorously using the numerical scaling of the entanglement entropy [93], which shows that this criti- with the massless Dirac fermion (see Fig. 10).
Let us now discuss the orange dashed line of Fig. 9, which describes an exact solution that becomes valid in the strong-coupling limit g 2 1. From the parameter correspondence (59), this regime corresponds to the strongly-interacting Hubbard model, where one expects to find super-exchange interactions between the fermions [94]. In this case, these superexchange can be described in terms of an orbital Ising model with ferromagnetic coupling J = −2/g 2 a, and subjected to a transverse magnetic field B = 2(m + 1/a). According to the exact solution of the transverse Ising model [95], the strongcoupling critical line J = 2B corresponds to g 2 = 1/2(ma+1). This line, and its mirror image, have been depicted by the orange dashed lines of Fig. 9, and shows a very good agreement with the numerical critical points of the GNW model at strong-couplings g 2 1. Since the strong-coupling mapping yields a transverse Ising model, we learn again that the corresponding continuum QFT at criticality is that of a Majorana fermion, which is corroborated again by the matrix-productstate scaling of the entanglement entropy yielding a central charge of central charge c ≈ 1/2 (see Fig. 10). Therefore, the condensed-matter mapping teaches us that the GNW lattice model has an exact solution in the strong-coupling limit, and both critical lines delimiting the Aoki phase lead to a continuum limit controlled by a Majorana-fermion QFT. These results show that condensed-matter methods can offer a useful and systematic tool to benchmark large-N methods applied to problems of asymptotically-free LFTs in a high-energy context. In future works, we will study leading order 1/N correc-tions to the present large-N approach, and see how fast they approach the exact and quasi-exact results for the phase diagram discussed in this section.
First of all, the bare tunneling must be inhibited by the gradient t ∆. Then, the inter-leg tunnellings of Fig. 7(b) (crossed black lines) can be laser-assisted by a Raman pair [98], which also leads to the energy imbalance terms (yellow loops). We set (i) the Raman frequencies to ∆ω 1 = (ε ↑ − ε ↓ ) + ∆ + ∆ε/2, and ∆ω 2 = (ε ↓ − ε ↑ ) + ∆ − ∆ε/2, where ∆ε is small detuning, (ii) the two-photon Rabi frequencies (phases) to Ω 1 = Ω 2 =: Ω (ϕ 1 = ϕ 2 =: ϕ), and (iii) the corresponding Raman wave-vectors to ∆k k k 1 · e x = ∆k k k 2 · e x = 0. In a rotating frame, the Raman-assisted tunneling arising from the corresponding v i, j σ ,σ (t) term contributes with which contains precisely the desired crossed tunnellings for ϕ = π/2, and the energy imbalance of Fig. 7(b). In order to engineer the horizontal tunneling of Fig. 7(b) (green lines), we shall make use of a third Raman pair, but this time far detuned from the atomic transition ∆ω 3 (ε ↑ − ε ↓ ). In this situation, when the corresponding laser intensities are weak, the Raman term leads to a crossed-beam ac-Stark shift that can be interpreted as slowly-moving shallow optical lattice that acts as a periodic modulation of the on-site energies where Ω σ is the two-photon ac-Stark shift for each of the hyperfine levels, which can be controlled by tuning the intensity and polarization of the lasers. We set (i) the Raman frequency in resonance with the gradient ∆ω 3 = ∆ (ε ↑ − ε ↓ ); (ii) the Raman wave-vector to ∆k k k 3 · e x = k L,x with respect to the static optical lattice; and (iii) the Raman phase ϕ 3 = 0. In a rotating frame, the atoms can absorb energy from this shallow moving lattice, such that the horizontal tunneling gets reactivated [99,100], according to where we have introduced the n-th order Bessel function of the first class J n (x). According to this expression, we can laser-assist the horizontal hopping with the desired signs of Fig. 7(b) by exploiting the state-dependence of the dressed tunneling rates, and setting This can be achieved, while simultaneously maximizing the dressed tunneling, by setting Ω ↑ = 3∆ ≈ 0.6Ω ↓ . Let us note that the cross-tunneling (62) will also get a multiplicative renormalization due to this periodic modulation (63), which will be proportional to J 0 ((Ω ↑ + Ω ↓ )/∆). This dressing is similar to the effect exploited for the so-called coherent destruction of tunneling [101]. To achieve the relation of the tunnellings of Fig. 7(b), one should modify the Rabi frequency of the Raman beams Ω, such that although we note that there might be other strategies to fulfill both constraints (65)-(66) simultaneously. Altogether, considering also Hubbard interactions, the correspondence between the cold-atom and the Gross-Neveu parameters is As announced at the beginning of this section, this scheme provides a slight simplification over the proposal for the Creutz-Hubbard model [71], which required the use of an intensity-modulated superlattice, instead of the shallow moving lattice (63) already implemented in experiments [100]. At this point, we comment on an interesting alternative that would simplify considerably the cold-atom scheme. As realized recently [73], a different choice of the gamma matrices γ 0 = σ x , γ 5 = σ y , simplifies considerably the tunneling of Eq. (60), since iγ 5 + rγ 0 = 2σ + for a Wilson parameter r = 1, where we have introduced the raising operator σ + = |↑ ↓|. Accordingly, the kinetic energy of the Wilson fermions can be depicted by the scheme of Fig. 7(c). Let us note that the BDI symmetry class can be readily understood by realizing that the synthetic ladder of this figure can be deformed into a single chain with dimerized tunnellings, and thus corresponds to the Su-Schrieffer-Hegger BDI topological insulator [102].
This representation was exploited in [73] to propose a coldatom realization of quantum electrodynamics with Wilson fermions in (1 + 1) dimensions (i.e. Schwinger model). In that case, one should introduce a bosonic species to simulate the gauge field, and exploit the spin-changing boson-fermion atomic scattering to obtain the gauge-invariant tunneling of the lattice gauge theory. In our case, the required experimental tools are already contained in our previous description and, more importantly, can be considerably simplified with respect to the above discussion. The vertical tunnellings of Fig. 7(c) can be obtained from a Raman pair with ∆ω 1 = (ε ↑ − ε ↓ ), whereas the diagonal tunnellings require another pair of Raman beams with ∆ω 2 = (ε ↓ − ε ↑ ) + ∆, but no additional periodic modulations would be required. Therefore, if no additional disorder is to be considered, which could depend on the particular symmetry class and choice of gamma matrices, this later approach should be followed for the cold-atom experiment, as it simplifies the experimental requirements for the quantum simulation of the GNW model.
Let us finally comment on another interesting alternative. The non-interacting Creutz ladder has been recently realized in multi-orbital optical-lattice experiments [103] that exploit two orbital states of the optical lattice to encode the legs of the ladder, and orbital-changing Raman transitions to implement the inter-leg tunnelings. It would be interesting to study the type of multi-orbital interactions [104] that can be generated in this setup, and the possibility of simulating directly the GNW model studied in this work.

IV. CONCLUSIONS AND OUTLOOK
In this work, we have described the existence of correlated symmetry-protected topological phases in a discretized version of the Gross-Neveu model. We have applied large-N techniques borrowed from high-energy physics, complemented with the study of topological invariants from condensed matter, to unveil a rich phase diagram that contains a wide region hosting a BDI topological insulator. This region extends to considerably strong interactions, and must thus correspond to a strongly-correlated symmetry-protected topological phase. We have shown that this phase, and the underlying topological invariant, can be understood in terms of the renormalization of Wilson masses due to interactions (i.e. dynamic mass generation due to a scalar fermionic condensate). This renormalization has been used to find a critical line at weak couplings that separates the topological insulator from a gapped phase that can be adiabatically deformed into a trivial product state (i.e. trivial band insulator). Moreover, we have shown that for sufficiently-strong interactions, a gapped phase where parity symmetry is spontaneously broken (i.e. Aoki phase) is formed due to the appearance of a pseudoscalar fermion condensate. The large-N prediction has allowed us to find the critical line separating the topological insulator from the Aoki phase by studying the onset of the pseudoscalar condensate, and show that it terminates at a tricritical point where all these three phases of matter coexist.
By using both Hamiltonian and Euclidean lattice approaches, we have been able to pinpoint important details that must be carefully considered when taking the time-continuum limit of the lattice approaches, such that standard methods of lattice field theories can be used to describe quantitatively the phase diagram of the Gross-Neveu-Wilson Hamiltonian. In particular, we have described how lattice artifacts can appear in the standard dimensionless formulation of the Euclidean field theory, and how the spurious time doublers, even when residing at the cutoff of the theory, can renormalize the bare parameters and introduce qualitative modifications to the layout of the phases. The results hereby presented will serve as the starting point for the application of other well-established Euclidean lattice techniques to explore the phenomenology of leading-order corrections that appear for finite N.
Motivated by the possibility of implementing a cold-atom quantum simulator of the Gross-Neveu-Wilson model for a single flavor N = 1, which has also been described in this work, we have benchmarked these large-N predictions by means of quasi-exact numerical methods based on matrix product states. In particular, we have shown that the single-flavor model, corresponding to a discretized version of the massive Thirring model, can also be mapped into a condensed-matter Hamiltonian of spinless fermions hopping on a two-leg ladder, and interacting via Hubbard-type couplings. This connection has allowed us to identify the phases of the Gross-Neveu-Wilson model, discussed above, with condensed-matter counterparts that include orbital paramagnets and ferromagnets, as well as a chiral-unitary topological phase. In this way, the matrix-product-state simulations can readily access a variety of observables to determine the position of the critical lines, which show a remarkable qualitative agreement with the large-N predictions that becomes even quantitative in the region where the continuum Gross-Neveu QFT is expected to emerge (i.e. weak couplings). These numerical simulations also prove that the symmetry of the large-N phase diagram holds for N = 1, and should then be maintained at all orders O(1/N α ). Beyond the matrix-productstate simulations, the aforementioned mapping has allowed us to import exact results for the Gross-Neveu-Wilson model in the regime of intermediate and strong couplings, which originate from quantum-impurity and quantum magnetism techniques in condensed matter.
Therefore, we believe that our work constitutes an example of the useful dialogue and exchange of ideas between the high-energy physics, condensed-matter, quantuminformation, and quantum optics communities, stimulating further cross-disciplinary efforts in the future. As an outlook, one can easily foresee that lattice field-theory techniques to study leading-order corrections to the large-N behavior will be very useful to elucidate the mechanism that induces strong correlations in the symmetry-protected topological phase of the Gross-Neveu-Wilson model. Likewise, quantum-information approaches might be useful to understand the entanglement content of those phases, making a connection to the lattice field-theory techniques. As already pointed out by the Euclidean lattice results, new phases of the Gross-Neveu-Wilson model can arise as one moves away from half-filling. It will be very interesting to explore the nature of these phases using some of the high-energy and condensed-matter techniques hereby discussed. We also note that the techniques hereby presented can be generalized to other lattice Hubbard-type models, not necessarily connected to well-known relativistic QFTs. In particular, it will be very interesting to apply them to the study of higher-dimensional models hosting topological phases of matter. In this context, Aoki phases have been identified in the limit of very-strong Coulomb interactions via strong-coupling techniques of lattice gauge theories [37,77]. These results have been used to conjecture the qualitative shape of the phase diagram in the regime of weak to intermediate inetractions [77], which is expected to be more relevant for the understanding of correlation effects in topological insulating materials.