Lattice regularisation and entanglement structure of the Gross-Neveu model

We construct a Hamiltonian lattice regularisation of the $N$-flavour Gross-Neveu model that manifestly respects the full $\mathsf{O}(2N)$ symmetry, preventing the appearance of any unwanted marginal perturbations to the quantum field theory. In the context of this lattice model, the dynamical mass generation is intimately related to the Coleman-Mermin-Wagner and Lieb-Schultz-Mattis theorem. In particular, the model can be interpreted as lying at the first order phase transition line between a trivial and symmetry-protected topological (SPT) phase, which explains the degeneracy of the elementary kink excitations. We show that our Hamiltonian model can be solved analytically in the large $N$ limit, producing the correct expression for the mass gap. Furthermore, we perform extensive numerical matrix product state simulations for $N=2$, thereby recovering the emergent Lorentz symmetry and the proper non-perturbative mass gap scaling in the continuum limit. Finally, our simulations also reveal how the continuum limit manifests itself in the entanglement spectrum. As expected from conformal field theory we find two conformal towers, one tower spanned by the linear representations of $\mathsf{O}(4)$, corresponding to the trivial phase, and the other by the projective (i.e. spinor) representations, corresponding to the SPT phase.


Introduction
Lattice field theory, and in particular, lattice gauge theory, has been among the most successful techniques to probe the non-perturbative behaviour of quantum field theories (QFTs), such as those appearing in the standard model. The accurate determination of the proton and neutron masses has been one of the most noteworthy triumphs resulting from this effort. The default approach is to apply Monte Carlo sampling to the path integral in discretised Euclidean spacetime [1][2][3][4][5].
In recent years, the use of tensor network methods has been proposed as an alternative [6], with the promise that these are able to access dynamical information and do not suffer from sign problems in the case of fermionic densities or far-from-equilibrium situations [7]. One can apply tensor renormalisation group techniques as an alternative to Monte Carlo sampling to the path integral in discretised spacetime [8][9][10][11][12][13][14][15][16][17][18]. Alternatively, one can target the wave functional using a tensor network ansatz and apply variational techniques using the field theory hamiltonian (where only the spatial dimensions are discretised) . This approach is also closely related to the various experiments and proposals for the analog or digital quantum simulation of lattice field theory using various platforms such as trapped ions, superconducting circuits or cold atoms in optical lattices (see Ref. [6] and references therein). Aside from preliminary explorations of Z 2 and U(1) gauge theories in (2+1) dimensions [21,26,28,34,43], most of the tensor network effort has so far been invested in QFTs in (1+1) dimensions, and in particular the λφ 4 model [8,9,17,18,20,22,24] and the Schwinger model, i.e. (1+1)-dimensional quantum electrodynamics, (as well as nonabelian generalizations thereof) [19, 23, 25, 27, 29-31, 33, 35-39, 41, 44]. These models are superrenormalizable, meaning that the coupling constant has a positive mass dimension and sets the energy scale. The relation between the lattice and continuum parameters is governed by a limited number of divergent diagrams, and observables converge like power laws in the lattice spacing a as the continuum limit is approached.
In this manuscript, we use lattice field theory and tensor network tools (numerical and analytical) to probe the non-perturbative properties of the Gross-Neveu (GN) model [45], a (1+1)-dimensional model of N massless but interacting fermion flavours, which shares several non-perturbative features with (3+1)-dimensional quantum chromodynamics (see Refs. [13,42] for tensor network studies of the closely related Thirring model, an integrable model for a single massive interacting fermion flavour). The GN interaction has a discrete chiral symmetry and is marginally relevant, (i.e. renormalisable and asymptotically free). The interaction term leads to spontaneous breaking of this chiral symmetry and, associated with this, dynamical mass generation. Here observables converge logarithmically slow as the continuum limit is reached. This increases the importance of symmetries prohibiting the presence of other marginally relevant perturbations that could spoil the already slow convergence. Indeed, it turns out to be crucial to meticulously construct the lattice Hamiltonian so as to maximally preserve the symmetries of the field theory, in order to reliably obtain the continuum limit.
Being a paradigmatic model, the GN model has been the subject of several numerical and theoretical studies. Theoretical studies have focused on determining the scattering matrix and full excitation spectrum [46][47][48], as well as a precise determination of the mass gap [49][50][51][52][53][54] using a variety of techniques, including thermodynamic Bethe ansatz, large N expansions and the variationally optimised renormalisation group. Most numerical lattice studies use Monte Carlo techniques on the Euclidean lattice, where the fermions are dealt with by replacing the four point interaction by a coupling to an auxiliary bosonic field (using a Hubbard-Stratonovich transformation) and integrating out the resulting quadratic fermion terms, leaving the calculation of the fermion determinant as a computational problem [55][56][57]. In particular, there has been interest in the phase diagram at finite temperature and chemical potential, and the possible existence of an inhomogeneous phase [58,59].
A lattice prescription of the kinetic term of the fermion model can be obtained using the Wilson prescription [60] or using the staggered formulation of Kogut and Susskind [61,62]. The former explicitly breaks the chiral symmetry, resulting in additive mass corrections that need to be compensated by a properly tuned bare mass term, in order to reach the continuum limit. Furthermore, the Wilson prescription also leads to Aoki phases [63,64], where reflection (parity) symmetry is broken and a pseudoscalar condensate is formed. Triggered by interest from the optical lattice community, the phase diagrams of this 'Gross-Neveu-Wilson' lattice model and its chiral extension in the limits N → ∞ and N = 1 were studied in recent publications [65,66], and feature both trivial, topological and symmetry broken Aoki phases.
The staggered formulation, on the other hand, exhibits remnant lattice symmetries which prohibit perturbative mass corrections. For the particular case of a lattice Hamiltonian (i.e. continuous time) in (1+1) dimension, this remnant symmetry corresponds to full translation invariance of the staggered model [62]. Spontaneous breaking of discrete chiral symmetry can then be related to Peierls dimerisation, so that the GN model with N = 2 also arises as a continuum description of polyacetylene. In particular, the GN model provides a good description of the resulting topological soliton (kink) that interpolates between the two ground states [67,68] and which is traditionally described as an explicit domain wall in the Su-Schrieffer-Heeger (SSH) model [69]. While the GN model as low-energy field theory in the context of the SSH model is well documented [70], the reverse direction where the SSH model is used as inspiration to construct a precise lattice regularisation of the GN model was, to the best of our knowledge, not considered.
The outline of this paper is as follows. Section 2 summarises the field theoretic description of the model. In Section 3 we construct the lattice model and discuss the dynamical mass generation and kink degeneracy from a condensed matter perspective. In section 4, a large-N mean-field solution is given and found to be consistent with the large N field theory. Section 5 introduces the symmetric uniform matrix product state (MPS) ansatz, which is then used in section 6 to numerically probe the low-energy behaviour of the model for N = 2. In section 7, we discuss the continuum limit from the point of view of entanglement. Finally, section 8 provides a concluding discussion and outlook.

