A higher-order topological twist on cold-atom SO(5) Dirac fields

Ultracold Fermi gases of spin-3/2 atoms provide a clean platform to realise SO(5) models of 4-Fermi interactions in the laboratory. By confining the atoms in a two-dimensional Raman lattice, we show how this system can be used as a flexible quantum simulator of Dirac quantum field theories (QFTs) that combine Gross-Neveu and Thirring interactions with a higher-order topological twist. We show that the lattice model corresponds to a regularization of this QFT with an anisotropic twisted Wilson mass. This allows us to access higher-order topological states protected by a hidden SO(5) symmetry, a remnant of the original rotational symmetry of the 4-Fermi interactions that is not explicitly broken by the lattice discretization. Using large-$N$ methods, we show that the 4-Fermi interactions lead to a rich phase diagram with various competing fermion condensates. Our work opens a route for the implementation of correlated higher-order topological states with tunable interactions that has interesting connections to non-trivial relativistic QFTs of Dirac fermions in $D = 2 + 1$ dimensions.

Symmetry plays a primary role in our most fundamental theories of Nature.So far, all forms of matter observed in the laboratory can be ultimately described by the standard model [1], a relativistic quantum field theory (QFT) that contains Dirac fermions locally coupled to both scalar and gauge bosons, and is invariant under the Lorentz transformations.The role of symmetry goes beyond this relativistic invariance, as the specific form of the local fermion-boson couplings is dictated by the invariance of the QFT under various groups of local, so-called gauge, symmetries [2].Additionally, there are also global symmetries that leave the standard model invariant.These symmetries can be spontaneously broken at a certain energy scale, such that the vacuum of this QFT displays a certain order parameter that is no longer invariant under the action of the symmetry.An example of great relevance in high-energy physics is chiral symmetry breaking [3], which leads to a so-called chiral condensate, and accounts for most of the mass of the matter surrounding us.In the early days of quantum chromodynamics (QCD) [1], various effective QFTs that account for mass generation by chiral symmetry breaking were explored, including 4-Fermi QFTs such as the Nambu-Jona-Lasinio [3][4][5][6] and Gross-Neveu [7][8][9][10][11] models.These QFTs describe self-interacting Dirac fermions with various forms of quartic interactions and, moreover, can be defined in different spacetime dimensions leading to clear analogies with higher-dimensional non-Abelian gauge theories [12,13].
Besides serving as effective models that capture some of the phenomenology observed at particle colliders qualitatively, low-dimensional Dirac QFTs also emerge in condensed matter and atomic molecular and optical (AMO) physics, providing a playground to test the QFT predictions quantitatively.Systems such as graphene [14,15], topological insulators [16,17], and ultracold atoms in optical lattices [18,19], are clear examples leading to low-dimensional Dirac matter [20,21].We remark that, away from D = 1+1 dimensions, it is not straightforward to find experimental setups which, in spite of being highly non relativistic, are accurately described by Lorentz-invariant Dirac QFTs at low energies.Moreover, in the presence of interactions, these QFTs can present a number of non-trivial properties for D < 4 [22][23][24] that are still the subject of active research.A particularly-relevant case is that of Dirac QFTs in D = 2 + 1 dimensions with a Gross-Neveu 4-Fermi term [25,26] or variants thereof [27], which can give rise to strong correlations and novel critical phenomena.This type of 4-Fermi models leads to different phase transitions that are no longer characterised by the above chiral condensate but, instead, require finding other symmetry-breaking channels with their associated order parameters.These generalised 4-Fermi models have been the subject of renewed interest in recent years (see the reviews [28,29] and references therein).
In this work, we are interested in models of synthetic Dirac QFTs that can be implemented with ultracold atoms [30].A rather unique possibility of these systems is that, in addition to the emerging Lorentz invariance, one can actually tailor the other symmetries, either local or global.These platforms can thus be used as analogue quantum simulators [31][32][33][34][35] for a specific QFT of interest, allowing us to test theoretical predictions in a controllable experimental environment, as well as to explore other regimes that cannot be accessed using current analytical or numerical methods, such as real-time dynamics or finite fermion densities [36,37].In particular, we focus on spin-3/2 neutral atoms at ultra-cold temperatures, as their s-wave scattering leads to 4-Fermi terms that naturally yield a large symmetry group of SO (5) transformations [38,39].When these cold atoms are loaded on standard optical lattices, one obtains a non-relativistic SO(5) Hubbardtype [40] model that is interesting in its own right [38,41,42], and leads to a rich phase diagram already in D = 1 + 1 dimensions [43][44][45][46].In fact, the role of SO (5) symmetry in condensed-matter systems is broader, as it also plays a key role in certain theories with competing magnetic and superconducting orders [47,48].
We show in our work that, by including additional Raman beams that interfere with the standing wave underlying the above optical lattice, we can enter a new regime in which the SO(5) 4-Fermi model corresponds to a specific discretization of a relativistic QFT of self-interacting Dirac fermions.More specifically, we demonstrate how this lattice regularization, which in general breaks explicitly the rotational symmetry, can actually preserve a hidden SO(5) π/2 rotation, and provide a neat route to explore correlation effects in higher-order topological insulators (HOTIs) [49][50][51].In particular, this lattice discretization corresponds to an anisotropic version of twisted-mass Wilson fermions which, as we show, leads to flat bands and strictly-localised zero-energy corner modes protected by the hidden SO (5) symmetry.These anomalous corner modes are a boundary manifestation of a non-trivial topological invariant in the bulk [52], which connects to the phenomenology explored for other lattice models of higher-order topological states.In these models, the study of the effects brought up by including interactions has seen an increase of activity very recently [53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70], but still remains largely unexplored in comparison to the situation in their first-order counterparts (see the reviews [71][72][73][74]).
We contribute to this line of research by exploring the ef-fect that the SO(5) 4-Fermi interactions have on the aforementioned Wilson-fermion HOTI.In connection to the fermion condensate for chiral symmetry breaking in high-energy physics, we show that our model accounts for a competition between various possible condensation channels and that, at sufficiently strong interactions, the HOTI phase gives way to a pseudo-scalar fermion condensate where the hidden SO(5) rotational symmetry that protects the HOTI gets spontaneously broken.We present a non-perturbative account of this phenomenon based on the large-N limit of this 4-Fermi QFT, where one considers N flavours of the Dirac fermions coupled by the quartic SO(5) interactions.By resuming the leadingorder Feynman diagrams for N → ∞, we calculate the effective potential, and perform a minimization that allows us to infer the values of various condensates.Moreover, this large-N techniques can be readily used to obtain an estimate for the many-body topological invariant, allowing us to chart the entire phase diagram of the model showing, in addition to the aforementioned condensates, correlated HOTIs and trivial band insulators.Since the model studied can be realised in experiments of spin-3/2 neutral atoms in Raman optical lattices, possible future experiments could test these predictions and their connection to non-trivial strongly-coupled QFTs.

II. EUCLIDEAN 4-FERMI FIELD THEORIES WITH A TWIST
In this section, we introduce our model of interacting Dirac fields in D = 2 + 1 dimensions, which is motivated by a specific Kaluza-Klein-like dimensional reduction.We also discuss a non-standard lattice regularization that will allow us to study the non-perturbative phenomena induced by fermionfermion interactions on higher-order topological groundstates.
A. Dimensional reduction and SO(5) → SO (3) Our model of self-interacting Dirac matter is built from a relativistic QFT of fermions with rotationally-invariant 4-Fermi interactions.As discussed in Appendix A, the partition function of this QFT can be written as a path integral [75] over two independent Grasmmann spinors ψ(x), ψ(x) [76], which represent the Dirac fermions in a 3-dimensional Euclidean spacetime with imaginary time x = (τ, x x x), where x x x = (x 1 , x 2 ).The path integral is expressed in terms of an Euclidean action that contains two terms S = S 0 + S int .The first one describes free Dirac fermions with two possible mass terms where µ ∈ {0, 1, 2} labels the spacetime coordinates, ∂ µ = ∂ /∂ x µ , and m 1 , m 2 are the corresponding bare masses, which will be latter connected to a mass twisting.Here, we have used Einstein's convention of repeated index summation, and natural units h = c = 1.The set of 5 gamma matrices {γ 0 , γ 1 , γ 2 , γ 3 , γ 5 } fulfill {γ a , γ b } = γ a γ b + γ b γ a = 2δ a,b 1 d s , which can only be satisfied by considering Grassmann fields with d s = 4 spinor components.We recall that, in an Euclidean metric, these gamma matrices are all Hermitian, and can be defined by various possible choices of tensor products of operators within the Pauli basis {1, σ x , σ y , σ z } [77,78].
Although the specific choice of gamma matrices is arbitrary at the level of the QFT (1), the implementation based on spin-3/2 cold-atom gases discussed in Sec.III B fixes our choice of the spacetime gamma matrices to In addition, the remaining gamma matrices are also fixed as As discussed in Appendix A, this set of matrices actually forms a reducible representation of the Clifford algebra in the underlying 3-dimensional spacetime.These gamma matrices can be used to define the generators of Lorentz transformations, which correspond to a spinor representation of the SO(3) rotations in the Euclidean metric.We remark that such Lorentz transformations could also be generated using an irreducible representation of the Clifford algebra, which would only require using two-component spinors [77,78].However, this would not allow for the two independent mass terms in Eq. ( 1), as the spacetime gamma matrices would already exhaust all the possible mutually anti-commuting Hermitian matrices, e.g.γ 0 = σ z , γ 1 = σ x , γ 2 = σ y .In this case, there can only be a single mass term m 0 ψψ that breaks parity [26].As we will see, having two anti-commuting masses (1) plays a key role in the phenomena discussed in our work.
As discussed in detail in Appendix B, a different perspective motivating the choice of the action (1) is that the reducible representation (2)-( 3) is the result of a Kaluza-Klein-type compactification.In the original Kaluza-Klein context [79], gravity and electrodynamics in D = 3 + 1 dimensions were shown to result from the compactification of an extra dimension in a 5-dimensional theory of pure gravity.In the current context (1), the situation is much simpler, as one only needs to consider a QFT of Dirac fermions in a higher 5-dimensional spacetime coupled to background fields that will induce the two mass terms [78,80].In the larger spacetime, SO(5) symmetry is manifest in the action, and can also be used to understand the structure of the Lorentz-invariant 4-Fermi selfinteractions, as we now discuss.Considering the reducible representation of the Clifford algebra, we can define SO(5)invariant 4-Fermi interactions that can be written as follows where we have introduced the Euclidean fermion current as J µ = iψγ µ ψ.This action can again be interpreted from the perspective of dimensional reduction, where the first five terms can be written as the squared norm of a higherdimensional vector whose components are fermion bilinears.Therefore, the norm of the vector is conserved under the SO(5) Lorentz transformations (see Appendix B).In addition, the last term in Eq. ( 4) is a scalar under these SO(5) spatial rotations, and thus also remains invariant.From the perspective of the higher-dimensional parent theory, these quartic terms correspond to a linear combination of the standard Gross-Neveu and Thirring [81,82] interactions.This is very different from the interactions allowed by an irreducible representation, which can all be reduced to a single Gross-Neveu quartic interactions.On the other hand, for a reducible representation, the physics of Gross-Neveu and Thirring QFTs in D = 2 + 1 turns out to be very different [83].Therefore, working with reducible gamma matrices leads to a richer dimensionallyreduced QFT with more interaction channels than those allowed by an irreducible representation.The structure of Eq. ( 4) will yield a competition of various fermion condensates with the aforementioned topological groundstates.