Gross-Neveu model in a nutshell
We first provide a short introduction to the GN model and its symmetries before porting it to the lattice. The Lagrangian density for the GN field theory reads where c labels the N different flavours or colours of fermions, ψ c is the two-component Dirac spinor for flavour c, / ∂ = γ 0 ∂ 0 + γ 1 ∂ 1 andψ c = ψ † c γ 0 , where {γ µ , γ ν } = 2g µν and g µν is the inverse metric tensor.
The model has an obvious SU(N ) flavour mixing symmetry that can be extended to an O(2N ) symmetry, which also includes the total U(1) particle number symmetry and the charge conjugation symmetry (which relates the two disconnected components of O(2N )). This O(2N ) symmetry is made explicit by rewriting the Dirac spinor in terms of its Majorana components. By choosing a specific set of gamma matrices where both γ 0 and γ 1 are strictly imaginary (so that β = γ 0 is imaginary and thus antisymmetric, whereas α = γ 0 γ 1 is real symmetric), the Majorana components correspond to the real and imaginary components of the Dirac spinor, i.e. ψ c = (λ 2c−1 + iλ 2c )/ √ 2 and thus This formulation shows the explicit invariance under λ m → O mn λ n for O ∈ O(2N ). Coleman's theorem for relativistic theories [71], related to the Mermin-Wagner theorem in condensed matter or statistical physics [72], guarantees that this continuous symmetry cannot be broken and is thus present in the spectrum of the theory.
The GN model has an additional Z 2 chiral symmetry that acts as ψ → γ 5 ψ and prohibits perturbative contributions to the condensate σ = c∈N ψ c ψ c , or thus, a perturbative mass term. Nonetheless, this Z 2 symmetry is spontaneously broken and gives rise to a non-perturbative mass scale. The effect that the ground state exhibits a dimensionful condensate, despite the absence of dimensionful parameters, other than the ultraviolet (UV) regulator scale, is known as dimensional transmutation. The required renormalization group (RG) invariant mass scale can be obtained from the β function where µ is the regulator scale. This mass scale Λ is only an infrared (IR) scale for asymptotically free theories (β 0 < 0). Furthermore, the mechanism by which it enters the IR theory (if any) must necessarily be non-perturbative. For the O(2N )-symmetric GN model, the first terms in the β function were calculated to be β 0 = − N −1 2π and β 1 = N −1 4π 2 [73,74]. The condensation of σ gives rise to a rich spectrum of massive particles [46,48,75,76]. Given the O(2N ) symmetry, some understanding of the representation theory of the corresponding Lie-algebra so(2N ) is useful to label this spectrum. There are two distinct fundamental (half) spin representations of dimension 2 N −1 , which are transformed into each other by conjugation or by application of a reflection element from O(2N ) (determinant −1). The other fundamental representations r = 1, . . . , N − 2 are tensor representations, with r = 1 the defining (vector) representation. We also refer to the spinor representations as projective representations, which generalise the concept of 'representations up to a phase' -as opposed to linear representations such as the tensor representations-to arbitrary groups.
The spectrum of the GN model contains both trivial and topological excitations, i.e. kinks that interpolate between the two vacuum states. Unlike in conventional (i.e. Isingtype) Z 2 symmetry breaking, where the kink from one vacuum to the other is unique, in the case of GN the kinks are of the Callen-Coleman-Gross-Zee type [77] and bind massless fermions. They transform according to the fundamental spinor representations [47]. This is similar to Jackiw-Rebbi kinks [78] and we will interpret this from a condensed matter perspective as the protected gapless edge modes on the interface between a trivial and SPT phase, when constructing the lattice model. Trivial elementary excitations are labelled by a principal quantum number n = 1, . . . , N − 2 and have a mass m n relative to the kink mass m K given by [48,75] For every n, there are multiplets of these excitations labeled by r = n, n − 2, . . . , ≥ 0, the linear fundamental representations of so(2N ). These multiplets are fermionic (bosonic) for r (and thus also n) odd (even). The elementary fermion corresponds to n = 1 and thus transforms according to the defining vector representation of O(2N ). An exact result for the mass of this elementary fermion was derived in Ref. [49], namely in terms of a specific choice for the RG-invariant scale Λ MS , known as the modified minimal subtraction scheme when using dimensional regularisation. Note that for N = 2, the elementary fermion is not stable and decays into two kinks, i.e. m 1 = 2m K . In that case, Eq. (2.6) is providing a definition for (twice) the kink mass m K . When using a different regularisation scheme, such as the lattice Hamiltonian introduced next, the coupling and its UV dependence differ. As a result, the RG independent scales Λ defined from Eq. (2.4) need to be matched between different regularisation schemes.
In what follows, we obtain Λ MS = 8 e Λ lat using an exact solution of the lattice model in the limit N → ∞. A more standard yet involved Feynman diagram calculation of the scattering matrix in Appendix A proves that this relation is valid for all values of N .

Lattice Hamiltonian
To construct a lattice regulated version of the GN Hamiltonian, we follow the staggered fermion formulation from Ref. [62]. While this procedure is well known, we review it with some detail, in order to properly motivate our lattice proposal for the GN interaction.
One interpretation of the staggering procedure, which is useful for what follows, is to discretise the two components of the Dirac spinor at positions differing by half a lattice spacing, i.e. ψ c,1 (na) → φ c,2n / √ a and ψ c,2 ((n + 1 2 )a) → φ c,2n+1 / √ a, with a the lattice spacing. Furthermore, we choose the matrix α = γ 0 γ 1 appearing in the kinetic term of the Dirac Hamiltonian off-diagonal. With this choice the free massless Dirac Hamiltonian only couples derivatives of the first component to the second component of the Dirac spinor, and vice versa. For such terms a symmetric finite difference approximation of the derivative 1 1 The staggered fermion formulation is often introduced with the two spinor components discretized at the same positions and an asymmetric finite difference scheme, e.g. forward for ∂ψ1 and backward for ∂ψ2. Explicitly shifting the discretisation positions of the two components leads to the same end result, but is somewhat more aesthetically pleasing, and also clarifies how to deal with terms where both components (not their derivates) appear on the same position, e.g. ψ † 1 ψ2.
leads to e.g.
Combined with the requirement that α is real to make the O(2N ) symmetry explicit (and thus easier to preserve in the lattice model) leads to α = σ x , and the resulting lattice Hamiltonian is given by with K the (N -flavor) tight-binding or hopping operator where we have introduced Majorana modes χ m,n as φ c,n = (χ 2c−1,n + iχ 2c,n )/ √ 2 to make the O(2N ) symmetry of this lattice operator explicit.
Before adding the GN four point interaction, let us first discuss how to add an explicit mass term. Having fixed α = γ 0 γ 1 = σ x , the field theory allows for any choice β = γ 0 = cos(θ)σ y + sin(θ)σ z . The original proposal of Susskind in Ref. [62] was β = γ 0 = σ z , which is then trivially discretised into a lattice mass term on the doubled lattice as with ∆ = ma the mass in dimensionless lattice units. This clearly indicates how one-site translations on the staggered lattice flip the sign of the mass term, and can thus be related to a lattice remnant of the discrete chiral transformation ψ → γ 5 ψ. However, this lattice term breaks the O(2N ) symmetry, as can be made explicit by rewriting it in terms of the Majorana components. The alternative choice β = σ y yields terms involving both components of the Dirac spinor on the same position, which can be discretised on our staggered lattice by averaging one of the two components over the two nearby positions. This gives rise to an alternative lattice mass term, which takes the form of a staggered hopping ∆ 2 n (−1) n K n,n+1 (3.5) and thus respects the O(2N ) symmetry. This term is well known from the SSH model, where the alternating hopping strengths result from dimerisation. On the lattice, these two mass terms, resulting from two different choices for β (and thus, ultimately, a different choice of basis for the Dirac spinor in the continuum) are not equivalent. From the periodic table of topological insulators and superconductors [79,80], it is well known that the SSH mass term preserves sublattice symmetry (class AIII or BDI), which gives rise to a protected topological invariant labeled by Z. Sublattice symmetry is also known as chiral symmetry in that context, but we refrain from using this terminology, as it is clearly different from the chiral symmetry of the field theory relevant to our study, and which is broken by either mass term.
Writing the Hamiltonian terms in momentum space after blocking two sites, they take the form similar to the field theory Hamiltonian but with d(p) a periodic function of the lattice momentum p ∈ [−π, +π) on the blocked lattice. A gapped model has nonzero d(p) for all p. Sublattice symmetry imposes that d(k) is confined to a two-dimensional space, and the topological invariant corresponds to the winding number of d(p) around the origin. Both the kinetic term in Eq. (3.2) and the SSH mass term in Eq. (3.5) only have non-zero d x and . They lead to a well-defined winding number, which is non-zero for ∆ < 0 (shifting the unit cell definition is equivalent to ∆ → −∆), thus indicating a symmetry-protected topological (SPT) phase protected by either sublattice symmetry or bond-centred inversion, where the latter is also defined for interacting systems. Susskind's mass term [Eq. (3.4)] corresponds to d z (p) = ∆ and breaks the topological invariant. In the field theory, the kinetic term has a single component (i.e. d x (p) = p if α = σ x ), and so either choice of β is equivalent. While the winding number is undefined as momentum space is unbounded, topological features still manifest themselves when considering a domain between positive and negative mass, which gives rise to gapless edge modes, as described by Jackiw and Rebbi [78]. Hence, the SSH mass term provides a more faithful lattice description of the massive Dirac field.
We can rewrite the mass term from Eq. (3.5) as with the three-site operator which plays the role of a local order parameter, i.e. the lattice equivalent ofψψ. Whereas K n,n+1 has non-zero expectation value even with respect to the ∆ = 0 ground state, Σ n,n+1,n+2 is an absolute measure for the mass condensate. With this, it has now become straightforward to formulate a lattice Hamiltonian for the GN model, Phase diagram of the lattice GN Hamiltonian from Eq. (3.9) with additional SSH mass term from Eq. (3.5). The point g = ∆ = 0 is the lattice realisation of the N -flavour free fermion conformal field theory, which plays the role of a UV fixed point. For negative/positive mass perturbations the ground state is in a symmetry protected/trivial gapped phase respectively. The GN model has no explicit mass term and corresponds to a first order phase transition between those phases, i.e. it has two gapped ground states with zero and non-zero topological invariant respectively. The continuum limit of the Gross Neveu model approaches the c = N point via the red arrows, while the black arrows denote the continuum limit of N massive fermions.
where the interaction coefficient was changed from g 2 /2 to g 2 /4 as we associated one interaction term with every site of the doubled lattice. Doing so, this model has singlesite translation invariance, which corresponds to the lattice remnant of discrete chiral symmetry, as well as O(2N ) symmetry. The projective nature of the O(2N ) action on the single-site Hilbert space, which is discussed in Sec. 5, in combination with translation invariance enables the application of the Lieb-Schultz-Mattis theorem [81]: this model cannot have a unique, gapped ground state. It is critical for g = 0, but we expect the interaction to be marginally relevant and lead to a symmetry broken state for g = 0. The Mermin-Wagner theorem excludes the O(2N ) symmetry to be broken, thus leading to dimerisation, the lattice manifestation of a mass condensate, as the most likely scenario. By adding an explicit SSH mass term, a two-dimensional phase diagram is obtained, depicted in Fig. 1, where the Hamiltonian in Eq. (3.9) (i.e. the ∆ = 0 line) can be identified as a first order phase transition between the trivial and non-trivial SPT phase. While this is somewhat similar to the Z 2 -symmetry breaking phase of the Ising model being a first order line between the the explicit symmetry-broken regimes with positive and negative longitudinal field, the topologically distinct nature of the phases at both sides of the transition results in symmetry fractionalisation in the kink excitations that interpolate between the two ground states at the first order line. We confirm that the low-energy behaviour of this lattice model does indeed replicate all features of the Gross-Neveu field theory in the next sections using a large N calculation (Section 4) and by constructing an MPS ansatz (Sections 5 and 6). The MPS ansatz is not only used for numerical simulations for N = 2, but also provides further insight into the symmetry structure of the excitation spectrum for general N .
permutation symmetry corresponding to exchanging the different flavours in combination with monogamy of entanglement [82,83] can be used to argue that the ground state |Ψ will be a product state over the different flavours, where each flavour is described by the same state: |Ψ = |φ ⊗N . The energy of this state is given by with k n,n+1 and σ n,n+1,n+2 the single-flavour versions of K n,n+1 and Σ n,n+1,n+2 respectively. The GN interaction splits into N terms which act on a single flavour, and N (N − 1) terms which act across different flavours, and are transformed into a product of expectation values due to our product state ansatz: correlations between flavours vanish for N → ∞. In order to obtain finite results, we take the limit N → ∞ while keeping λ = g 2 (N − 1) fixed. As g 2 itself goes to zero, the self-interaction of the flavours vanishes. Minimising the energy with respect to φ yields, after adding a Langrange constraint for the normalisation, a self-consistent eigenvalue problem for the state |φ : where the Lagrange parameter E φ can be interpreted as the energy of a single flavor, and Assuming a dimerised solution, we set φ|σ n,n+1,n+2 |φ = (−1) n σ 0 . The state |φ of a single flavor is now determined as the ground state of the quadratic mean-field (i.e.

Hartree-Fock) Hamiltonian
in which λσ 0 plays the role of an SSH mass. The mean-field Hamiltonian is diagonalised by blocking the lattice and going to momentum space, which gives rise to single particle energies ε(p) = ± 4 sin 2 (p/2) + λ 2 σ 2 0 cos 2 (p/2). (4.5) with, as before, p ∈ [−π, +π) the lattice momentum on the blocked lattice. One can verify that the self-consistency condition for σ 0 is equivalent to minimising with e Ψ the energy density associated with the sites of the blocked lattice. The value of σ 0 is thus determined by the condition 1 λ = π −π dp 2π cos 2 (p/2) This can be further expanded as with K and E the complete elliptic integral of the first and second kind, respectively. An asymptotic expansion for small λσ 0 , which is the dimensionless mass and should go to zero to recover the continuum limit, yields As a result, the effective fermion mass is in the large-N limit given by such that m 1 is indeed proportional to the RG-invariant mass scale Λ that was introduced in Section 2, with µ = a −1 . In particular, by comparing to the N → ∞ limit of the exact result in Eq. (2.6), i.e. m 1 = Λ MS , we are lead to conclude that if we define then Λ lat = 8 e Λ MS and thus we should recover (1 + O(g 2 )) (4.12) in the continuum limit g → 0. However, the N → ∞ solution is in itself not sufficient to support this conclusion, as other N dependent scale factors might appear. A careful comparison between Λ lat and Λ MS using the fermion-fermion scattering amplitude at finite N leads to the same result, as explained in Appendix A. Henceforth, we omit the lattice spacing a, as it appears trivially in length or mass scale quantities and does not directly affect the distance to the continuum limit. So we stop differentiating between dimensionless lattice and field theory quantities.