B. Anisotropic mass twisting of Wilson fermions and the explicit breakdown of SO(3) symmetry
After setting up the continuum QFT ( 1)-( 4), we now discuss how higher-order topological [49][50][51] and trivial groundstates may arise as the result of a non-standard lattice regularization.We introduce this regularization by the use of Wilson fermions [84] for the more standard first-order topological insulators [85,86], and then highlight the "twist" that is required to connect to higher-order topology.
Let us start by noting that all of these regularizations are related to the problem of fermion doubling [87,88] in lattice field theories [89].In particular, we shall consider a nonzero lattice spacing a along the spatial directions, such that contains an even number of lattice sites per axis, leading to N s = N 1 N 2 as the total number of lattice sites (see Fig. 1(a)).The spatial derivatives in Eq. ( 1) are substituted by finite differences, while the Euclidean time remains continuous τ ∈ R.This asymmetric treatment of the spacetime will allow us to make a direct connection with the Hamiltonian approach to field theories in section III.The fermion doubling can be understood by writing the free action (1) after this regularization which, in momentum space, reads where k = (k 0 , k k k) is the three-momentum, an we have introduced the short-hand notation k := 1 2π .In the expression above, we have defined k0 = k 0 ∈ R using the zerotemperature limit of the Matsubara frequencies [75].Additionally, the spatial components of the momentum k k k are related to the corresponding crystal momenta We recall that, as a consequence of the lattice regularization, the crystal momenta are quantised within the first Brillouin zone k k k ∈ BZ, namely k j = −π/a + 2πn j /aN j with n j ∈ Z N j .

By a Taylor expansion
around any of the four Dirac points one finds that the long-wavelength action stemming from Eq. ( 5) coincides with that of a massless Dirac fermion (1).Altogether, the long-wavelength infra-red (IR) features are governed by N D = 4 Dirac fermions instead of one.Typically, one refers to ℓ ℓ ℓ ∈ {(0, 1), (1, 0), (1, 1)} as the spurious fermion doublers, each of which has a different emergent chirality γ 5 ℓ ℓ ℓ = (−1) ℓ 1 +ℓ 2 γ 5 .The presence of these spurious doublers can change the physics considerably, specially when adding further interactions such as those in Eq. ( 4).The idea of K. Wilson to cope with fermion doubling [84] was to add a momentum-dependent Wilson mass, which sends these doublers to the ultra-violet (UV) cutoff scale of the QFT, such that they become very massive and are not expected to interfere with the low-energy physics of the remaining Dirac fermion ℓ ℓ ℓ = (0, 0).In Appendix C, we show that the standard Wilson regularization amounts to setting m 2 → m 2 = 0 and m 1 → m 1 (k k k) in Eq. (C1), which leads to two copies of a Chern insulator [90][91][92][93][94].In this sense, the use of a reducible representation of the gamma matrices is rather trivial, and we could have obtained the same topological features using an irreducible representation and two-component spinors.
Note, however, that we have only used one of the mass terms in Eq. ( 5), and followed the Wilson regularisation scheme verbatim.In this work, we explore a different nonstandard regularization which, although having a similar effect on the doublers, leads to very different manifestations of topology.We consider an anisotropic twisted Wilson mass regularization S = S TM 0 + S int , in which the free part is expressed in momentum space as This contains the above anisotropic twisted Wilson mass where r ∈ R is the analogue of the Wilson parameter in the standard regularization of App. C. The interacting part of the action S int is defined by considering the 4-Fermi terms in Eq. ( 4), but discretizing the spatial coordinates such that the integral becomes a shorthand for d 3 x := a 2 ∑ n n n∈Λ s dτ.Twisted-mass Wilson fermions have been previously considered in lattice gauge theories [95][96][97][98][99][100][101][102][103] (see [104] for a detailed account), although the twisting procedure is very different from the one considered in our work.In order to understand the differences, let us describe mass twisting in a broader context.We start by focusing on the standard situation, which involves even-dimensional spacetimes such as As discussed in Appendix A, the usual mass twisting follows from an axial rotation of angle θ on the standard Dirac mass term mψψ, leading to a pair of anticommuting mass matrices proportional to m 1 = m cos θ , m 2 = m sin θ (A7).The new mass term for θ ̸ = 0 breaks parity (B6) explicitly, unless the rotation/twisting angle is promoted to a pseudo-scalar axion field θ → θ (τ, x x x) [105][106][107].Let us now consider this mass twisting for a Wilson-fermion lattice regularization, which upgrades the parity-invariant mass to a momentum-dependent one m 1 → m 1 (k k k) similar to Eq. (C1), but accounting for the larger spatial dimensions.This Wilson term leads again to a different mass for each of the Dirac points, and can lea to a non-zero topological invariant.From the perspective of first-order topological insulators [92], a non-zero twisted mass m 2 > 0 breaks both parity and timereversal, making the topological invariant no longer quantised [108] which to connects to axion electrodynamics [109].
All of these effects require having an odd number of Dirac points with negative Wilson masses which, in the context of gauge theories such as QCD, would shift the value of the vacuum theta angle.When the goal is to recover the continuum limit of QCD, it is sensible to avoid this possibility unless one is using a discretization based on domain-wall fermions with an extra spatial dimension [110][111][112].Accordingly, the works on lattice QCD with Wilson fermions typically work in regions of parameter space that are close to a critical line in order to recover a continuum limit, but making sure that these topological contributions to the theta angle are zero.The motivation to include a mass twisting is then completely different from the above discussion.When trying to improve the lattice regularization to achieve a faster convergence to the continuum [113][114][115], the approaches based on Wilson fermions [116,117] can take advantage of the mass twisting [96,97].In fact, as first shown in [98], by working with a maximal twist angle θ = π/2, one can find an automatic O(a) improvement that only requires controlling a single parameter of the model.
If we now move back to our D = 2 + 1 dimensions, as noted above, an irreducible representation of the gamma matrices would forbid considering a mass twisting.On the other hand, our reducible representation (2)-(3) permits additional mass terms (1), yielding a a 3-dimensional version of the mass twisting.Still, we remark that this twisted Wilson mass would be very different from our anisotropic mass twisting in Eq. (9).As discussed in [104], the usual mass twisting for Wilson fermions can be transformed to a physical basis where, rather than the bare masses, one rotates the Wilson term responsible for the momentum dependence of the mass.Exploring angles different from the full twist, e.g.θ = π/4, would bring us closer to the type of Wilson mass twisting of Eq. ( 9).Yet, there is a fundamental difference, the Wilson mass twisting considered in lattice gauge theories is always isotropic.We are instead considering that the dependence on momentum of the twisted Wilson masses (9) is highly anisotropic: m 1 (k k k) only depends on k 1 , and m 2 (k k k) on k 2 .In turn, as will become clear below, this means that the real-space Wilson term is anisotropic; it is different when the fermions tunnel to neighbouring sites along each of the two spatial directions.The other important difference is that, instead of considering a mass twisting that combines ψψ and ψiγ 5 ψ, we are admixing ψiγ 3 ψ and ψiγ 5 ψ terms, both of which are parity invariant in D = 2 + 1 dimensions.This will be very important to get a model with a hidden SO(5) rotational symmetry, which will protect the higher-order topological state.
In the following section, we will show in detail how the groundstate of the free lattice action can be characterised by a topological invariant which is, however, distinct from that (C6) of the Chern insulators discussed above.Indeed, as shown in Sec.III, the groundstate in this case corresponds to a higher-order topological insulator (HOTI).The bulkboundary correspondence [118] leads to a boundary manifestation that differs from the edge states of Chern insulator, as we find zero-energy states that are only localised in the corners of the spatial lattice.To make this connection clearer, we start by introducing a Hamiltonian version of this QFT.