Matrix product states
For finite N , correlations between the different flavours cannot be ignored, and the ground state of our lattice model is a fully correlated quantum state, both in the spatial and in the flavour direction. We now try to approximate this ground state using a MPS ansatz, which is known to capture the quantum correlations in low-energy states of gapped local Hamiltonians for quantum spin chains [84]. The MPS ansatz associates with every site n of such a spin chain a 3-leg tensor A(n) of size D n−1 × d × D n , with d the local Hilbert space dimension of the physical index. The left (respectively right) virtual index of size D n−1 (D n ) is then contracted with the right virtual index of the previous (left virtual index of the next) tensor, resulting in a correlated state whose bipartite entanglement for a cut between site n and site n + 1 is upper bounded by log(D n ), independent of the system size (in accordance with the area law for entanglement entropy in one-dimensional systems [85]). By defining a unit cell, i.e. a periodic n dependence in the tensors A(n), we can describe quantum states directly in the thermodynamic limit.
Our lattice model is easily translated into a spin chain using a Jordan-Wigner transformation where we introduce a linear ordering in the flavour direction c = 1, . . . , N , and thus associate N qubits or spins with each site, so that the local Hilbert space dimension is 2 N . We keep these N qubits together (as opposed to treating them as N individual sites) in order to preserve the O(2N ) symmetry and to be able to capture it in the MPS ansatz. Using this particular ordering in the Jordan-Wigner transformation, the generators of the associated Lie algebra so(2N ) transform into a sum of one-site operators so that the resulting symmetry transformations act on-site. The local Hilbert space of a site can be identified with the direct sum of the two fundamental so(2N ) spinor representations (each of which has dimension 2 N −1 ) with opposite total fermion parity.
In order to construct SO(2N ) symmetric MPS 2 , the virtual indices of the tensors should also carry representations of the group [86,87] and the local tensors A(n) should intertwine the representation on the right virtual index with the tensor product of the representations on left virtual and physical index. The representation on the physical index, i.e. the direct sum of the two fundamental spinor representations, is projective 3 . Therefore, if the right virtual index is associated with a linear representation (i.e. a direct sum of tensor representations of SO(2N )), then the left virtual index should also be projective (and thus be composed of spinor irreducible representations), and vice versa. We are thus naturally led to a two-site unit cell. This is the MPS manifestation of the Lieb-Schultz-Mattis theorem [88]: MPS represent finitely correlated (and thus gapped) states, and cannot be simultaneously invariant under translation symmetry and an on-site symmetry whose physical action is projective. As the on-site symmetry is continuous and cannot be broken, we thus propose an ansatz for the ground state with a two-site unit cell, in line with the expected dimerisation: By working with the representations of the Lie-algebra so(2N ), we effectively only impose SO(2N ) symmetry. However, we discuss the role of the additional mirror symmetry that extends SO(2N ) to O(2N ) for the particular case N = 2 which was used in our simulations, and find that it is unbroken. 3 The spinor representations of (S)O(N) are linear representations of the universal covering group, known as (S)Pin(N ).
We can also formulate MPS-based ansätze for elementary excitations on top of the ground state [89][90][91], as well as kink excitations that interpolate between the two ground states. Both topologically trivial excitations and kinks can be created by modifying a single tensor (which has an effect on an extended region) and building a proper momentum superposition (unlike in semiclassical studies where the kinks or solitons are localised in real-space). The ansatz for kinks, where the two different ground states (corresponding to a one-site shift of the unit cell) surround the new tensor, is diagrammatically represented as with once again p ∈ [−π, +π) the momentum on the blocked lattice. The new tensors (labelled B i ) carry an additional leg corresponding to the Hilbert space of the irreducible representation of the excitation that is being targeted. In the case of kinks, as depicted here, the two virtual legs of B tensors need to be identical and thus carry the same representations, and it follows automatically that physical symmetry sector of the kink states needs to be a spinor representation. These constructions are well-known in the case of half-integer spin chains [92], where they also correspond to the renowned result that the elementary excitations in e.g. the half-integer spin Heisenberg chain are spinors [93,94]. The analogous construction for topologically trivial excitations on top of a single ground state illustrates that these are labelled by linear irreducible representations of SO(2N ). For an in-depth review on these excited states and their implementation we refer to Ref. [91].
The variational ansatz for excitations gives rise to an energy-momentum dispersion relation, e.g. E K (p) for the kink state |K p , from which we can extract a range of mass scales related to its value, inverse curvature and higher derivatives at p = 0, where the dispersion relation has its minimum. Indeed, by fitting to the dispersion relation for small values of the lattice momentum p, we obtain two different mass scales 4 . Alternatively, expecting relativistic invariance, we can rewrite this expansion for the square of the energy as and thus interpret the ratio m K,1 /m K,2 as an effective speed of light, which should go to one if Lorentz invariance is obtained in the continuum limit.
In a theory near a relativistic continuum limit, information about the particle masses is also encoded in the ground state, more particularly in the spectrum of inverse correlation lengths. This is true for trivial excitations, by writing two-point correlation functions using the Källén-Lehmann representation [95]. To extract the kink mass from a correlation function, one needs to study the correlation of string operators. In the MPS language, the corresponding inverse correlation length can easily be extracted by studying the mixed transfer matrix, made from two different ground states. In this particular case, these two ground states are related by a one-site shift, and we define with ρ the spectral radius (largest magnitude eigenvalue). The right hand side of this equation gives an inverse length scale, and requires a (dimensionless) velocity to give an energy scale. Again we assume this velocity to be one as we approach the continuum limit, and we will directly compare m K,3 as another estimate of the kink mass. Note that the eigenvectors of the mixed transfer matrix, due to the different nature of the two legs on which they act, also transform according to spinor representations, reflecting the spinorial nature of the kinks to which the mixed transfer matrix is related. One could also calculate the leading eigenvalue in the trivial transfer matrix as an estimate for the fermionic correlation length an hence inverse fermion mass. While (inverse) correlation lengths converge slowly as a function of the MPS bond dimension, a scaling theory for their behaviour was recently developed based on a parameter δ that quantifies the level spacing in the logarithmic eigenvalue spectrum of the transfer matrix (which should become continuous in the infinite bond dimension limit). We refer to Refs. [17,96] and use these techniques to extrapolate reported mass values to the infinite bond dimension limit.

Simulation of the N = 2 mass gap
As a first application of our MPS simulations, we compute the different estimations of the mass gap m K of the kink that interpolates between the two vacua for N = 2. Hereto, we first find the optimal MPS representation of the ground state using the "variational algorithm for uniform matrix product states" (VUMPS) [97], which directly minimises the energy (density) (i.e. variationally) in the thermodynamic limit.
Our implementation of this algorithm, as well as the algorithm for computing the dispersion relation using the excitation ansatz, can be found in "MPSKit.jl" [98], an open source package for MPS algorithms using the scientific programming language Julia. This package builds upon "TensorKit.jl" [99], a lower level open source package for representing and manipulating tensors with arbitrary (abelian and non-abelian) symmetries.
Specifically for N = 2, we enforced the tensors to be representations of SO(4), or rather its universal cover Spin(4). This group is equivalent to SU(2)×SU(2) and irreducible representations are labeled by a tuple of two SU(2) quantum numbers, i.e. half integers or integers. The resulting representation is a projective (i.e. spinor) representation of SO(4) if only one of both quantum numbers is a half-integer. With both quantum numbers integer or half-integer, a linear (i.e. tensor) representation of SO(4) is obtained. This symmetry can easily be understood by considering the two sets of generators, corresponding to rotations in flavour space (odd fermion subspace of single occupancy) and in some pseudospin space (even fermion subspace of zero or double occupancy), exactly as in the Hubbard model at half filling [100]. The disconnected part of O(4) (or its double cover, Pin (4)) is generated by a fermionic particle-hole transformation on one of the fermion flavours, and has the effect of interchanging the two SU(2) factors. It thus results in degeneracies between sectors (j 1 , j 2 ) and (j 2 , j 1 ), i.e. whenever j 1 = j 2 , these two representations will always come together as (j 1 , j 2 ) ⊕ (j 2 , j 1 ), where the direct sum constitutes a proper (linear or projective) representation of O (4). While this extra symmetry is not enforced on our MPS representation, it does seem to be perfectly preserved in the ground states we find numerically. Exploiting the Spin(4) ∼ = SU(2) × SU(2) symmetry has enabled us to push the bond dimension to D ≈ 4000. Using the exact equation for the mass gap [Eq. (2.6)], the relation between Λ MS and Λ lat proven in Appendix A, and the fact that the kink mass is half the fermion mass for N = 2, we obtain that the dimensionless mass of the elementary kinks for N = 2 should approach the value m K = 8 e e π 1 √ 2π ge −π/g 2 (6.5) in the continuum limit. Fig. 2 depicts four extrapolated mass scales as a function of the coupling. The blue dots, referred to as m K,3 , i.e. the inverse correlation length extracted from the mixed transfer matrix, are relatively cheap to compute but require careful extrapolation towards δ = 0. This extrapolation is illustrated in Fig. 3. Here, we present a handful of linear extrapolations for each coupling. One fit considers only the 5 highest bond dimension simulations, other extrapolations discard the highest bond dimensions and use the 5 next best bond dimensions. The resulting variation in the extrapolated mass gap gives rise the error bars depicted in Fig. 2. In a similar manner, the largest correlation length extracted from the normal transfer matrix (topologically trivial sector) should be determined by the fermion mass, and half of this value should also provide an estimate of the kink mass in the continuum limit. It is depicted by the green dots in Fig. 2. Note that for large values g 4 (which are not relevant for the continuum limit), the fermion mass (as extracted from the inverse correlation length), is less than twice the kink mass. Hence, at those values of the coupling constant, the lattice model is likely to exhibit a stable particle in the topologically trivial sector.
The pink and red dots in Fig. 2 correspond to the gap and inverse curvature obtained from the excitation ansatz, referred to as m K,1 and m K,2 before. These points converge more quickly with bond dimensions and require little extrapolation, yet are more costly to obtain. We have only calculated these points for a selection of couplings. The observation that the value of the gap, its curvature and the inverse correlation length coincide clearly shows an emergent Lorentz symmetry with speed of light equal to one, as intended. To further illustrate this, we plot the kink dispersion relation at g = 1.4 in Fig. 4 and compare the dispersion to the relativistic prediction E 2 p = m 2 + p 2 . Even for relatively large lattice momenta up to p ≈ π/2 the correspondence is good. Note that the mass here is already of the order of 0.4 in lattice units, and we are thus already quite far from the proper continuum limit. For even larger values of g 2 , where the mass becomes of the order of one in lattice units, deviations between the different mass scales can be observed, as expected.
To compare our mass data to the proposed continuum limit of Eq. (6.5), it is useful to consider the right panel of Fig. 2 where log(m K ) is displayed as a function of 1/g 2 . For sufficiently small couplings the logarithm is dominated by −π/g 2 resulting in a linear relation that is ideal for fitting. A fit of the form log(m K ) = − π g 2 +log(g)+log (C) is shown in both panels and we find C ≈ 1.121 for the inverse correlation length data m K,3 . This value compares well with the expected result C = 8 e e π 1 √ 2π = 1.092, but shows a small overshoot of about 3%. Note that the logarithmic contribution to the fit complicates the fitting process and makes the fit parameter very sensitive to data in the small g regime, where the masses become extremely small and thus hard to pinpoint exactly. Nonetheless, we do conclude that this data provides ample evidence for the continuum limit of our lattice model being well described by the GN field theory, as intended.  Figure 3. Extrapolation of the inverse correlation length of the mixed transfer matrix (topological sector), corresponding to the mass estimator m K,3 , as a function of δ, the spacing in the logarithmic spectrum of the transfer matrix [96]. For every value of the coupling g 2 we show a handful of linear extrapolations towards δ = 0 each taking different points into consideration. For example at g 2 = 0.56 we have highlighted the 5 points with highest bond dimension and the corresponding extrapolation in green. Another extrapolation where we discarded the 4 points with highest bond dimension is highlighted in red. These extrapolate to slightly different masses which allows us to estimate the error on the extrapolation that is shown in Fig. 2.  7 Entanglement structure of the groundstate Another advantage of tensor network representations of quantum states is that they give full access to the entanglement structure of the state, which is an interesting concept in its own right as it provides an fresh perspective into the quantum correlations of the state and has thus received a lot of attention lately. In particular, the half-space reduced density matrix defines the entanglement Hamiltonian H E via which for a relativistic field theory is also known as the modular or Rindler Hamiltonian (corresponding to an accelerating observer). For a conformal field theory (CFT), the modular Hamiltonian can be mapped back to the original Hamiltonian using a conformal (logarithmic) transformation. In a gapped theory close to a CFT in the ultraviolet (UV), we expect the entanglement structure, which is anyway determined by UV modes of the theory 5 , to follow the CFT prediction closely up to length scales of the correlation length. The logarithmic mapping should thus transform the modular Hamiltonian onto the CFT Hamiltonian on a finite system with length approximately given by the logarithm of the correlation length in the system. This argument was recently formalised for CFTs perturbed by a relevant interaction by Cho, Ludwig and Ryu [101]. In what follows we we will first calculate the prediction for the entanglement spectrum for general N . We will then check that the prediction matches our simulations for N = 2 despite the marginal nature of our interaction term.