III. COLD-ATOM HAMILTONIAN FIELD THEORY
In this section, we present the Hamiltonian of the above Euclidean field theory (8), which will be useful when discussing the HOTI, and a possible cold-atom implementation.
A. The Creutz-Hubbard multi-layer model Since our discretization keeps the imaginary time continuous, one can also describe the system through a Hamiltonian lattice field theory by rotating back to real time tau → −it.In the Hamiltonian formulation [1,119], one works with field operators instead of Grassmann variables.We thus define ) in terms of fermionic creation-annihilation operators defined on the lattice sites (see Fig. 1 (a)), which are supplemented with the following equal-time anti-commutation relations In a Minkowski spacetime, the adjoint is no longer independent but, instead, related to the creation operators Ψ n n n = Ψ † n n n γ 0 .The Hamiltonian operator H governing the dynamics of these fields can be found from the partition function discussed in Appendix A, recalling that the basis of fermionic coherent states is used to write the partition function Z = Tr{e −β H } as a path integral over the Grassmann fields ψ n n n , ψ n n n [75].The identified operator can be written as the sum of two terms H = H TM 0 + H int .The free term H TM 0 is obtained by the discretization of the spatial derivatives of Eq. ( 1) in terms of finite differences, which leads to tunnelling terms between nearest neighbours.Additionally, the anisotropic twisted Wilson mass is also realized by including tunnelling terms that give momentumdependence to the local masses in (1) according to Eq. (9).Altogether, this leads to the quadratic lattice Hamiltonian where we have introduced the tunnelling matrices These depend on the following tunnelling strengths as well as on the Dirac α-matrices defined as Note that the above matrices are still expressed in terms of products of the Euclidean gamma matrices in Eqs. ( 2)-(3).In Minkowski spacetime, it is customary to work with { γa , γb one recovers the standard conventions for the Hamiltonian lattice field theory of Dirac fermions [1].Using the standard definitions of the Dirac α and β matrices, we also find where the later coincide with the expressions in Eq. ( 14).
In the above lattice Hamiltonian (11), we have also introduced the mass matrices which are expressed in terms of the following parameters In addition to this quadratic Hamiltonian (11), the 4-Fermi terms in Eq. ( 4) lead directly to the quartic interactions (18) In the next section, we will see how this specific interaction emerges naturally when considering spin-3/2 Fermi gases tightly-confined by optical potentials [38,39].We will argue that this connection fixes the choice of the α and β matrices to those in Eqs. ( 14)- (15), and thus forces the choice of gamma matrices (2)-(3) in our 4-Fermi QFT in Eqs. ( 1) and (4).
In the Hamiltonian formulation, the discussion of the SO(5) invariance must be revisited in light of the definition of the adjoint operator below Eq. (19).In fact, the Euclidean SO(5) Lorentz invariance must now be described in terms of SO (1,4) Lorentz transformations, where the boosts do not admit a unitary spinor representation [1].One could define a completely analogous Kaluza-Klein compactification, where the above Hamiltonians arise from a 5-dimensional parent model regularised on a lattice.In analogy to the Euclidean action, the continuum limit of the dimensionally-reduced model is expected to recover the lower-dimensional SO(1,2) invariance.On the other hand, we are not only interested in the continuum limit, but also in the HOTI phases of the theory where one can go beyond this continuum emergent symmetry.From this perspective, we should look for other transformations, including discrete spatial transformations, which correspond to exact symmetries of the full lattice model in (11) and (18).
As depicted in Fig. 1 (a), the non-interacting Hamiltonian (11) can be interpreted as a multi-layer fermionic model with both intra-and inter-layer tunnellings.An aspect that will be important in our analysis below is that there are certain background π-fluxes that dress the tunnelling along certain plaquettes involving the inter-layer synthetic dimensions (see Fig. 1 (b)-(c)).These fluxes lead to flat-band regimes, and a generalization of the so-called Aharonov-Bohm cages [120], which become very useful to understand the bulk-boundary correspondence.In this regard, our model can be considered as a higher-dimensional multi-layer version of the Creutz ladder [121][122][123][124][125][126][127].Moreover, the quartic interactions (18) are purely local, and can thus be interpreted as a Hubbard-like interaction.For D = 1 + 1 dimensions, the Hubbard interaction maps exactly onto a Gross-Neveu quartic term [128][129][130][131], although one could also use a bosonic species to mimic an auxiliary field that carries the Gross-Neveu interactions [132].For D = 2 + 1 dimensions, when working with an irreducible representation and two-component spinors, the Hubbard interaction maps again into a Gross-Neveu quartic term [93,94].In the current reducible case, where we have four-component spinors, the Hubbard-type interaction is richer and contains both inter-layer density-density interactions among all pairs of spinor states, as well as spin-changing collisions that involve effective inter-layer pair tunnellings (see Fig. 1 (d)).

B. Spin-3/2 atoms and 4-Fermi interactions
In this section, we present the details of how spin-3/2 fermionic atoms can naturally lead to the 4-Fermi interactions in Eq. ( 18).This is an example of the unique opportunity emphasised in the introduction: the possibility of tailoring local and global symmetries that connect to interesting models of high-energy physics.This brings us closer to the field of quantum simulators [32][33][34][35]; controllable quantum manybody systems that behave according to a specific model of interest [31].In the context of quantum simulators for highenergy physics (see the reviews [133][134][135][136][137][138][139][140][141]), there have been several proof-of-principle experiments showing the quantum simulation of relativistic QFTs [18,19,, including lattice gauge theories [119,179].The case of gauge theories is particularly demanding in terms of the required resources, as the tailored symmetries must be local and, ultimately, non-Abelian, requiring the introduction of additional gauge degrees to allow for this local symmetries.On the other hand, for synthetic Dirac matter with quartic interactions, the requirements are in principle milder, as one restricts to global and spacetime symmetries, including non-Abelian ones, but one can dispense with the extra gauge degrees of freedom.
Let us consider a gas of fermionic neutral atoms that can be tightly confined by optical potentials in a square lattice We emphasise that the physical lattice spacing is set by half the wave-length λ L /2 of the laser that leads to the optical-lattice potential, which is kept fixed in the experiment.This physical lattice spacing will not be mapped onto the lattice-field-theory spacing a, which must be sent to a → 0 to recover the continuum limit.Another difference is that, in second quantization [180][181][182], the atoms are described by dimensionless operators f † n n n,σ ( f n n n,σ ) that create (annihilate) an atom in the position specified by n n n, and in the internal electronic internal state given by σ ∈ S σ .The set S σ generally depends on the particular type of atom, and its specific isotope, which can also control the bosonic/fermionic nature of the operators.We will be interested in the fermionic case, where In the tight-binding regime, the system is thus described by a spin-conserving Hamiltonian, and its non-interacting part can be written [180][181][182] in second quantization as where t j is the standard nearest-neighbour tunnelling coupling along the j ∈ {1, 2} axis, and we set t1 = t2 =: t.In order to find the set of internal states S σ , we need to consider the atomic energy level structure, focusing in particular in the groundstate manifold.For instance, in the case of 6 Li Alkali atoms, we have principal quantum number n = 2, total orbital angular momentum L = 0 and spin S = 1/2 which, in spectroscopic notation leads to the 2 2 S 1/2 groundstate manifold.The total nuclear spin is I = 1, which leads to a couple of hyperfine levels with total angular momentum F ∈ {1/2, 3/2}.If we focus on the lower-energy state F = 1/2, the set of internal states is given by the two possible Zeeman sublevels M F = ± 1 2 , namely S σ = {− 1 2 , + 1 2 } =: {0, 1}.At sufficiently cold temperatures, the scattering of the dilute Fermi gas is dominated by s-wave collisions between pairs of 6 Li atoms, which mostly contribute [180] to a Hubbard densitydensity interactions that can be written as where n n n n,σ = f † n n n,σ f n n n,σ is the number operator, and is the Hubbard coupling strength [30].Here, k L = 2π/λ L is the laser wavevector, and E R = k 2 L /2m a is the recoil energy of the 6 Li atoms of mass m a .In this expression, we have introduced the opticalpotential depths along the different axes, which will be constrained to V 0,z ≫ V 0,x ,V 0,y such that the dynamics takes place within the xy plane.Finally, a key quantity in the Hubbard coupling is the s-wave scattering length a 0 , which only depends on singlet scattering channel [30].At cold temperatures and within the lowest hyperfine multiplet, the inter-atomic potential is rotationally invariant within the total angular momentum of the colliding pair F F F t = F F F 1 + F F F 2 , which in this case leads to a couple of channels with F t ∈ {0, 1}.Due to Pauli exclusion principle and the effective contact interactions for the s-wave channel, only the singlet case F t = 0 is allowed, which is described by the above s-wave scattering length a 0 .We note that the Hubbard interaction in Eq. ( 21) has a global SU(2) symmetry and, moreover, its strength can be controlled via Feshbach resonances by e.g.applying a magnetic field [183].
This type of interaction (21) would suffice to make connections to Gross-Neveu interactions for an irreducible representation of the gamma matrices [93,94,129], since the spinor components are only two.In this work, however, we are interested in reducible representations with a larger number spinor components S σ = {0, 1, 2, 3}, where a larger non-Abelian symmetry appears in the interactions.A well-known example of large non-Abelian global symmetries in the scattering appears for other atomic species, such as 87 Sr Alkalineearth atoms.In this case, there are two valence electrons, and the groundstate manifold has principal number n = 5, and vanishing total spin and orbital angular momentum S = L = 0, leading to the manifold 5 1 S 0 .For vanishing J = 0, there is no hyperfine splitting due to the nuclear spin, such that F = I = 9/2, and we get a single multiplet with N = 10 Zeeman sub-levels M F ∈ {−9/2, −7/2, • • • , 9/2}.Since there is no hyperfine coupling, the atoms all interact with an s-wave scattering length that is independent of the nuclear features and, thus, equal for all of the N = 10 sub-levels U σ 1 σ 2 = U 0 ∀σ 1 ̸ = σ 2 .Hence, Eq. ( 21) has an exact SU(N) symmetry [184,185].When considering also the long-lived 5 3 P 0 level, one gets more flexibility, leading to the so-called two-orbital SU(N) Hubbard models [186].When considering a mixedspecies Fermi gas with a couple of alkaline-earth atoms, the inter-and intra-orbital scattering preserve the SU(N) symmetry.Provided that one can control their corresponding scattering lengths via Feshbach resonances, there are specific conditions where the two-orbital SU(N) interactions would connect to the Gross-Neveu term between N fermion flavours.
In this article, however, we are interested in a specific type of 4-Fermi term that goes beyond the Gross-Neveu couplings (18).These interactions, even for a single fermion flavor, have a non-Abelian SO(5) symmetry.As realised in [38,39], it turns out that there is an exact SO(5) symmetry in the theory of s-wave scattering when working with spin-3/2 alkali gases, similar to the case of 6 Li for the F = 3/2 hyperfine multiplet.The only caveat is that we should consider other atomic species in which the F = 3/2 multiplet S σ ∈ {−3/2, −1/2, 1/2, 3/2} =: {0, 1, 2, 3} corresponds to the lowest-energy level, as the scattering of the higher-energy hyperfine levels can otherwise lead to processes that bring the atoms into the lower hyperfine multiplet [187].There are various possible Alkaline-earth atoms, such as the fermionic isotope 132 Cs, which fulfill this condition and have an F = 3/2 low-energy multiplet.To the best of our knowledge, experiments with Cesium have been reported only for the bosonic 133 Cs isotope, e.g.[188,189].Other possibilities would be to work with the atomic species 9 Be, 135 Ba, 137 Ba, and 201 Hg.
Following our discussion above, the total angular momentum of a colliding pair could be F t ∈ {0, 1, 2, 3} in this case, where Pauli exclusion principle forbids the oddmomentum channels.We thus have a pair of scattering lengths for the singlet a 0 and quintet a 2 channels.The contact interactions can be written in terms of projection operators on these two total angular momenta [187], namely which can be controlled by the individual coupling strengths As shown in [38], these interactions can be rewritten as a linear combination of 4-Fermi terms with a certain definition of gamma matrices.In order to connect to our previous discussion, we first need to define field operators with the right units, such that the spinor operator is , where a is an effective lattice spacing that still needs to be connected to the microscopic cold-atom parameters.Using the corresponding Clebsch-Gordan coefficients ⟨σ 1 σ 2 |F t , M F t , and the specific α-β Dirac matrices ( 14)-( 15), we find where we have introduced the vector α α α = (α 1 , α 2 , α 3 , α 5 ), and the individual coupling constants This brings us already really close to the desired 4-Fermi term in Eq. ( 18), which would require tuning g 2 = g2 , which would require setting a 0 = −13a 2 /11, we would recover exactly the interaction term that combines Thirring, Gross-Neveu and squared mass terms.Let us emphasise, however, that the coupling proportional to g2 , which is proportional to the temporal component of the fermion current J 2 0 , shall not play any role in the phase diagram of the fermionic QFT.Therefore, even if one cannot adjust the scattering lengths of the singlet and quintet channels, the physics will be completely equivalent, at least under the half-filling conditions explored in our work.In fact, as will also be clearer in our discussion below after discussing the cold-atom implementation of the quadratic part, there is a hidden SO(5) symmetry corresponding to a 90 o rotation, regardless of the relative value of these two couplings.This symmetry is responsible for protecting the higher-order topological state discussed in Sec.IV A.

C. Raman lattices and twisted Wilson fermions
Let us now discuss how to realise the twisted Wilson mass regularization in the cold-atom system.In order to obtain the tunnelling structure required by Eq. ( 11), we need to extend Eq. ( 20) by considering additional Raman beams that also assist spin-changing tunnelling processes against certain energy offsets provided by the Zeeman effect of an external magnetic field.This experimental scheme falls within the so-called Raman optical lattices [157,[190][191][192][193][194], which have been exploited as quantum simulators of synthetic spin-orbit coupling [195][196][197].In previous works [93,94,131,198], we highlighted the potential of Raman optical lattices for the quantum simulation of Gross-Neveu-type QFTs with a standard Wilson discretization .For square lattices, these proposals connect to the recent realization of Chern insulators in [157].The goal of this section is to present a Raman-lattice scheme for four spinor components that could serve as a quantum simulator of our SO(5) Dirac field theory regularised with the anisotropic twisted Wilson mass.
Let us first focus on tunnellings that change the atomic states corresponding to the spinor components S σ ∈ {0, 1, 2, 3} by one unit, i.e. σ → σ ′ = σ + 1, where we follow the conventions of Ref. [94].Along the x (y) axis, these tunnellings can be assited by adding Raman beams along the y (x) axis polarized in the x (z) direction [Fig.2(a)].Due to the difference in polarization, the latter, together with the z (x)-polarized standing wave responsible for the standard optical lattice along the x (y) axis, give rise to two-photon spinchanging Raman processes.The key observation is that, due to the different spatial periodicity of the standard lattice and these Raman process, this spin-changing terms cannot contribute with on-site terms, but drive instead an assisted tunnelling.It can be shown that, in the tight-binding limit, this configuration gives rise to spin-changing tunnellings [190][191][192][193] that read where t σ σ +1 j is the corresponding Raman-assisted tunnelling along the j-th axis.Here, we have introduced , where φ σ σ +1 j is the relative phase between the standing wave and the Raman beam.We have also introduced δ σ σ +1 j = ω S − ω σ σ +1 j − (ε σ − ε σ +1 ), which is the corresponding detuning for the two-photon Raman transition, with ω S and ω σ σ +1 j the frequencies of the standing wave and the Raman beam, respectively, and ε σ is the electronic energy for the level σ , which are controlled by the external magnetic field [see the two-photon transitions in Fig. 2(b) and (c)].
In order to realize the lattice field theory described by Eq. ( 11), we need to combine these spin-changing processes.In particular, we need to connect the states (0, 1) and (2, 3) both in the x and y directions, and choose the proper phases φ σ σ ′ j , which can be checked by inspecting the structure of the T 1 and T 2 matrices in Eq. (12).Additionally, we need processes that flip the spinor twice in the y direction, connecting the states (0, 2) and (1, 3).These can be obtained in a similar fashion by using instead y-polarized Raman beams in the x direction [Fig.2(a)], leading to the same expression as in Eq. ( 25) for the case σ → σ ′ = σ + 2, where φ σ σ +2 2 is now the relative phase between two Raman beams and 2 the frequency of the y-polarized Raman beam in the x direction [Fig.2(d)].Let us now detail how the relative phases need to be tuned for each different tunnelling process.
We first focus on T 1 (12), this is, the tunnelling processes along the x axis.The corresponding spin-changing terms can be obtained by using two Raman beams with the same polarization, as explained above, connecting the pairs (0, 1) and (2, 3).The latter can be produced from the same laser source using acusto-optical modulators to generate beams with different detunings δ σ σ ′ j and phases φ σ σ ′ j .Here, we choose in particular δ 01 x = −δ 23 x =: δ x and φ 01 x = φ 23 x = π.After performing the following gauge transformation and rescaling it can be easily checked that the above configuration gives rise to the terms (11), with the following parameters, m1 = δ x /2,t = 1/2a, r = t/t, where we take t := t 01 x = t 23 x .We remark that the required M 1 mass matrix ( 16) is generated by the Ramanassisted tunnelling once we move to the above rotating frame.
Let us now consider the spin-changing processes associated to T 2 (12), for which we also add a spin-independent gradient along the y axis, which can be implemented e.g. by accelerating the optical lattice in that direction [199].This gradient serves a two-fold purpose.First, for ∆ ≫ t, it suppresses the spin-conserving tunnelling in the y direction, which is absent in Eq. ( 11).Additionally, it allows us to tune the relative phase between the two terms in Eq. ( 25), as required by T 2 .Specifically, for each of the four spin-changing pairs of terms involved in T 2 , we employ two Raman beams, and choose the values of the detunings as follows: For every pair of levels, this allows each Raman beam to independently assists one single spin-changing process in Eq. ( 25).The phases of these processes can now be chosen freely, which can be seen by transforming first to the interaction picture with respect to the gradient term in Eq. ( 27), and then applying the rotating-wave approximation in the limit of large detunings.In particular, if we take 11) after applying again the transformations in Eq. ( 26), where we take t 01 y = t 23 y = t and t 02 y = t 13 y = t.Finally, the mass term a 2 ∑ n n n Ψ † n M 2 Ψ n can be obtained by driving transitions between the (0, 1) and (2, 3) spinor pairs using microwave drivings with a Rabi frequency that gives the remaining micorsocopic parameter m2 = Ω y /2.
Let us finally note that the relevant dimensionless parameter that appear in the phase diagrams to be discussed in the rest of the article correspond to where we recall that the Wilson parameter is controlled by the ratio of the tunnellings r = t/t, and the fact that we have neglect the g2 interaction as it will play no role (see the discussion below Eq. ( 24)).The important feature of this mapping is that all of the relevant parameters can be tuned independently in the experiments.As noted previously, the continuum limit does not require sending λ L → 0, but actually working in the vicinity of possible critical lines in parameter space (m 1 a, m 2 a, g 2 /a), which we start to explore below.

IV. HIGHER-ORDER TOPOLOGICAL INSULATORS (HOTIS) AND HIDDEN SO(5) SYMMETRY
In this section, we discuss the regions of parameter space where the free twisted-mass Wilson regularization can lead to HOTIs.In addition of presenting exact expressions for the zero-energy corner modes and the associated topological invariants, we will also discuss the symmetry responsible for protecting these topological states.Although the original SO(5) invariance of the 4-Fermi interactions is explicitly broken by the lattice discretization, we find a hidden discrete SO (5) rotation that protects the groundstate topology.

A. Flat bands and zero-energy corner modes
We now discuss the topological features of the half-filled groundstate in the absence of interactions g 2 = 0.By performing a Fourier transform , where the single-particle Hamiltonian can be written as , 2}, yields a band structure with four energy bands ε q,± (k k k) = ±ε(k k k) that display a two-fold degeneracy This expression allows one to realise that, by setting m1 = m2 = 0 and t = t (i.e.fixing r = 1 and m 1 a = m 2 a = −1), the energy bands become totally flat ε q,± (k k k) = ±2t.This is the consequence of the Aharonov-Bohm phases depicted in Fig. 1, which lead to a destructive interference forbidding the tunnelling of bulk fermions more than two sites apart along any of the spatial directions (see, for instance, the vanishing amplitude of the two black-grey paths in Figs. 1 (b) and (c)).This can be understood as the phenomenon of Aharonov-Bohm caging [120], and finds its minimal manifestation at the corners of the multi-layer, where a single fermion with certain amplitudes over the various spinor components must remain localised.Considering that we have a total of N s = N 1 N 2 lattice sites, with N j sites per spatial direction, we find that the states corresponding to such localised solutions on the corners , and u u u R = (N 1 , N 2 ) have zero energy (see Fig. 3).These zero modes can be expressed as follows As argued in the following section, these corner states are the boundary manifestation of certain topological invariants [200,201] in the bulk bands of the system, which can lead to a topological quadrupole in analogy to the Bernevig-Benalcazar-Hughes (BBH) model [50,51] of HOTIs, in particular to second-order topological insulators.When moving away from the flat-band limit, e.g. by increasing the masses m1 = m2 > 0 or by switching on the quartic interactions g 2 > 0, the perfect Aharonov-Bohm interference will disappear, and these zero modes will no longer be perfectly localised at the corners but, instead, start to spread within the bulk of the system.In particular, a certain localization length will emerge, which characterises the exponential decay of the corner-state amplitude as one moves towards the bulk.These corner states, which remain pinned to zero energy until the bulk energy gap is closed, are an example of anomalous boundary states, which only exist in the presence of the bulk.Moreover, they are protected by a certain hidden symmetry.We previously argued that, in the long-wavelength limit around one of the Dirac points (7), the Euclidean action (8) recovers the SO(3) rotational symmetry of the Lorentz group.In the Hamiltonian formulation, this corresponds to the SO(1, 2) group of boosts, and two-dimensional spatial rotations of angle θ within the xy plane.Let us emphasise that this is a property of the long-wavelength limit, since the lattice regularised action (8) is not invariant under arbitrary rotations On the other hand, there may exist other spatial symmetries that are exact for the full lattice model, including the 4-Fermi interactions (4).For instance, as discussed in more detail in Appendix B, there are two parity symmetries which, at the level of the fermionic field operators, act as follows These transformations correspond to mirror symmetries that take either , and clearly commute with the single-particle Hamiltonian in Eq. ( 29): The composition of the two parities corresponds to the lattice inversion, which is a symmetry of the full Hamiltonian corresponding to a specific SO(1, 2) rotation with angle π.
We remark that the above symmetries need not exhaust all possibilities, as there may be additional spatial symmetries that are not connected to the SO(1, 2) group.The possibility of finding such hidden symmetries becomes clear when inspecting the form of the interacting term (18) which, in fact, allows for generic SO(5) rotations.In the Hamiltonian formulation, these admit a unitary representation, and do not include boosts but only spatial rotations.By simple inspection, it is clear that the first quartic-term ∑ n n n (Ψ † n n n Ψ n n n ) 2 that appears in Eq. ( 18) will be a scalar under any such unitary transformation.In contrast, the remaining terms can be rewritten as ∑ n n n N N N 2 n n n , which involves the norm of a 5-component vector of fermion bilinears Accordingly, this part of the 4-Fermi interaction is also invariant under SO(5) rotations N N N n n n → RN N N n n n ′ , where R is an orthogonal matrix R t R = 1 with detR = 1, and n n n ′ is the corresponding two-dimensional rotation of the lattice.
Although the twisted-mass Wilson regularization (11) breaks explicitly this arbitrary SO(5) symmetry, there can exist specific rotation angles θ and matrices R that correspond to an exact discrete symmetry of the full model.In particular, let us focus on a discrete π/2-rotation that transforms the spatial coordinates as n n n = (n 1 , n 2 ) → n n n ′ = (−n 2 , n 1 ).It is important to emphasise that the action of this rotation on the Dirac spinor Ψ n n n → S R Ψ n n n ′ will not be the same as that of the corresponding SO(1, 2) Lorentz rotation Λ, namely z }, and one can then check that invariance of the lattice model H TM 0 → H TM 0 would require a momentum-independent mass, which is no longer the case with the Wilson-mass regularisation (9).In order to find a π/2-rotation that leaves the model invariant we can, instead, look for a different rotation S R within the above SO(5) group of rotations.This SO(5) invariance requires that the Dirac α and β matrices in Eqs. ( 14)-( 15) transform as We note that the first row of Eq. ( 34) coincides with the transformations that would be obtained from the π/2-rotation of the SO(1, 2) Lorentz group.On the other hand, the second row describes different transformation laws that are crucial to attain invariance of the twisted Wilson mass (11).One can check that Eq. ( 34) has the following solution where S = exp{i π 4 (1 − σ z )} is known as the phase-gate in quantum computing [202], which maps the eigenvectors of σ x onto those of σ y , namely S |± x ⟩ = ± y .Combining these transformations with the action of the rotation on the crystal momenta, we find that the Bloch Hamiltonian ( 29) is indeed invariant when . At the level of the zero modes (31), the spatial part of this hidden SO(5) transformation respects the set of corners, and one says that these anomalous boundary states are protected by this hidden spatial symmetry.Since they only have support on a region of codimension 2, they are the boundary manifestation of symmetry-protected HOTI groundstates.We note that this symmetry can be interpreted as the multi-layer counterpart of the C 4 symmetry of the BBH model [50,51].
Before closing this subsection, we remark that the unitary transformation on the spinors Ψ n n n → S R Ψ n n n ′ , together with the transformation of the Dirac α and β matrix (34), can be used to show that the vector of bilinears simply transforms as It is then clear that its norm is conserved and, thus, the quartic interactions (18) are also left invariant under this hidden SO( 5) rotation H int → H int .Altogether, this proves that the full lattice model is invariant H = H TM 0 + H int → H. Since the set of all four corners is also invariant under a π/2 rotation, we expect that the anomalous corner states will be protected by this symmetry and, thus, robust when varying the microscopic parameters of the model and switching on interactions.The only possibility to get rid of them is by closing the bulk energy gap, which would signal a quantum phase transition to a trivial band insulator or to a groundstate with symmetrybroken long-range order.We will explore these possibilities in the sections below but, first, let us provide a bulk perspective of the higher-order topology associated to these corner modes.