General result
Anticipating that the entanglement spectrum corresponds to the CFT spectrum on a finite system with open boundary conditions, we thus compute the spectrum for N free fermions (the UV fixed point of our model) on a finite interval of length L. As before, we denote with λ m (m = 1, . . . , 2N ) the 2N Majorana fields, and with λ m,1 and λ m,2 their two spinor components. After partial integration, the CFT Hamiltonian is given by Conformal boundary conditions can be of the type to which we refer as the Ramon type. In both cases, the prefactors where chosen such that the Majorana modes obey their usual anti-commutation relations To construct a Fock space we need to define normal fermionic modes. For k = 0 we can defineφ m (k) =λ m,1 (k) + iλ m,2 (k), which again transform under the fundamental (i.e. vector) representation of SO(2N ). The 2N zero zero modesλ m,2 (0) are grouped into N additional fermions α c with c ∈ 1, . . . , N . In terms of these operators and after proper normalisation so that Tr(e −2πH E ) = 1, the resulting entanglement Hamiltonian is given by for the Neveu-Schwarz boundary conditions and for Ramon boundary conditions. Due to the zero modes, all eigenvalues of H (R) will be at least 2 N fold degenerate. In particular, the ground state will be an SO(2N ) scalar in the Neveu-Schwarz case and a direct sum of the two fundamental spinor representations in the Ramon case. More generally, as higher excited states are obtained by acting with the vector operators φ † m on those ground states, all eigenspaces of H (NS) E will transform as tensor representations, whereas all eigenspaces of H (R) E will transform according to spinor representations. It is thus straightforward to relate these two towers of eigenvalues to the entanglement spectrum obtained across cuts corresponding to virtual bonds with linear representations (for NS) and with projective representations (for R).