B. Higher-order topological invariants
According to our current understanding of the bulkboundary correspondence of symmetry-protected topological phases [118], the anomalous zero modes are a boundary manifestation of a topological band structure in the bulk.For the present model, this phase should be a HOTI od second order with a certain non-vanishing topological invariant.As shown in [200,201], this invariant can be expressed as the product of two winding numbers, which will allow us to find the phase transitions between the HOTI and a trivial insulator for g 2 = 0.
At the level of the twisted-mass free Hamiltonian (11), one sees that β H 0 (k k k)β = −H 0 (k k k).This corresponds to the AIII class in the classification of topological insulators under global symmetries, and shows that the band structure (30) always comes in pairs of positive-negative energies.The idea now is to find an alternative unitarily-equivalent representation of the Dirac matrices that intertwines the effect of this AIII symmetry with the dimensionality of the problem.This is, once more, only allowed by the fact that we are working with a reducible representation of the Clifford algebra.By applying the following unitary transformation U, built again from the S-gate, one finds that the β -matrix transforms as where β 1 = σ z and β 2 = σ y .On the other hand, the Dirac α matrices transform according to the following expressions which can be used to show that the transformed Bloch Hamiltonian ( 29) has the following tensor-product structure Here, we have defined the following single-particle Hamiltonians that only depend on the kinetic energy and Wilson mass along a specific direction in momentum space One can check that each of these Hamiltonians has an individual AIII symmetry β 1 h 1 (k 1 )β 1 = −h 1 (k 1 ), and , which guarantees that the full Hamiltonian fulfills the desired transformation β H 0 (k k k)β = −H 0 (k k k).Each of these terms (40) can be understood as a two-band model corresponding to an un-twisted Wilson regularization of a (1 + 1)-dimensional Dirac-fermion QFT [129].This tensor-product construction also highlights that the corner states are not equivalent to the edge states of AIII topological insulators of one-dimensional chains, arranged along the boundaries of a square lattice.It is really the two-dimensional bulk that is required to host and protect these corner modes.
The lower-dimensional band structures of h j (k j ) have a pair of Dirac points corresponding to the projection of the previous Dirac points (7) onto the respective axis k D,ℓ j = k k k D,ℓ ℓ ℓ • e j , each of which presents a different mass In analogy to our discussion of the Chern insulator (C6) an the mass matrix (C4), we can now define two mass matrices each of which contains the information about each of the twisted Wilson masses at the corresponding projection of the Dirac point.The Berry connection for each of these projections is We represent the topological invariant in Eq. ( 45) as a function of the bare twisted masses m 1 and m 2 .In the green inner square m j a ∈ (−2r, 0), the topological invariant is non-trivial e iπW 1 W 2 = −1, and the groundstate of the twisted-Wilson lattice model corresponds to a higher-order topological insulator (HOTI).The shaded line with m 1 = m 2 represents the regime where the model has a hidden SO (5) symmetry that protects the corner modes.The white region where m j a / ∈ (−2r, 0) for one or both masses corresponds to a trivial band insulator (TBI) with trivial topological invariant e iπW 1 W 2 = +1.
the negative-energy modes that would be filled in the corresponding lower-dimensional groundstate.One can define a Chern-Simons form [80] associated to this Berry connection or, equivalently, a so-called Zak's phase [203], which plays the role of the above Chern number (C6) in this reduced dimensionality.We find that each of these invariants is again non-trivial when an odd number of the projected Dirac points have a negative twisted Wilson mass As is well-known for standard first-order topological insulators [204], these topological invariants CS j = W j /2 are proportional to the winding number W j of the mappings d , where d d d j (k j ) is the vector of coefficients of the individual Hamiltonians (40) in the Pauli basis h j (k j ) = d d d(k j ) • σ σ σ .As discussed in Reference [80], one can define a topological invariant that is also gauge invariant by considering the Wilson loop associated to such winding number.In this way, by simply multiplying the winding numbers together and exponentiating them, we obtain a topological invariant for HOTIs with AIII symmetry (39), which reads A non-vanishing invariant e iπW 1 W 2 = −1 signals the nontrivial topology of the bulk, and must have a boundary manifestation in the form of corner states.Indeed, for m 1 a = m 2 a = r, which corresponds to the previous flat-band limit m1 = m2 = 0 when r = 1, we find that e iπW 1 W 2 = −1 in the bulk, while the zero-energy corner states are those of Eq. (31).Away from this flat-band limit, and while m j a ∈ (−2r, 0) in both directions, the groundstate is still a HOTI with e iπW 1 W 2 = −1, and the anomalous boundary states remain exponentially localised to the corners.The non-interacting HOTI in parameter space corresponds to the square displayed in Fig. 4. For mass parameters in the inner square, the groundstate is a HOTI whereas, outside this region, it is trivial.
In the next section, we will explore the fate of this HOTI as the fermion self-interactions g 2 are increased.Let us note that, in our D = 3-dimensional spacetime, considering the role of the quartic interactions by simple power counting can be misleading.Indeed, the dimensions of the Dirac field is whereas the coupling strength has units of length [g 2 ] = L and, thus, a negative energy dimension.Naive power-counting arguments would then suggest the 4-Fermi interaction is perturbatively irrelevant, such that the QFT in D = 2 + 1 would be nonrenormalizable.On the other hand, it is well-known that the relevant/irrelevant nature of the couplings can be modified after resummation in the large-N limit [205], whereupon the Thirring-like interaction Even if the corresponding 4-Fermi field theories are not perturbatively renormalizable, they become 1/N renormalizable [8,[23][24][25][26].The goal of the following section is to analyse the effect of the 4-Fermi interactions on the HOTI, including also the mass-squared interactions ∑ j=3,5 (Ψ † n n n α j Ψ n n n ) 2 , using the large-N expansion.

V. CORRELATED HOTIs, FERMION CONDENSATES AND QUANTUM PHASE TRANSITIONS
In this subsection, we discuss the effect of the 4-Fermi interactions in detail.We argued previously that the HOTI groundstate is protected by a hidden SO(5) symmetry and, thus, should be robust under symmetric perturbations unless those are sufficiently strong such that the bulk energy gap closes allowing for a change of the topological invariant, or if a certain symmetry-breaking phase transition takes place.We also showed in Eq. ( 36) that the 4-Fermi interactions preserve this hidden symmetry, such that one expects the HOTI phase to become a correlated HOTI as one increases the coupling strength g 2 > 0. Eventually, when the interactions are sufficiently strong, there may be a symmetry-breaking phase transition at some g 2 c , which paves the way for the appearance of new phases of matter typically referred to as fermion condensates in the QFT literature.The goal of this section is to explore the possible condensates allowed by the rich SO(5) structure of the self-interactions, and provide a quantitative account about which of the fermion condensates is expected to form at which point in parameter space (m 1 a, m 2 a, g 2 /a).

A. Auxiliary fields and the hidden symmetry
To accomplish this goal, we shall return to the Euclidean formulation of our SO(5) Dirac matter in Eqs. ( 1) and ( 4), where we can make use of controlled approximations such as the large-N expansion [205].To present a non-perturbative, yet tractable, account of the 4-Fermi interactions, one can generalise the Euclidean action in Eqs. ( 1) and ( 4) to N flavours The free part of the regularised Wilson twisted-mass action S TM 0 (8) simply becomes a sum of the corresponding actions for each of the fermion flavours.In contrast, the 4-Fermi term does couple the different flavours where, to have a consistent N → ∞ limit, one must rescale the coupling strength as where the gamma matrices for N flavours appearing the in the above bilinears should be understood as γ a → 1 N ⊗ γ a .We will assume this in all the expressions below, which leads to an additional U(N) symmetry in flavour space.
The first step of the large-N approximation is to introduce auxiliary fields via a Hubbard-Stratonovich transformation [206,207], such that the partition function becomes quadratic in the Grassmann spinors.We need to introduce six real bosonic fields a µ (x), σ 1 (x), σ 2 (x), and π(x), such that the partition function can be exactly rewritten as 2g 2 (aµa µ +σ j σ j +π 2 ) .
(48) As a consequence of this transformation, we get the following action coupling between the fermionic and bosonic fields (49) Except for the π-field, the rest can be incorporated in the free action S TM 0 by the following substitution which shows that the auxiliary fields a µ (x) act as an effective gauge-like potential, the components of which admix by the Lorentz transformations, and couple to the fermions minimally.We stress, however, that the free action (48) of these auxiliary fields is not gauge invariant, but simply a masslike term that becomes very heavy in the large-N limit.On the other hand, the σ j (x) are scalar fields that couple to the fermions via a pair of twisted Yukawa-type couplings.To be more accurate, we should consider the Euclidean formulation of the two parity transformations (32), which amount to Since the 4-Fermi interactions are invariant under these parity transformations, we know that S ′ int → S ′ int , which require the auxiliary gauge-like fields to transform as whereas the auxiliary σ fields transform as Therefore, we see that the σ fields are all parity even, the a 0 (x) component of the gauge-like field is parity even, while the spatial components a a a(x) are parity odd.
We have so far left out the discussion of the π(x) field, as it gives rise to a new term that was not present in S TM 0 .Considering the parity transformations in Eq. ( 51), we see that in order to recover parity invariance, the π(x) field should transform as and, thus, corresponds to a pseudo-scalar field that is parity odd.According to Eq. ( 49), this pseudo-sclalar auxiliary field couples to the fermions via the standard Yukawa coupling.Let us finally discuss the hidden SO( 5) symmetry ( 35) which, for the Euclidean Grassmann fields, corresponds to where we recall that S R is the unitary matrix given by Eq. (35).Considering that the Euclidean gamma matrices transform as one can readily see that the original Euclidean action Eqs.( 1) and ( 4) is also invariant under this hidden rotational symmetry (55) when m 1 = m 2 .Once the auxiliary fields are introduced, this hidden symmetry implies that the vector field should transform as whereas the scalar and pseudo-scalar fields must fulfill One can check that these pair of equations ( 57) and ( 58) define a transformation on a vector of auxiliary fields φ φ φ (x) := (a 0 (x), a 1 (x), a 2 (x), σ 1 (x), σ 2 (x), π(x)) → Oφ φ φ (τ, x 2 , −x 1 ), with O t O = 1 and detO = +1.We could then say that, within the Euclidean formulation where the field and the adjoint are independent Grassmann fields, the hidden symmetry (55) can be interpret as a specific SO(6) rotation of the auxiliary fields and fermion bilinears.We note that the interacting part of the Euclidean action has indeed a larger symmetry under arbitrary SO(6) rotations when expressed in term of auxiliary fields, which gets broken down by the lattice regularization of the free fermions.It is only for the specific π/2 rotation above that one recovers invariance of the full Euclidean action.However, it must be noted that this transformation is a rotation of the auxiliary fields about the axis a 0 (x) → a 0 (x ′ ).
If one considers the adjoint definition in the Hamiltonian approach, a 0 (x) should always remain invariant, such that the hidden symmetry reduces to the previous SO( 5) rotation.
The idea of the large-N method in the present context is to assume that these scalar fields will be homogeneous in the groundstate of the interacting theory φ φ φ (x) = Φ Φ Φ, ∀x, and try to determine the regime in parameter space (m 1 a, m 2 a, g 2 /a) where some of them achieve a non-zero vacuum expectation value.This is the region of parameter space where the groundstate supports a specific combination of fermion condensates Each of these fermion condensates is responsible for breaking a particular symmetry, and may even change completely the QFT that governs the continuum limit in the vicinity of such a symmetry-breaking phase transition.This is the case of the vector condensate, which is proportional to a fermion current A µ ∝ ⟨J µ ⟩, and would thus forbid recovering Lorentz invariance even in the continuum limit.Such condensates have been identified in related two-band models [93,208].
The vacuum expectation values of the scalar and pseudoscalar condensates play a different role.In fact, the scalar condensates Σ 1 ∝ ⟨ψiγ 3 ψ⟩ and Σ 2 ∝ ⟨ψiγ 5 ψ⟩ are generally non-zero except for a particular line in parameter space m 1 a = m 2 a = −r, which is a consequence of the twisted Wilson mass regularization.These scalar condensates do not break any of the parity symmetries (32), but contribute with a renormalization of the bare twisted masses m j → m j + Σ j which, as will be discussed below, can change abruptly the value of the topological invariant (45).When the bare masses m 1 = m 2 =: m, the two scalar condensates take equal values Σ 1 = Σ 2 =: Σ, such that the hidden SO(5) protecting symmetry ( 58) is not broken, and one can still talk about the symmetry-protected HOTI.By increasing the interaction g 2 , as discussed below, one of the possibilities is that the values of Σ will change and lead to an interaction-induced quantum phase transitions between the correlated HOTI, and a trivial band insulator with no corner states and a vanishing manybody topological invariant.
Before closing this section, we comment on the remaining fermion condensate Π ∝ ⟨ψψ⟩.Although in even spacetime dimensions, this condensate is parity even and associated to the breakdown of chiral symmetry, in our even-dimensional QFT it plays a different role.As discussed above, this condensate is odd under any of the parities (32), and a finite vacuum expectation value would imply the spontaneous breakdown of parity.The possibility of finding such condensates in the standard Wilsonian lattice regularization of Dirac QFTs was initially considered by S. Aoki [209], and it is typically referred to as an Aoki condensate in lattice gauge theories.In our current model of HOTIs, rather than parity breaking, it is more important to consider the hidden SO(5) symmetry, which is spontaneously broken by a non-zero value of this π condensate (58).We note that the formation of vector condensates A A A ̸ = 0 0 0 would also break spontaneously the protecting symmetry in light of Eq. ( 58).As shown below, understanding the competition of the different condensation channels is the key to understand the phase diagram of our HOTI.

B. Large-N condensates and the effective potential
In this subsection, we describe the results of the aforementioned large-N technique to chart the phase diagram of the interacting HOTI.As advanced previously, there are various possible fermion condensates characterised by different vacuum expectation values, which could be obtained in the N → ∞ limit by solving a set of gap equations.As discussed in the context of lattice gauge theories [210], these gap equations are, however, only valid for non-vanishing values of the vacuum expectation values Φ a ̸ = 0. On the other hand, we are also interested in the competition of the HOTI with a trivial band insulator, where the symmetry-breaking fermion condensates vanish and there is no spontaneous symmetry breaking.In order to explore the whole phase diagram, we need to go beyond the gap equations and obtain the large-N effective potential V eff (Φ Φ Φ), the minimum of which will provide the values of the auxiliary fields in Φ Φ Φ for any specific point in parameter space determining, in particular, which of the possible symmetry-breaking condensates prevail.
The large-N effective potential can be obtained diagrammatically by considering the Feynman diagrams with a single fermion loop and an increasing number of external lines describing the auxiliary fields.Any other one-particle irreducible diagram with more fermion loops and internal auxiliary lines contributes with a higher order in 1/N, and can thus be neglected when N → ∞.In the standard calculation of the chiral-invariant Gross-Neveu QFT, one can show that only diagrams with an even number of external auxiliary legs can give a non-zero contribution [205].On the other hand, for our twisted Wilson mass regularization, one cannot apply the same arguments, and must also take into account the diagrams with an odd number of external legs.A similar situation can also be found in Chern-insulator models with a standard Wilson discretization although, there, one finds that these odd terms are zero due to the vanishing of the corresponding integrals [93,94].For the present model, this is not the case, and both the even and odd Feynman diagrams have a non-zero contribution that must be accounted for (see Appendix D).
The resummation of these Feynman diagrams yields the leading-order quantum radiative corrections δV eff (Φ Φ Φ) to the classical potential of the auxiliary fields As we have assumed that these auxiliary fiels are homogeneous, all the relevant information is included in this effective potential, which plays a crucial role in determining the parameter regimes where the groundstate displays non-zero vacuum expectation values Φ Φ Φ ̸ = 0.In fact, the classical part of the potential predicts a zero vacuum expectation value Φ Φ Φ = 0, and it is the contribution of radiative corrections δV eff (Φ Φ Φ) for g 2 > 0 that can change the minimum of Eq. ( 60) allowing for condensation.As discussed in Appendix D, this resummation can be accomplished to all orders of g 2 , which allows us to address non-perturbative effects in the phase diagram of the model.We find that the quantum correction can be expressed as where we recall that k k k is the regularised spatial momentum in Eq. ( 6), and we have introduced the shifts in the twisted Wilson masses (9) stemming from the additive renormalisations that one finds for a non-zero values by the scalar condensates We also recall that the integral symbol is a short-hand notation (5) for the spatial mode sum and the zero-temperature limit of the Matsubara sum.As discussed in Appendix D, the temporal component of the pseudo-vector field A 0 does not contribute to the effective potential and, thus, cannot condense.The occurrence of other condensation channels will depend on the minimum of V eff (A A A, Σ Σ Σ, Π).For m 1 = m 2 , this effective potential can be easily seen to be invariant under the hidden SO(5) rotation that takes Let us note that the above expressions with the specific lattice regularization are in fact more general, and would also apply to continuum QFTs with the same 4-Fermi interactions.This would simply require substituting k k k → k k k and m j (k k k) → m j , recovering in this way the underlying Lorentz invariance.On the other hand, the expressions could also be used with a discretized Euclidean time by substituting k 0 → sin(k 0 a)/a with k 0 = −π/a + 2π(n 0 + 1/2)/aN 0 with n 0 ∈ Z N 0 .From the perspective of the cold-atom quantum simulator, one is interested in the continuum-time limit and the Hamiltonian field theory.In this case, one can actually perform the integral over k 0 ∈ R, and express the radiative corrections as follows where we have introduced a short-hand for the spatialmomenta sum within the Brillouin zone k k k = 1 (N s a) 2 ∑ k k k∈BZ .In addition, we have generalised the energy dispersion relation in Eq. (30) to account for the possible non-zero vacuum expectation values of the fermion condensates Expression ( 63) has a very simple interpretation, the large-N radiative corrections are given by the total shift of singleparticle energy levels in the filled bands, i.e. those with nega-tive energies forming the Dirac sea, once some of the condensates form Φ Φ Φ ̸ = 0 0 0. For the scalar condensates Σ j , such corrections appear as soon as the interactions are switched on g 2 > 0.
The only exception is the straight line at m 1 a = m 2 a = −r, in which the scalar condensates vanish by symmetry arguments.For the vector A A A and pseudo-scalar Π condensates, the situation is completely different.A non-zero value of these condensates would break the hidden SO(5) symmetry in Eqs. ( 57)-( 58), which can only happen spontaneously for a sufficiently-strong coupling g 2 > g 2 c (m 1 , m 2 ).In order to find out which of the condensates prevails, we need to minimize the full effective potential which is an unconstrained multi-parameter non-linear minimization problem that must be addressed numerically, as we detail in the following subsection.

C. Self-energy and correlated HOTIs
In the previous subsection, we have discussed the procedure to find the values of Φ Φ Φ ⋆ , which will allow us to detect the symmetry-breaking condensates and localise the critical lines that separate them from the HOTI and the trivial band insulator.On the other hand, we would also like to predict the flow of the critical lines separating the HOTI from the trivial band insulator as the coupling increases g 2 > 0. Since topological phases cannot be distinguished by a local order parameter, we also need to calculate the topological invariant (45) away from the non-interacting free theory.An approximation that has already been used for the many-body topological invariants of standard topological insulators [211][212][213][214][215], deals with the so-called topological Hamiltonian.
Many-body topological invariants can be defined via the two-point Green's functions G(x 1 − x 2 ) = ⟨T {Ψ † (x 1 )Ψ(x 2 )}⟩ [216,217], where the expectation value is calculated over the groundstate with non-zero interactions.Following the prescriptions of quantum-many body physics within condensed matter [75], by going to momentum space, the inverse of the Green's function can be expressed as where H 0 (k k k) is the single-particle Hamiltonian, the N-flavour version of Eq. ( 29) in our case, and Σ s (ik 0 , k k k) is the socalled self-energy.This self-energy contains all the oneparticle irreducible "tadpole" contributions to the fermion propagator arising from intermediate scattering processes in which particle-antiparticle pairs are virtually created from the groundstate.Within our large-N theory, we have precisely calculated those at leading order in N by introducing the auxiliary fields.In fact, the above condensates can be readily used to approximate this self-energy as which has no momentum dependence Σ s (ik 0 , k k k) = Σ s (0, 0 0 0) since we have assumed the condensates to be homogeneous.As discussed in [211][212][213][214][215], the static contributions to the self-energy Σ s (0, k k k) can be used to define the so-called topological Hamiltonian.To consider the same symmetry class as the non-interacting one, we set A A A = 0 0 0 and Π = 0, such that the hidden SO(5) rotation is preserved.In this case, the only contribution to the topological Hamiltonian stems from the scalar condensates Within this large-N approximation, the many-body topological invariant can be expressed in terms of the two Chern-Simons forms (69) with the twisted Wilson masses renormalised by the two scalar condensates (62), namely The full topological invariant of the correlated HOTI can be approximated, within the large-N limit, as follows Accordingly, the large-N solution obtained by minimizing the effective potential ( 65) can be readily used to extract Σ 1 , Σ 2 in the parameter regime where the hidden SO(5) symmetry is still preserved, and localise the critical surface that separates the correlated HOTI from the trivial band insulator.We recall again that, in order to have the symmetry protection, we need to consider the regime where m 1 = m 2 and Σ 1 = Σ 2 , which will cut this critical surface determining a critical line.

D. Large-N results and phase diagram
We start by making the assumption that only the pseudoscalar auxiliary field condenses Π = ⟨π⟩ ̸ = 0, which we expect describes the leading symmetry-breaking channel among all other condensates, and should set in at a certain value of the interaction strength g 2 > 0. This assumption will allow us to derive a set of simple gap equations that can be later used as a reference in the minimization of the full effective potential (65) that contains all possible condensation channels (63).As with gap equations in other models [93,94,131], this formalism is applicable whenever the order parameter is non-zero Π > 0 [210].As advanced in the previous sections, the σ condensates typically attain non-zero values for any point in parameter space (m 1 a, m 2 a, g 2 /a), except for the line at fixed twisted mass m 1 a = m 2 a = r.We hence consider that Σ Σ Σ = ⟨σ σ σ ⟩ ̸ = 0 0 0. Therefore, the main assumption in the following gap equations is that the vector condensate vanishes A A A = ⟨a a a⟩ = 0 0 0 (this will be justified in Fig. 7 below).The gap equations are derived by solving the set of nonlinear equations given by the stationary point ∂ Π V eff | A µ =0 = FIG. 5. Interacting higher-order topological phase diagram I: We consider parameter space (m 1 a, m 2 g 2 /a), and represent the different phases in the twisted-Wilson lattice model with 4-Fermi interactions.The blue and magenta curves (74) delimit the inner green area (73) at fixed-g 2 slices that has a non-trivial many-body topological invariant e iπW 1 (Σ Σ Σ)W 2 (Σ Σ Σ) = −1 for an odd number of fermion flavours.The red lines are obtained by solving the gap equations to localise the boundary of the parity-breaking Π condensate ( 71)- (72).The black lines demonstrate that this pseudo-scalar condensate actually grows from the four corners of the non-interacting HOTI phase at g 2 = 0 (see Fig. 4), forming the shape of an Eiffel tower that rests on the correlated HOTI.All numerical evaluations employed a spatial lattice with N s = 32 2 sites.
Setting the Wilson parameter to r = 1 henceforth, we find where we have assumed that Π > 0 to simplify the equation.In addition, we get the following two equations from the derivatives with respect to the scalar condensates where we have used Eq. ( 71) to simplify them, getting an expression that relates to the twisted masses m 1 , m 2 .Equations ( 71)-( 72) can be solved with Π = 0 + to map out the boundary of the phase hosting a pseudo-scalar condensate.This critical region is a surface in parameter space (m 1 a, m 2 a, g 2 /a), and is shown in red in Fig. 5.This surface is spanned by a number of trajectories that represent solutions to the gap equations for Π = 0 + , for which we fix a different value of the twisted masses renormalised by the scalar condensates (m j + Σ j ).These trajectories are plotted in Fig. 5 as a collection of red dashed curves; the resulting network visualises a closed surface which descends from strong coupling down to g 2 /a ≈ 0.7.There are a couple of interesting comments: (i) The volume inside the red surface describes the spontaneous symmetry-broken phase with a pseudo-scalar Here the lines where the gap functions f 1 = 0 (blue) and f 2 = 0 (magenta) are plotted for increasing values of the coupling strength g 2 /a ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7}, which grow as one moves towards the center of the square.The HOTI phase is contained within a star-shaped region that is delimited by groups of four of these lines, corresponding to the same value of the coupling strength g 2 /a.condensate.This condensate breaks any of the parities, and corresponds to the so-called Aoki phase found in other lattice field theories [209].As emphasised above, the more important thing is that this pseudo-scalar condensate also breaks spontaneously the hidden SO(5) symmetry responsible for the protection of the HOTI.Hence the higher-order topological invariant an the corner modes cannot coexist with this pesudo-scalar condensate.(ii) The volume that contains this SO(5) breaking condensate is centered around the symmetry line m 1 a = m 2 a = −1 in which the scalar condensates vanish Σ Σ Σ = 0 0 0, and gets more and more compressed as the interaction strength becomes large g 2 /a ≫ 1.On the other hand, for g 2 = 0, the gap equation ( 71) becomes singular at the four corners (m 1 , m 2 ) ∈ {(0, 0), (0, −2), (−2, 0), (−2, −2)}, and we have found that this phase extends all the way down to weak coupling in four very sharp spikes.Let us note that this is a large-N prediction, and a different method should be used to determine the extent of this phase for finite N.
As already remarked above, the gap equations ( 71)-( 72) are only valid for Π > 0. On the other hand, the original question that motivated the study was to see the extent of the HOTI as one increases interactions, which would require exploring the region with g 2 > 0 for which Π = 0, more particularly, the weaker-coupling regime beneath the Aoki phase.To explore it, we will directly minimize the effective potential V eff (0 0 0, Σ Σ Σ, 0) at a specific (m 1 a, m 2 a, g 2 /a).As discussed in Eq. ( 70), the many-body topological invariant for the HOTI in the large-N limit is non-trivial e iπW 1 (Σ Σ Σ)W 2 (Σ Σ Σ) = −1 when and we have an odd number of flavours.The correlated HOTI phase will then be contained in a region bounded by contours along which a gap function f j (k j a, Σ j a) = m j a + (1 − cos k j a) + Σ j a (with j = 1, 2) vanishes at either the origin k j = 0, or at the zone edge k j = π/a.The procedure is then to search for solutions of the non-linear equations using the value of the scalar condensates Σ 1 , Σ 2 obtained by numerically minimisation of V eff (0 0 0, Σ Σ Σ, 0).In practice, we define a circle centred at the line of symmetry where Σ Σ Σ = 0 0 0 by defining the twisted masses as (m 1 , m 2 ) = (−1/a, −1/a) + m(cos θ , sin θ ), for m > 0 and θ ∈ [0, 2π).We then search for the roots of Eq. ( 74) by scanning first in θ and, subsequently, in m.The resulting contours along which a gap function vanishes f j (k j , Σ j ) = 0 are plotted as blue (magenta) lines in Figs. 5 and6 for various values of the coupling strength g 2 .Each pair of blue (magenta) lines for a fixed coupling strength corresponds to the renormalised twisted mass proportional to iγ 3 (iγ 5 ) satisfying f 1 (0, Σ 1 ) = 0 or f 1 (π/a, Σ 1 ) = 0 ( f 2 (0, Σ 2 ) = 0 or f 2 (π/a, Σ 2 ) = 0).The region of the correlated HOTI phase is enclosed within the areas inside the four intersecting lines, and is depicted in Fig. 5 by a shaded green area that connects to the green square in the non-interacting limit (Cf.Fig. 4), and projected onto the (m 1 , m 2 )-plane in Fig. 6.As g 2 increases, the borders of this region curve inwards, and the HOTI shrinks until it roughly coincides with the lateral extent of the Aoki phase at the critical coupling g 2 /a ≃ 0.7 (note the lower surface of the Aoki phase has convex curvature, confirmed in Fig. 7).We remark that the hidden SO(5) symmetry responsible for the protection of the anomalous corner modes takes place within the m 1 = 2 area of the green volume of Fig. 5. On the other hand, all the empty region surrounding both the green and red volumes corresponds to a trivial band insulator, in which the topological invariant is trivial and the pseudo-scalar condensate vanishes.All the critical surfaces predicted within our large-N methods correspond to higher-order quantum phase transitions, either topological or symmetry-breaking ones.
Once we have discussed the phase diagram of Fig. 5 in detail, we should check for the consistency of the assumptions we made about the competing condensates by minimising the full effective potential V eff (A A A, Σ Σ Σ, Π).We recall again that one can set A 0 = 0 as discussed in Appendix D, but must consider the other competing condensation channels along lines of fixed (m 1 , m 2 ).Fig. 7 shows the resulting condensates for two choices of (m 1 , m 2 ) at different distances from the line of symmetry m 1 = m 2 = −1/a.The Π condensate signalling the SO(5)-breaking Aoki phase rises from zero at a critical coupling g 2 c /a ≈ 0.7, and the opposite signs of the scalar condensates Σ 1 = −Σ 2 correspond to (m 1 , m 2 ) lying on opposite sides of the symmetric line (−1/a, −1/a).Crucially, the current condensate A 1 remains zero throughout, justifying our previous assumption where we set A A A = 0 0 0. Solutions with A A A ̸ = 0 0 0, as found for an interacting Chern-insulator twoband model in [93,94], could only be found in our lattice model by artificially constraining Π 0. Otherwise, we find that the SO(5)-breaking Aoki condensate always describes the leading condensation channel Π > 0 as one increases the coupling strength.Further from the line of symmetry, as shown by the dashed lines in Fig. 7  same, but with values of the scalar condensates Σ i that become larger in magnitude; this time at strong-enough coupling the trajectory actually re-enters the SO(5)-symmetric phase, reflected by the renewed vanishing of Π and the kinks in Σ i (g 2 ).however, rather than re-entering into the HOTI phase, one goes into a trivial band insulator where, even if the SO(5) symmetry is preserved, the topological invariant is trivial and the zero states are no longer localised at the corners of the system.