Numerics for N = 2
For N = 2, the lowest excited state of the entanglement Hamiltonian with NS boundary conditions, i.e. k=1 in Eq. 7.8 is a (1/2, 1/2) quartet (the SO(4) vector representation) with gap : Where L = log(κξ kink ) is the typical length scale of the system after the comformal mapping. The free fit parameter κ takes care of setting the UV scale. In Fig. 5 we show the inverse gap   Fig. 3 ). The blue line is the CFT prediction 2 π L and describes the extrapolated data well for κ ≈ 2.83.
Filling the k = 1 mode twice results in an energy π L sextet (the antisymmetric rank-2 tensor representation of SO (4)). There are two possibilities to obtain energy 3 2 π L , namely by filling the mode k = 2 (a vector representation), or filling the k = 1 mode with 3 particles (an antisymmetric rank-3 tensor representation, which is equivalent to the vector representation for SO(4)). For the difference between the ground state energy in the two sectors we find which can be shown to converge to π 4L for sufficiently large L (or thus, exponentially large correlation lengths). Here, we have assumed that the same length parameter L can be used for the two types of boundary conditions (NS or R). These eigenvalues, and a few more, for both H (NS) and H (R) , relative to π/L, are depicted in the right panel of Fig. 6. In the left panel we show the corresponding ratios obtained trough extracted MPS data. The gaps − 1 2π log λ i /λ (N S) 0 are rescaled by twice the gap in the NS sector which according to Fig. 5 is approximatly π/L. As anticipated, towards the continuum limit (i.e. large L) all ratios converge to the predicted value. Note, however, that significant lattice effects for smaller L are present, as the logarithmic mapping in the definition of the modular Hamiltonian makes it exponentially harder for the entanglement spectrum (as compared to e.g. the excitation spectrum) to correspond with the continuum limit. In particular, additional degeneracies between different SO(N ) representations which are predicted by the CFT result are not exactly reproduced away from the continuum limit. Finally, the predicted value for the gap between the groundstate energies in the Raymond and Neveu Schwarz sector (cfr. Eqs. 7.11 and 7.10) is also shown; it matches well with the data even for smaller L where this is not expected.
To further highlight this last peculiar observation, we use the gap between the dominating entanglement eigenvalues in the trivial sector to predict the gap for the topological sector, by eliminating L between Eq. (7.11) and Eq. (7.10). Fig. 7 compares this predic- L=log( kink ) (1, 0) (0, 1) (1, 0) (0, 1) . We obtain excellent agreement, even for large couplings far away from the continuum limit, where the CFT prediction is no longer expected to hold. Indeed, beyond g 2 4, the correlation length ξ is less than a lattice site.
Finally, Fig. 8 shows the scaling of the total bipartite entanglement entropy along trivial and topological cuts of the MPS as a function of the logarithmic length scale  Fig. 3. The curves represent the entropies corresponding to the entanglement spectra predicted from entanglement hamiltonians Eqs. 7.8 and 7.9. We find good agreement except for very large L; this is due to ill converged MPS results, as was also the case in Fig. 2.
. L = log(κ fit ξ kink ). The curves are the predicted entropies calculated from reduced density matrices corresponding to Eqs. 7.8 and 7.9. Again there is good agreement apart from large L where the MPS results are not sufficiently converged.

Discussion and outlook
We have constructed a lattice regulated version of the GN model that preserves the full O(2N ) symmetry, and has a lattice remnant of the discrete chiral invariance. This prescription differs from the typical regularisation using Wilson fermions (i.e. the Gross-Neveu-Wilson model, which has recently received attention from the cold atoms community), and also takes a different prescription for the mass term as proposed in the original staggered formulation by Susskind. This prescription is well known in the context of the SSH model, but a Gross-Neveu type interaction resulting from it had, to the best of our knowledge, not been considered.
By studying this lattice model in the limit of large N -where mean field theory becomes exact-as well as at N = 2 using MPS simulations, we have established that its low energy behaviour replicates all the features (and in particular degeneracies) expected from the field theory. At the same time, we argued how the resulting lattice model lies at the first order phase transition between a trivial and topological insulator (according to symmetry class BDI in the ten-fold way), and much of the degeneracies in both the excitation and entanglement spectrum can be reinterpreted from that perspective. At the quantitative level, we observed that the non-perturbative behaviour of this marginally relevant interaction makes it especially challenging to accurately probe the continuum limit. The mass remains very small for a significant range of the coupling constant, and then shoots up quickly, so that the regime where MPS can probe the behaviour of the continuum limit is rather small.
As the spectrum of massive particles becomes more interesting for larger values of N , it would be interesting to also study the model in this regime. However, our results on the entanglement structure indicate why using MPS simulations for larger values of N is non-trivial. The entanglement structure, and in particular the entanglement entropy, is dominated by the UV CFT of the model, which is that of N massless free fermions. We thus anticipate a linear scaling of the entanglement entropy in the number of fermion flavours, which translates to an exponential scaling in N of the required MPS bond dimension, in order to obtain similarly accurate results. It is an interesting question whether exploiting the full (S)O(2N ) symmetry of the model could help to overcome this exponential scaling. However, this first requires that the necessary representation data (Clebsch-Gordan coefficients and/or 6j-symbols) of (S)O(2N ) are computed, as these are less readily available for general N .
Other potentially interesting directions of further research concern the phase diagram for finite values of the temperature and chemical potential, which are also within the scope of MPS simulations [13,22,30,36,38]. There is active interest in the possible existence of an inhomogeneous phase at sufficiently large values of the chemical potential [58,59]. This interest is again spurred by the similarity of the GN model with QCD. Due to the sign problem, probing the QCD phase diagram with lattice Monte Carlo at moderate densities and with realistic values of the quark masses (the regime interesting for heavy ion experiments) is near impossible. While an inhomogeneous phase in the GN phase diagram would result in breaking of translation invariance, a continuous symmetry of the field theory, there might be arguments to believe that the Coleman-Mermin-Wagner theorem does not apply and the GN model could indeed exhibit such a phase. Coleman explicitly assumed relativistic invariance in his version of the theorem, which is broken by the chemical potential, whereas more general arguments against continuous symmetry breaking rely on the specific dispersion relation and the counting of the would-be Goldstone bosons that restore the symmetry, which is non-trivial when breaking spacetime symmetries. It would be interesting to study if MPS techniques can shed a new perspective on this question, though infinite MPS simulations would also need to choose a particular unit cell and would also struggle with incommensurate filling fractions.
A final extension, which we explore in a future publication 6 , is to apply the discretisation scheme presented in this paper to the chiral extension (with full continuous chiral symmetry) of the Gross-Neveu model. Preliminary results indicate that the resulting lattice model has emerging continuous chiral symmetry along a critical line in the phase diagram that corresponds to a deconfined quantum phase transition. To relate Λ MS to our Λ lat , one needs to match the dimensional regularisation scheme to our lattice regularisation scheme. In particular, we require the first coefficient c 1 in the expansion The standard strategy to obtain such a matching is to compare results for a physical quantity, which by definition should be independent of the particular renormalisation scheme. We consider the two-fermion scattering S-matrix in the large energy/momentum regime, where perturbation theory is reliable, and the physics is well described by weakly interacting massless Dirac fermions. Notice that for the lattice regularised version 'large energy/momentum' E means Λ lat E 1/a, with the first inequality assuring the weakly interacting regime and the latter inequality assuring the QFT continuum regime. In Fig. 9 we display the different Feynman diagrams that contribute up to one loop to this scattering process. Notice that, as in the original Gross-Neveu paper [45], it is convenient to decompose the quartic term g 2 2 ψ a ψ a 2 in the QFT (2.1) (or − g 2 4 Σ 2 n,n+1,n+2 in the lattice Hamiltonian (3.9)) by introducing a Hubbard-Stratanovich field σ with trivial propagator −i and interactions: −gσψ a ψ a in the QFT and g √ 2 σΣ n,n+1,n+2 on the lattice (A.2) For reference, we first briefly discuss the computation in the MS scheme. To get the Feynman rules one first needs the free-field propagator (see e.g. [102]) using relativistic notation (with e.g. Figure 9. The Feynman diagrams contributing to the two fermion scattering S-matrix up to second order.
The Feynman rules then read: With the conventions of [102], writing S = 1 + iT , the tree-level diagram (a) then gives for iT : with the second factor K 1 on the first line arising from the projection on the particular fermion polarisations of the in-and out-modes that we consider, while the third factor K 2 (second line) arises from the colour conservation and momentum conservation. Notice that we do not consider the crossing diagrams for the outgoing legs, which gives terms in iT proportional to a different colour pre-factor δ a b δ b a . Using the standard machinery of dimensional regularisation and 'Diracology', we then obtain for the second diagram (b) in the MS scheme: with K = K 1 × K 2 the same polarisation factor and colour/momentum conservation factor as for (a). Furthermore for the diagrams (c) and (d) we find: Finally one can verify that the individual (logarithmic) UV divergencies in diagrams (e) and (f) cancel out when summed together. (e)+(f) is therefore scheme independent and plays no role in the matching. Notice that these separate UV divergencies in (e) and (f) would require a marginal counter term ∝ (ψ a γ µ ψ a )(ψ b γ µ ψ b ) which is prohibited by the full O(2N ) symmetry.
Collecting the different terms together we finally find up to order g 4 (neglecting scheme independent terms and terms ∝ δ a b δ b a ): for the S-matrix in the MS scheme.
Let us now turn to the computation of the same diagrams, but now with our lattice Hamiltonian (3.9). As a first ingredient we consider the free-field propagator. From the free Hamiltonian (3.2) one easily shows: 0| T (φ m (t) † φ n (u)) |0 = θ(t − u) with ω(k 1 ) = 2 sin |k 1 |, the particle energies (in lattice units) corresponding to the momenta k 1 (and we have omitted the colour indices). Using the Fourier representation of the stepfunction θ(t) we can rewrite the expression for the propagator above as: Notice that here we are not blocking the staggered sites, and the two Dirac spinor components now transpire in the two different branches −π/2 ≤ k 1 < π/2 and π/2 ≤ k 1 < 3π/2 of the (angular) spatial momentum k 1 of our single component fermions. In particular, with the identification of the physical (lattice-independent) momentum p 1 = 2k 1 a (A. 10) we have φ 2n+1 ≈ φ 2n for small momenta |p 1 | 1/a, corresponding to the Fourier transformed right-handed Dirac component ψ R in the Weyl representation, while the identification gives φ 2n+1 ≈ −φ 2n for |p 1 | 1/a, corresponding to the Fourier transformed left-handed Dirac component ψ L . One can easily verify that these identifications give the correct positive energy k 0 > 0 poles (in lattice units) in the first term of our propagator (A.9) k 0 ≈ a|p 1 |, both for the right-moving fermions (p 1 > 0) and left-moving fermions (p 1 < 0), while the second term of our propagator gives the proper poles for the anti-fermions.
We are now ready to write down the Feyman rules for our Hamiltonian (3.9) with decomposed interaction term (A.2) : Here for the expression of the fermion-propagator we summed the two terms in (A.9). Also notice the extra momentum structure in the vertex vis-á-vis the vertex rule (A.4) in the MS-scheme. For small physical momenta p 1 (see (A.10) and (A.11)) this structure simply expresses the fact that theψψ term couples the right-handed Dirac components (k 1 ≈ 0) to the left-handed components (k 1 ≈ π).
For the tree-level diagram (a) in Fig. 9 we then find (replacing p → k and q → h): (a) = ig 2 × cos(k 1 ) − cos(k 1 ) cos(h 1 ) − cos(h 1 ) Taking into account the different normalisations (e.g. k|k = 2πδ(k − k ) here, and p|p = 4πE p δ(p−p ) in the QFT computation) one can show that this reduces to the treelevel MS result (as it should) for small momenta |p 1 a| 1, either on the right branch from Eq. (A.10) or on the left branch of Eq. (A.11). Notice that for a non-vanishing scattering in this continuum limit, we either need k 1 , h 1 ≈ 0 (right-movers) and k 1 , h 1 ≈ π (leftmovers) or the other way around. Since for the remainder we are only interested in the continuum QFT limit, we can effectively set K 1 = 4 and anticipate that k 1 − k 1 = π + ∆ 1 , with ∆ 1 1.
It is easy to show that we get the same term from the diagram (d), while as we explained above we can forget about the (e) and (f) diagrams for the matching as these diagrams will be finite and universal for a manifest O(2N ) symmetric regularisation scheme like ours. Collecting all the relevant diagrams we then find for the QFT S-matrix in our lattice regularisation: and if we plug this result in the definition (2.4) for Λ, we obtain: by which we have shown explicitly that the matching result based on the large N mean-field computation of section 4, generalises to any finite N , in particular to N = 2.