VI. CONCLUSION AND OUTLOOK
In this work, we have presented a non-standard lattice regularization of Dirac QFTs based on a new type of Wilson fermions that have an anisotropic twisted Wilson mass.We have shown that the anisotropic twisted Wilson mass is responsible for the occurrence of HOTIs that display zeroenergy corner modes, and a non-vanishing topological invariant in the bulk.We have discuss a cold-atom implementation of this lattice field theory that exploits Raman optical lattices, and spin-3/3 Fermi gases of alkaline-earth atoms.Interestingly, the s-wave scattering of this atoms leads to a SO(5) invariant 4-Fermi interaction, which leads to an interesting competition of HOTI phases and various fermion condensates.We explored the full phase diagram of the model using large-N techniques, and argued that the correlated HOTI phase eventually gives way to a parity-breaking fermion condensate.
Since the microscopic parameters of the model can be independently tuned in the proposed cold-atom experiment, it woudl be very ineteresting that the predictions presented in this work could be tested in future experiments.We note that the recent experimental work on the quantum simulator of Chern insulators using 87 Sr Fermi gases in Raman optical lattices [157] is very promising in this direction.If one can use these methods for a different alkaline-earth atoms such as 132 Cs, controlling the four spin states as discussed in our work, the experiment would directly probe the SO(5) self-interacting Dirac QFT that has been the subject of our work.This experiment, together with related theory works [218,219], have also shown that it is possible to perform certain measurements to infer the value of a topological invariant similar to the Chern number.It would be interesting to explore if such methods can also be adapted to the spin-3/2 Fermi gas, and used to infer the value of our higher-order topological invariant.Other than that, further studies are required in order to propose other measurement schemes.For instance, in order to retrieve the order parameter associated to the fermion condensates discussed in this work, some sort of spin-resolved imaging by either illuminating the gas and processing its shadow, or using quantum gas microscopes, should be required.These methods should be combined with microwave transition and spin-selective techniques to infer the atomic densities corresponding to the order parameters.
From a more theoretical perspective, it would be very interesting to explore finite-temperature and finite-density phases in this model of SO(5) interacting Dirac matter.In particular, by moving away from half filling, one should consider other possible condensation channels that likely include superconducting and inhomogeneous orders, or even new exotic orders that go beyond the Landau symmetry breaking paradigm.The study of those drawing further connections between highenergy physics, condesed matter, and AMO physics, will contribute to the growing interest in this interdisciplinary line of research [220][221][222][223][224][225][226][227][228].
leads to a crucial difference with respect to Minkowski spacetime, and forbids defining the adjoint as ψ(x) = ψ † (x)γ 0 .The field and its adjoint are mutually anti-commuting Grassmann spinors with an even number of components d s , and one postulates that they transform under spacetime rotations as One then finds that the Euclidean action (A1) is invariant under SO(D), which also requires using the transformations of the rank j < D Clifford elements as tensors under SO(D) This shows that the gamma matrices transform as a vector under SO(D), such that ψψ and ψγ a ∂ a ψ are scalars, and the Euclidean action is invariant under SO(D).
For even spacetime dimensions D = 2n, the highest-rank element of the Clifford algebra (A3) also plays an important role.It can be used to define an additional Hermitian matrix that anti-commutes with all the spacetime gamma matrices, and is thus left invariant under any SO(2n) rotation In D = 4, this is typically called the chiral gamma matrix γ 5 = −γ 0 γ 1 γ 2 γ 3 , which can be used to decompose the Dirac spinor into left-and right-handed two-component spinors, the so-called chiral Weyl fermions [1].Alternatively, in any even spacetime dimension, this gamma matrix can serve to propose a twisting of the scalar mass where m 1 = m cos θ , and m 2 = m sin θ .The anti-commuting mass terms can be expressed as which respects Lorentz SO(2n) invariance according to Eqs. (A4) and (A6), and can be seen as the result of an axial rotation ψ → exp{i θ 2 γ ⋆ }ψ, ψ → ψexp{i θ 2 γ ⋆ }.On the other hand, the second one breaks explicitly the parity symmetry since For odd spacetime dimensions D = 2n + 1, this γ ⋆ matrix plays the role of the gamma matrix for the new spatial direction γ 2n = γ ⋆ .Therefore, the product of all spacetime gamma matrices is trivial γ 0 γ 1 • • • γ 2n ∝ 1 d s , and the SO(2n + 1)invariant Dirac action for free Dirac fields can only take the form of Eq. (A1).Therefore, only the standard mass term m ψψ can be considered which, as discussed in the following section would break the invariance under the corresponding parity transformation.In the following subsection, we will explain how to go beyond these limitations when considering a reducible representation of the Clifford algebra.Motivated by the experimental situations discussed in the main text, we can also consider reducible representations of the Clifford algebra for a specific spacetime dimension.In this appendix, we will consider odd dimension D = 2n + 1, and understand a reducible representation of the Clifford algebra as a consequence of an effective dimensional reduction.We consider 2n+3 dimensions initially, where the spinor dimension is doubled with respect to the irreducible one of D = 2n + 1, and we get two additional gamma matrices γ 2n+1 , and γ ⋆ = γ 2n+2 .We will label the higher-dimensional spacetime coordinates with latin indexes a ∈ {0, 1, • • • , 2n + 2}.The SO(2n+3)-invariant action is that of Eq. (A1), and we will focus in the massless case m = 0, namely We can rewrite this action by separating the contribution of the two extra spatial dimensions B2) where the index µ ∈ {0, • • • , 2n} is restricted to the lower number of dimensions D = 2n + 1 in the reduced spacetime.
As discussed in [80], the dimensional reduction is inspired by the Kaluza-Klein compactification [79], and proceeds in two steps.In the first one, the x 2n+2 spatial direction is compactified to a circle x 2n+2 + r = x 2n+2 with a very small radius r → 0. Considering that the Grassmann fields are periodic in the spatial direction, the corresponding momentum p 2n+2 = −i∂ 2n+2 ∝ ℓ 2n+2 /r gets quantised in terms of the integers ℓ 2n+2 ∈ Z N 2n+2 , one readily sees that only the quantum number ℓ 2n+2 = 0 plays a role in the low-energy physics as r → 0. From the perspective of the non-compact dimensions, one gets a tower of very heavy Dirac fields, and focusing on low energies amounts to a truncation of such high-energy modes [78].In the presence of an additional scalar field σ ⋆ (x) that is minimally coupled to the fermions ∂ 2n+2 → ∂ 2n+2 + iσ ⋆ (x), and assuming that this scalar field is homogeneous σ ⋆ (x) = m ⋆ , the dimensional reduction leads to an effective low-energy action action that reads which is now invariant under Lorentz transformations in the reduced spacetime, such that SO(2n + 3) → SO(2n + 2).
By comparing this dimensionally-reduced action to that of Eq. (A1), we see that the sigma fields play a similar role to the mass terms when they are homogeneous.The main difference is that we have more freedom in the definition of parity, and these two mass terms can open a gap in the parity-symmetric case.For odd spacetime dimensions, parity must be understood as a transformation that reverses only an odd number of the spatial directions, for example Pψ(τ, x x x) = γ 2n ψ(τ, (B6) Such a transformation can be defined for any other spatial axis, or an combination of an odd number of them.We note again that this effective action (B4) is SO(2n + 1) invariant in the reduced spacetime, but there is a higher SO(2n + 3) invariance if one considers rotations in the full spacetime with the two additional compactified dimensions.The important point is that, as a remnant of the compactified dimensions, the spinors inherit the dimensionality given by the corresponding representation of the higher-dimensional Clifford algebra (see Appendix A).This enlarged number of spinor components and bigger symmetry group can play an important role once we introduce interactions.Indeed, we can also consider introducing SO(2n + 3) invariant quartic interactions, which can be obtain by considering the transformations of fermionic bilinears built from the Clifford algebra elements (A5).In particular, we can add an SO(2n + 3)-invariant 4-Fermi term to the action S int = d D x g 2 2 −(ψψ) 2 − (ψiγ a ψ)(ψiγ a ψ) . (B7) The first term corresponds to the so-called Gross-Neveu interaction [7], and is a scalar under the SO(2n + 3) Lorentz transformations.The other 2n + 3 terms are quartic interactions corresponding to the so-called Thirring term [81], which is a vector-vector interaction that is also invariant under SO(2n + 3) Lorentz transformations.On the other hand, if one rewrites these terms as (ψiγ a ψ) 2 = (ψiγ µ ψ) 2 − (ψγ 2n+1 ψ) 2 − (ψγ ⋆ ψ) 2 , and reinterprets them from the perspective of the dimensionally-reduced spacetime, only the (ψγ µ ψ) 2 terms correspond to the squared magnitude of a vector under SO(2n + 1), whereas the two additional terms (ψiγ 2n+1 ψ) 2 and (ψiγ ⋆ ψ) 2 are scalars.We emphasise that these additional terms are only allowed by the fact that we are working with a reducible representation of the gamma matrices, which are allowed by the larger number of spinor degrees of freedom in the enlarged spacetime.In the main text, we have referred to the action S 0 + S int in Eqs.(B3) and (B7) for n = 1 as our model of Dirac matter with SO(2n + 3)=SO(5) 4-Fermi interactions.We discuss possible lattice regularizations that allow us to discuss higher-order topological phases and competing symmetry-breaking condensates as one increases the strength of the quartic interactions.This regularization requires rewriting the dimensionally-reduced masses in Eq. (B4) in terms of twisted Wilson masses, as discussed in the main text.

Appendix C: Standard Wilson fermions and Chern insulators
In this Appendix, we present the details for a standard Wilson discretization of the reducible Dirac QFT in Eq. ( 1).As noted in the main text, this regularization amounts to the introduction of a momentum-dependent shift m 1 → m 1 (k k k) of one of the masses in Eq. ( 5), whereas the other mass is set to zero m 2 → m 2 = 0.The Wilson mass depends on a real parameter r as follows which can be understood as the consequence of a finitedifference discretizations of terms involving higher-order spatial derivatives [84].We note that we have used an overline in the function in order to differentiate this standard Wilson mass from the twisted Wilson mass in Eq. ( 9) For this regularised model in D = 2 + 1 dimensions, one actually finds that it corresponds to two copies of the squarelattice version [91][92][93][94] of Haldane's quantum anomalous Hall effect [90], leading to a Chern insulator.The full band structure consists of four energy bands with a two-fold degeneracy ε q,± (k k k) = ±ε(k k k) for q ∈ {1, 2}, where The groundstate is then obtained by filling all the negative energy states |gs⟩ = ∏ k k k∈BZ ∏ q=1,2 ε q,− (k k k) , and one clearly sees from the above dispersion that the Wilson term leads to a different mass for each Dirac point (7), namely As outlined above, by setting m 1 = 0, one sees that the spurious doublers become very heavy with a mass of the order of the lattice cutoff m 1 (k k k D,ℓ ℓ ℓ ) ∝ r/a, ∀ℓ ℓ ℓ ̸ = 0 0 0. On the contrary, the fermion at the origin of the BZ remains massless; m 1 (k k k D,0 0 0 ) = 0. Making a long-wavelength expansion around this point k k k → k k k D,0 0 0 + k k k for |k k k| ≪ Λ c now yields a longwavelength action that coincides with Eq. (1) for m 2 = 0. We remark that, although the lattice discretization breaks explicitly the invariance under SO(3) Lorentz transformations, one recovers it in the continuum limit around k k k D,0 0 0 .Let us finally discuss the connection of the groundstate of this lattice QFT to standard first-order topological insulators, in particular, to the so-called Chern insulators.For the explicit choice of the gamma matrices (2)-( 3), there is a block structure that can be exploited to find a basis in which the problem reduces to a pair of decoupled Chern insulators.Indeed, one can prove that the groundstate corresponding to the above Dirac sea can be characterised by a non-vanishing Chern number for the principal U(1) bundle associated to the filled bands [229].This topological invariant can be expressed as the integral of the Berry curvature F i j b,q (k k k) = ∂ k j A i q (k k k) − ∂ k i A j q (k k k), where the connection is A i q (k k k) = −i ε q,− (k k k) ∂ k i ε q,− (k k k) [230].One can then show that, for our Wilson-fermion QFT, the Chern number is Ch = 1 4π ∑ q dk i ∧ dk j F i j b,q (k k k) = (C6) In light of Eq. (C3), we thus see that, whenever m 1 a ∈ (−4r, −2r) ∪ (−2r, 0), there is a non-vanishing Chern number Ch = ±2, signalling that we have two copies of the standard Chern insulator Ch = ±1, each corresponding to the squarelattice version [91][92][93][94] of Haldane's quantum anomalous Hall effect [90].Even if there is no net external magnetic field piercing the spatial lattice, the system displays a quantised Hall conductance that is related to the Chern number as in the standard quantum Hall effect [231].The bulk-boundary correspondence links these topological invariants to the appearance of circulating edge states localised at the spatial boundaries, which are in fact low-dimensional versions of Kaplan's domain-wall fermions in lattice field theories [110][111][112].
Appendix D: Calculation of the effective SO(5) potential Remarkably, all resummations required for the calculation of the effective potential V eff (Φ Φ Φ) are already encountered in the "vanilla" Gross-Neveu model in 2 + 1 dimensions.We can rewrite this QFT in terms of an Euclidean Lagrangian containing just a single bosonic auxiliary field where φ (x) will be latter identified with the various components of the completing condensates Φ Φ Φ in our problem.In the large-N limit for a condensate Φ = ⟨φ ⟩, V eff then contains the sum of all diagrams with a single fermion loop and n external φ legs [205] (see Fig. 8) which contribute with The constant A k in these expressions is defined such a way that one recovers A nn = 2 n−1 .Calculation of the sum in (D2) proceeds by considering a resummation of three different cases.

4(−1)
where we have made use of the index q = n/2.(b) n even and k even: When the integer k is allowed to be even, we reindex using n = 2q and k = 2ℓ to find (D8)

I. Introduction 1 II. Euclidean 4 -
Fermi field theories with a twist 2 A. Dimensional reduction and SO(5) → SO(3) 2 B. Anisotropic mass twisting of Wilson fermions and the explicit breakdown of SO(3) symmetry 3 III.Cold-atom Hamiltonian field theory 5 A. The Creutz-Hubbard multi-layer model 5 B. Spin-3/2 atoms and 4-Fermi interactions 6 C. Raman lattices and twisted Wilson fermions 8 IV.Higher-order topological insulators (HOTIs) and hidden SO(5) symmetry 9 A. Flat bands and zero-energy corner modes 10 B. Higher-order topological invariants 11 V. Correlated HOTIs, fermion condensates and quantum phase transitions 13 A. Auxiliary fields and the hidden symmetry 13 B. Large-N condensates and the effective potential 15 C. Self-energy and correlated HOTIs 16 D. Large-N results and phase diagram 16 VI.Conclusion and Outlook 18 I. INTRODUCTION

e 1 -flux tunneling along π e 2 FIG. 1 .
FIG. 1. Multi-layer scheme for the lattice regularization of the SO(5) 4-Fermi QFT: (a) In the Hamiltonian formulation, the spatial directions are discretized on a square lattice, and the spinor components can be represented as four-layer model.The yellow region in the inset (b) describes the discretization of the kinetic and twisted Wilson mass along the x axis.The layers only get coupled in pairs, and can be interpreted as a distribution of rhomboid plaquettes pierced by π-fluxes corresponding to a pair of decouple Creutz ladders.The blue region in the inset (b) describes the discretization of the kinetic and twisted Wilson mass along the y axis.This tunnelling couples all four layers, and also has underlying π-fluxes which, together with those of bf (b) lead to flat bands.The purple region in the inset (d) represents the SO(5) 4-Fermi interactions, which include Hubbard density-density interactions among all possible inter-layer pairings (upper panel), and also spin-changing collisions that can be interpreted as inter-layer pair tunnelling processes.

FIG. 2 .
FIG.2.Raman optical lattice: (a) The atom cloud (light red sphere) is subjected to a 3D optical lattice, created using three pairs of counter-propagating laser beams (blue arrows) with mutually orthogonal polarizations (green arrows) and frequencies ω S .The atoms can be then confined into a 2D plane by increasing the potential depth in the z direction.Additionally, travelling Raman beams (orange arrows) with frequencies ω σ σ ′ j and appropriate polarizations (red arrows) generate spin-changing processes, as discussed in the main text.Finally, a magnetic field B ext is applied in the z direction to split the hyperfine atomic energy levels, and a gradient obtained by lattice acceleration creates and energy difference ∆ between nearest-neighbour sites in the y direction.The panels (b), (c) and (d) depict the two-photon Raman transitions giving rise to spin-changing processes, where ε σ is the electronic energy of |σ ⟩, and δ σ σ ′ j the corresponding detunings.

FIG. 3 .
FIG.3.Anomalous corner states in the flat-band regime: For a finite multi-layer with parameters m 1 = m 2 = r/a and r = 1, the Hamiltonian displays two-fold degenerate flat bands ε q,± (k k k) = ±2t and four zero-energy modes(31) that are strictly localised to the corners of the lattice.The colouring of the sites at the corners corresponds to the complex phases of Eq.(31).

FIG. 4 .
FIG. 4. Non-interacting higher-order topological phase diagram:We represent the topological invariant in Eq. (45) as a function of the bare twisted masses m 1 and m 2 .In the green inner square m j a ∈ (−2r, 0), the topological invariant is non-trivial e iπW 1 W 2 = −1, and the groundstate of the twisted-Wilson lattice model corresponds to a higher-order topological insulator (HOTI).The shaded line with m 1 = m 2 represents the regime where the model has a hidden SO(5) symmetry that protects the corner modes.The white region where m j a / ∈ (−2r, 0) for one or both masses corresponds to a trivial band insulator (TBI) with trivial topological invariant e iπW 1 W 2 = +1.

FIG. 6 .
FIG. 6. Interacting higher-order topological phase diagram II:Here the lines where the gap functions f 1 = 0 (blue) and f 2 = 0 (magenta) are plotted for increasing values of the coupling strength g 2 /a ∈ {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7}, which grow as one moves towards the center of the square.The HOTI phase is contained within a star-shaped region that is delimited by groups of four of these lines, corresponding to the same value of the coupling strength g 2 /a.
Appendix B: Dimensional reduction and 4-Fermi interactions

FIG. 8 . 2 ) 2 −
FIG.8.Large-N Effective Potential: Example of a Feynman diagram with n external zero-momentum auxiliary legs (dashed lines), and a single fermion loop (solid circle), yielding the main contribution to the effective potential V eff in the large-N limit.