The rotor Jackiw-Rebbi model: a cold-atom approach to chiral symmetry restoration and quark confinement

Understanding the nature of confinement, as well as its relation with the spontaneous breaking of chiral symmetry, remains one of the long-standing questions in high-energy physics. The difficulty of this task stems from the limitations of current analytical and numerical techniques to address non-perturbative phenomena in non-Abelian gauge theories. In this work, we show how similar phenomena emerge in simpler models, and how these can be further investigated using state-of-the-art cold-atom quantum simulators. More specifically, we introduce the rotor Jackiw-Rebbi model, a (1+1)-dimensional quantum field theory where interactions between Dirac fermions are mediated by quantum rotors. Starting from a mixture of ultracold atoms in an optical lattice, we show how this quantum field theory emerges in the long-wavelength limit. For a wide and experimentally-relevant parameter regime, the Dirac fermions acquire a dynamical mass via the spontaneous breakdown of chiral symmetry. Moreover, we study the effect of both quantum and thermal fluctuations, which lead to the phenomenon of chiral symmetry restoration. Finally, we uncover a confinement-deconfinement quantum phase transition, where meson-like fermions fractionalise into quark-like quasi-particles bound to topological solitons of the rotor field. The proliferation of these solitons at finite chemical potentials again serves to restore the chiral symmetry, yielding a clear analogy with the quark-gluon plasma in quantum chromodynamics, where this symmetry coexists with the deconfined fractional charges. Our results show how the interplay between these phenomena could be analyse in realistic atomic experiments.

Understanding the nature of confinement, as well as its relation with the spontaneous breaking of chiral symmetry, remains one of the long-standing questions in high-energy physics. The difficulty of this task stems from the limitations of current analytical and numerical techniques to address non-perturbative phenomena in non-Abelian gauge theories. In this work, we show how similar phenomena emerge in simpler models, and how these can be further investigated using state-of-the-art cold-atom quantum simulators. More specifically, we introduce the rotor Jackiw-Rebbi model, a (1+1)-dimensional quantum field theory where interactions between Dirac fermions are mediated by quantum rotors. Starting from a mixture of ultracold atoms in an optical lattice, we show how this quantum field theory emerges in the long-wavelength limit. For a wide and experimentallyrelevant parameter regime, the Dirac fermions acquire a dynamical mass via the spontaneous breakdown of chiral symmetry. Moreover, we study the effect of both quantum and thermal fluctuations, which lead to the phenomenon of chiral symmetry restoration. Finally, we uncover a confinement-deconfinement quantum phase transition, where meson-like fermions fractionalise into quark-like quasi-particles bound to topological solitons of the rotor field. The proliferation of these solitons at finite chemical potentials again serves to restore the chiral symmetry, yielding a clear analogy with the quark-gluon plasma in quantum chromodynamics, where this symmetry coexists with the deconfined fractional charges. Our results show how the interplay between these phenomena could be analyse in realistic atomic experiments.
Introduction.-Quantum field theory (QFT) provides a unifying framework to understand many-body systems at widely different scales. At the highest energies reached so far, the standard model of particle physics explains all observed phenomena by means of a relativistic QFT of fermions coupled to scalar and gauge bosons [1]. At lower energies, non-relativistic QFTs of interacting fermionic and bosonic particles form the core of the standard model of condensedmatter physics [2], which explains a wide variety of phases via Landau's seminal contributions of spontaneous symmetrybreaking (SSB) [3] and quasiparticle renormalisation [4]. In the vicinity of certain SSB phase transitions, the quasiparticles governing the long-wavelength phenomena can be completely different from the original non-relativistic constituents [5], and even be described by relativistic models analogous to those of particle physics. In fact, it is the careful understanding of this quasiparticle renormalisation, which yields the very definition of a relativistic QFT [6][7][8], and sets the basis for the non-perturbative approach to lattice gauge theories [9,10].
More recently, the range of applications of relativistic theories has been extended to much lower energies, as they also appear in experiments dealing with the coldest type of quantum matter controlled in a laboratory: ultracold neutral atoms [11][12][13][14][15][16] and trapped atomic ions [17][18][19]. In contrast to condensed matter and high-energy physics, coldatom/trapped-ion experiments deal with quantum many-body systems that can be accurately initialised, controlled, and measured, even at the single-particle level, turning Feynman's idea * daniel.gonzalez@icfo.eu of a quantum simulator (QS) [20,21] into a practical reality [22,23]. One of the unique properties of cold-atom QSs is the possibility of controlling the effective dimensionality of the model in an experiment. This is particularly important in a QFT context, where interactions tend to be more relevant as the dimensionality is lowered [6], bringing in an increased richness in the form of non-perturbative effects. Moreover, the reduced dimensionality sometimes captures the essence of these non-perturbative effects, characteristic of higher-dimensional non-Abelian gauge theories, in a much simpler arena. Some paradigmatic examples of this trend are the axial anomaly of the Schwinger model [24,25], the strong-weak duality of the Thirring model [26,27], asymptotic freedom and dynamical mass generation in the Gross-Neveu model [28], and the fractionalization of charge by solitons in the Jackiw-Rebbi model [29]. Therefore, the first QSs of QFTs are targeting models in low dimensions [30][31][32][33][34][35].
We note that the flexibility of these platforms offers an exceptional alternative: rather than using the QS to target a QFT already studied in the realm of high-energy physics or condensed matter, one can design the QS to realise new QFTs which, although inspired by phenomena first considered in these disciplines, lead to partially-uncharted territory and give an alternative take on long-standing open problems in these fields. For instance, despite the huge success of the standard model of particle physics, the absence of fractionallycharged quarks from the spectrum still presents unsolved questions in quantum chromodynamics (QCD), such as understanding the specific microscopic mechanism for the confinement of quarks into mesons/hadrons with integer electric charges [36]. One related problem that remains open is the nature of the confinement-deconfinement transitions at finite temperatures/densities [37] that lead to phases with isolated quarks and gluons-as well as its relation with the restoration of chiral symmetry [38]. Unfortunately, gauge theories in (1+1) dimensions are all confining regardless of the Abelian or non-Abelian nature of the gauge group. In higher dimensions, deconfinement is usually driven by four-body plaquette interactions, which are challenging to implement in cold-atom experiments [30,39]. Therefore, rather than looking for QSs of gauge theories, one can exploit the aforementioned flexibility of QSs to design new QFTs where characteristic phenomena of higher-dimensional non-Abelian gauge theories emerge in strongly-correlated phases. In this article, we follow this route, and identify a simple lattice model in (1+1) dimensions that regularises a relativistic QFTs where the interplay between dynamical mass generation and charge fractionalization leads to confinement-deconfinement transitions of quarklike quasi-particles, the mechanism of which can be neatly understood at the microscopic level. Although various mechanisms of confinement in QSs have been discussed in the literature [40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55], to the best of our knowledge, our work identifies for the first time a confinement-deconfinement transition between fractionally-charged quasi-particles. Moreover, the feasibility of our QS proposal with state-of-the-art cold-atom experiments indicates the possibility of experimentally observing this transition, together with other QCD-like phenomena, such as chiral symmetry restoration at finite densities.
The lattice model.-To motivate the nature of our model, we note that imposing non-linear constraints in QFTs is also a source of non-perturbative phenomena that resemble the phenomenology of non-Abelian gauge theories. The O(N) nonlinear sigma model, where a vector field is constrained to take values on the (N − 1)-sphere, also displays asymptotic freedom and dynamical mass generation [10,56], although the latter cannot be accompanied by SSB in (1+1) dimensions [57]. Remarkably, the O(3) non-linear sigma model with an additional topological term arises as the long-wavelength description of Heisenberg antiferromagnetic spin chains [58,59]. For leading antiferromagnetic correlations, these systems are effectively described by quantum rotor models, which consist of particles rotating in the surface of a sphere, such that their angular momentum competes with the interactions that favour a collective orientation [58][59][60]. This motivates our study of constrained QFTs involving rotors as mediators of interactions between Dirac fermions.
Let us introduce a (1+1)-dimensional lattice model that can display SSB of a chiral symmetry, leading to dynamical mass generation, which allows, as we will see, the emergence of confinement. We consider a Hamiltonian lattice field theory of fermions c i , c † i and spins S S S i , both residing at the sites x i = ia of a 1D chain of N s sites and length L = N s a (see Fig. 1). The fermions hop between neighbouring sites with tunnelling strength t across a potential landscape set by a spin-fermion coupling g g g = ge z , and determined by the spin along the quantisation axis  [58,59], we show that long-range antiferromagnetic order can take place even in (1+1) dimensions, as it is the result of the breakdown of a discrete chiral symmetry ( Fig. 2(a)). These properties can be neatly understood in the long-wavelength limit, where we find an effective Jackiw-Rebbi-type QFT with rotor fields playing the role of the self-interacting scalar field: a rotor Jackiw-Rebbi model. Interestingly, in contrast to the standard Jackiw-Rebbi model [29], the SSB does not take place at the classical level, but requires genuine quantum effects that lead to the non-perturbative dynamical generation of a fermion mass. We also explore how different rotor profiles can arise in the groundstate by varying the fermion density, which either leads to fractionally-charged or fermionbag quasi-particles, and proof their stability even in the ultimate quantum limit of spin S = 1/2. In fact, we show that the fermion-bag quasi-particles can be understood as confined pairs of fractional charges, resembling the confinement of fractionally-charged quarks in mesons that occurs in the standard model of particle physics ( Fig. 2(b)). Interestingly, we find that a confinement-deconfinement transition can be controlled by tuning a single microscopic parameter, and that this quantum phase transition is associated to the restoration of chiral symmetry by soliton proliferation.
Similarly to lattice gauge theories, in our lattice model (1), fermion-fermion interactions are mediated by bosonic fields. FIG. 2. Phase diagram of the rotor Jackiw-Rebbi model: (a) In the figure, we represent qualitatively the different phases that appear in the half-filled vacuum in terms of the interaction g and the longitudinal field h , for fixed values of the transverse field h t and temperature T . Chiral symmetry is spontaneously broken in the shaded region (χ-SSB), where the fermions develop a dynamical mass and the spins display Néel long-range order, as depicted in the upper panel of (b). This region is surrounded by a chiral-symmetric phase with a longitudinal paramagnet for the spins, such that the interacting massless fermions form a Luttinger liquid (LL). (b) Within the ordered region, we find two different quasi-particle regimes, separated in the figure by a dashed line. In the first one, deconfined topological defects in the spins bound repulsive quark-like fermions with fractional charges. In the second one, the quasi-particles attract each other, forming meson-like fermion bags with integer charge. For a finite doping density, the two regimes are separated by a first-order confinement-deconfinement phase transition, which coincides with a chiral symmetry restoration due to the proliferation of defects.
However, contrary to the former, our model does not possess gauge invariance, a challenging feature to simulate with atomic resources [30]. We show how this simplification allows one to implement the lattice model using state-of-the-art cold-atom QSs in a large regime of realistic experimental parameters. Our results thus show how non-perturbative high energy phenomena, such as charge confinement, could be investigated using minimal experimental resources. Let us start with the proposed scheme for the cold-atom QS.
Cold-atom quantum simulation.-The lattice Hamiltonian (1) can be realised with a Bose-Fermi mixture of ultracold atoms confined in spin-independent optical lattices. We will focus on the S = 1 2 case, as it presents fewer requirements and, moreover, all the relevant phenomena already appear in this limit of maximal quantum fluctuations. We discuss here a generic situation where the bosonic atoms populate a single internal state, whereas the fermionic atoms can be in two possible states, here labelled as σ ∈ {↑, ↓}. In App. A, we present a detailed account for a particular choice of atomic species: 87 Rb atoms in their absolute groundstate for the bosons, and 40 K atoms in two Zeeman sub-levels of the hyperfine groundstate manifold for the fermions. The bosonic and fermionic species will be individually trapped by a single one-dimensional optical lattice, such that the relative strength can be adjusted by tuning the corresponding lattice wavelength. The particular choice of atoms allows for a much deeper lattice for the fermionic species (see App. A), such that the dynamics can be described by a Bose-Fermi Hubbard Here, µ i (µ i,σ ) is the local chemical potential for the bosons (fermions), U (U ↑↓ ) is the on-site Hubbard interaction due to s-wave collisions between a pair of bosons (fermions), and U b↑ ,U b↓ stem from the corresponding fermion-boson scattering. In addition, the Hamiltonian for the internal degrees of freedom of the fermionic atoms reads which includes the atomic energies for the fermions ε ↑ , ε ↑ , and a local term e.g. a radio-frequency ω d that induces oscillations between the two fermionic hyperfine states with the Rabi frequency Ω d (see Fig. 3(b)). We consider an average filling of one fermion per lattice site, effectively reducing the fermionic Hamiltonian in a rotating frame to a sum of on-site terms, whose dynamics can be expressed in terms of spin- 1 2 As discussed in App. A, the strength of the local external driving leads to the transverse field h h h t , and its detuning from resonance The key of the present scheme is the possibility of controlling the different Hubbard interactions in Eq. (2) with a single magnetic field using suitable Feshbach resonances. As detailed in App. A, the difference between the bose-fermi scattering channels readily implements which can be adjusted using a single Feshbach resonance. Finally, in the limit U t, the Hamiltonian coincides with the targeted lattice model (1) by means of a Jordan-Wigner transformation [65]. In this limit, it becomes clear that we use hardcore bosonic atoms to simulate the fermion fields, and spinfull fermionic atoms to simulate the rotor fields. We note, however, that the relevant physics survives away from this hardcore limit, e.g. we explicitly show a spontaneous breaking of chiral symmetry at a finite value of U/t, giving rise to the generation of a dynamical mass for strongly-correlated bosons (see App. A).
In this Appendix, we also present the detailed calculations of all the relevant parameters for the Rb-K mixture, demonstrating a wide range of experimentally tunability. We would like to emphasise that the above ingredients have been realised in current cold-atom experiments, and the QS of the model (1) is thus realistic. We believe that the main experimental challenge will likely be the preparation of extended regions in the lattice with the proper filling for the bosons and fermions.
In the sections below, we give a detailed account of the results summarised above about the equilibrium states and low-energy quasi-particles of the Hamiltonian (1). When discussing these particular phenomena, we will comment on specific probing techniques and the temperatures that would be required to observe them. Let us start by describing these effects by means of an effective long-wavelength theory.
The rotor Jackiw-Rebbi QFT.-We now derive the continuum Hamiltonian QFT that captures the low-energy properties of the lattice model (1) at half filling, this is, with fermionic density ρ = 1 N s ∑ i c † i c i = 1/2, which will be considered as the vacuum state in the following. As customary in various (1+1)-dimensional systems [66], one starts by linearising the dispersion relation to obtain a Tomonaga-Luttinger model [67,68]. In this case, this amounts to a longwavelength approximation around the Fermi points ±k F = ±π/2a, such that c i ≈ e −ik F x i √ aψ + (x) + e +ik F x i √ aψ − (x) is expressed in terms of slowly-varying right-and left-moving fermion fields ψ ± (x). As detailed in Appendix B, this can be recast in terms of the staggered-fermion discretization of lattice gauge theories [69], where one identifies the Fermi velocity v F = 2ta as the effective speed of light c = v F . To proceed with the continuum limit, the spin operators are also expressed in terms of slowly-varying fields S S S i ≈ cos(k N x i )Sn n n(x) + a (x) with k N = π/a, corresponding to the Néel n n n(x) and canting (x) fields [58,59]. Performing a gradient expansion and neglecting rapidly-oscillating terms yields H = dxH (x) with where Ψ(x) = Ψ † (x)γ 0 is the adjoint for the spinor Ψ(x) = (ψ + (x), ψ − (x)) t , the gamma matrices are γ 0 = σ x , γ 1 = −iσ y , and the charge-density is j 0 (x) = Ψ(x)γ 0 Ψ(x). Additionally, we have introduced the coupling g g g s = gSe z , and ∂ 1 = ∂ /∂ x.
This relativistic QFT (6) can be understood as a Jackiw-Rebbi-type model [29], in which a Yukawa term couples the fermion-mass bi-linear Ψ(x)Ψ(x) to the Néel field n n n(x), instead of the standard Yukawa coupling to a scalar field φ (x). Additionally, the rotor dynamics is determined by the precession of the canting field under g g g s j 0 (x)/S − h h h, instead of the more familiar λ φ 4 term of the Jackiw-Rebbi QFT. Hence, this precession includes the back-action of the matter field onto the mediating fields via the charge density. In the large-S limit, and in phases dominated by Néel correlations, n n n(x) and (x) represent, respectively, the orientation and angular momentum of a quantum rotor lying on the unit sphere, such that this QFT (6) can be understood as a rotor Jackiw-Rebbi (rJR) model. In contrast to the standard rotor model [60], neighbouring rotor fields are not coupled via O(N)-symmetric interactions, but rotor-rotor couplings will instead be generated through their coupling to the Dirac fields, and viceversa. As shown below, the lack of a continuous symmetry in Eq. (6) plays a crucial role, and makes the physics of the rJR model very different form these O(N) counterparts.
If the spins order according to a Néel pattern n z (x) = 0, fermions will acquire a mass by the SSB of the discrete chiral symmetry Ψ(x) → γ 5 Ψ(x) with γ 5 = σ z , and n z (x) → −n z (x). However, in contrast to the JR model [29], this chiral SSB cannot be predicted classically by looking at the symmetrybroken sectors of the scalar field, i.e. the double-well minima of the classical potential in the λ φ 4 theory. In our case, the bare rotor term only describes precession of the rotor angular momentum, and does not include a collective coupling of the rotors that would induce Néel order already classically. Instead, in closer similarity to the Gross-Neveu model [28], chiral SSB shall occur by a dynamical mass generation that can only be accounted for by including quantum effects (i.e. rotorfermion loops). On the other hand, there are important differences, as the Gross-Neveu model uses an auxiliary Hubbard-Stratonovich field [70], whereas our rotor fields represent real degrees of freedom with their intrinsic quantum dynamics. As discussed below, this will lead to crucial differences for the chiral SSB, the quasiparticle spectrum, and confinement.
Dynamical mass generation and large-S limit.-Using a coherent-state basis [71,72], as discussed in App. C, the partition function Z = D[Ψ, Ψ, n n n, ]e −S E can be expressed in terms of an action S E = d 2 xL (x x x) with the Lagrangian Here, x x x is a 2-dimensional Euclidean space with imaginary time x 0 = cτ = c(it), and x 1 = x, and the Euclidean gamma matrices areγ 0 = γ 0 ,γ 1 = −iγ 1 . The Dirac and adjoint spinors are Grassmann-valued fields Ψ(x x x) = (ψ + (τ, x), ψ − (τ, x)) t , Ψ(x x x) = (ψ − (τ, x), ψ + (τ, x)), and the charge density is the Grassmann bi-linear j 0 (x x x) = Ψ(τ, x)γ 0 Ψ(τ, x). Likewise, the spins lead to constrained vector fields which, in the large-S and Néel-dominated limits, correspond to the position n n n(x x x) and angular momentum (x x x) of a collection of rotors with n n n(x x x) · n n n(x x x) = 1, n n n(x x x) · (x x x) = 0.
In the continuum limit, Since the Lagrangian (7) is quadratic in Grassmann fields, the fermions can be inte-grated out to obtain an effective action for the rotor fields S eff = − log D[Ψ, Ψ]e −S E . In analogy to the Gross-Neveu model [28], where dynamical mass generation can be derived by assuming a homogeneous auxiliary field, we consider n n n(x x x) = n n n, (x x x) = , ∀x x x ∈ (0, β ] × (0, L]. As shown in App. D, the functional integral leads to S eff = β L g g g 2 − h h h · − (g g g s · n n n) 2 4πc log Λ c g g g s · n n n 2 + 1 , where we have introduced the UV cutoff Λ c = 2t set by the bandwith of the bare fermion dispersion on the lattice (1).
Exploiting the analogy to the Gross-Neveu model [28], where one deals with N flavours of Dirac fermions and finds the non-perturbative dynamical mass generation in the N → ∞ limit, we will send S → ∞ to determine the non-perturbative phase diagram quantitatively. Since phases not governed by Néel correlations can also appear, and we are interested in possible quantum phase transitions transitions thereof, we must relax the first constraint in Eq. (8). By considering the explicit construction of the rotor fields in terms of the spin coherent states, we find the Néel and canting fields The fields are thus parametrized by θ A , θ B ∈ [0, π], each of which represents the angle of the spin coherent state associated to the A (odd sites) and B (even sites) sub-lattice, pointing along the great circle contained in the xz plane (see Fig. 1). One can check that the second constraint in Eq. (8) is readily satisfied n n n · = 0, whereas the first one will only be recovered in Néel-dominated phases θ A ≈ −θ B = ±π/2, where n n n·n n n ≈ 1. By inspecting the effective action (9), one finds that it is proportional to S eff = β LSV eff (θ A , θ B ), such that the large-S limit is obtained through the saddle-point equations ∇ ∇ ∇ θ V eff | θ θ θ = 0 0 0. The features of the phase diagram can be understood in two complementary regimes: (a) For h t /t h /t, g/t, the saddle-point equations will be solved by θ A = θ B = π, such that all spins are maximally polarised in the direction of the leading transverse field |g tP = ⊗ i |S, S x,i , where |S, m α,i is the common eigenstate of S S S 2 i , S α i with eigenvalues S(S + 1), and m ∈ {−S, −S + 1, · · · , S}, respectively. This state can thus be understood as a disordered transverse paramagnet, with all spins aligned along the equator of Fig. 1. Since n z = 0, there is no mass generation, and the chiral symmetry remains intact. Therefore, the fermionic sector will correspond to a metallic Luttinger liquid.
the saddle-point equations are solved by θ A = θ B = π/2, or θ A = θ B = −π/2, respectively. These correspond to disordered longitudinal paramagnets with groundstates |g + P = ⊗ i |S, S z,i , or |g − P = ⊗ i |S, −S z,i , with all spins pointing towards the north or south poles of Fig. 1, respectively. Once again, since g ± P | n z |g ± P = 0, there is no mass generation and the fermions are described by a massless Luttinger liquid. On the contrary, if h − < h < h + , the saddle points correspond to θ A = −θ B = ±π/2, which yields two Néel antiferromagnets |g ± N = ⊗ i∈A |S, ±S z,i ⊗ i∈B |S, ∓S z,i , in which chiral symmetry is spontaneously broken, yielding g ± N | n z |g ± N = ±1 respectively. In both cases, the spins of Fig. 1 alternate between the north and south poles, and the Dirac fermions acquire a mass dynamically, accompanied by a so-called scalar condensate Since the spinor components ψ ± (x) correspond to the longwavelength excitations around momenta ±k F = ±π/2a, the scalar condensate of the continuum QFT leads to a periodic modulation of the lattice density c † i c i = 1 2 + (−1) i Σ 0 a. This phase is reminiscent of a charge-density wave insulator in electron-phonon systems. However, in these systems, Peierls' argument shows that the 1D metal is always unstable towards the insulator [73], instead of showing different phases separated by quantum critical lines, as occurs for our rotor Jackiw-Rebbi QFT (11). Let us remark once more that, while in the standard JR model [29] chiral SSB can be understood by means of perturbation theory about the classically SSB sectors in g 1, the mass is dynamically generated in our case, and cannot be understood perturbatively, which can be appreciated by the fact that the log(St/g) dependence in Eq. (11) cannot be Taylor expanded for small g.
Quantum fluctuations and Ising universality class.-Let us now benchmark the above large-S calculations with numerical results based on a Hartree-Fock (HF) [74] self-consistent mean-field method, and a matrix-product-state (MPS) [75,76] formulation of the density-matrix renormalization group (DMRG) [77]. This serves a two-fold purpose: on the one hand, both methods are directly applied to the discretised model on the lattice (1), and can thus be used to identify the parameter regime where the continuum QFT predictions (9) are recovered. On the other hand, the quasi-exact MPS method gives direct access to corrections of the large-S predictions for finite values of S ∈ { 3 2 , 1, 1 2 }, testing the dynamical mass generation in the regime of large quantum fluctuations. Likewise, we can adapt the HF method to non-zero temperatures T and chemical potentials µ in order to explore the role of thermal fluctuations and finite densities (see Appendix E). Figure 4 contains our results for the zero-temperature halffilling phase diagram as a function of (gS/t, h S/t) for various values of the transverse field h t S/t. In the background, we represent the Néel order parameter, n z = 2 N s ∑ x e z · n n n(x), obtained by averaging over the HF estimate of the Néel field n n n(x) = 1 2S ( S S S 2i − S S S 2i−1 ). In Fig. 4(a), one can see how the HF predicts an intermediate region, here depicted in red, displaying antiferromagnetic long-range order n z ≈ 1 due to the SSB of the discrete chiral symmetry. In order to test the validity of our QFT predictions based on the phenomenon of dynamical mass generation, we benchmark these numerical results against the critical lines of Eq. (11), which are depicted In the presence of finite quantum fluctuations, the HF solution departs from the correct ground state, although it provides a good approximation even for small values of S in the case of small fluctuations. (c) The AF region shrinks as h t S/t increases, and it does it faster for smaller values of S. However, this phase can always be found, even in the limit S = 1/2, provided that gS/t and h S/t are sufficiently large.
as solid black lines in Fig. 4(a). In analogy to the large-N limit of other strongly-coupled QFTs [70], we must rescale the coupling to obtain physical results for S → ∞, such that gS remains finite. From the comparison of the HF and large-S results, we understand that it is the regime of gS t (i.e. couplings much smaller than the UV cutoff), where we can recover the continuum QFT from the lattice discretization, as typically occurs in asymptotically-free lattice field theories.
Note that, in fact, the lattice theory already agrees with the continuum predictions for intermediate couplings gS ∼ 0.25t, which is a sensible fraction of the UV cutoff, and signals the wide validity of the aforementioned scheme of dynamical mass generation. It is worth mentioning, however, that the SSB mechanism that yields the Néel phase is valid for an even wider range of parameters. Indeed, around any of the critical lines of Fig. 4(a), there will be an effective continuum QFT, albeit with different renormalised parameters that require us to rely on numerical methods or experimental QSs. This wider parameter regime is useful in light of the cold-atom realization presented above, which will be able to probe antiferromagnetic correlations in a regime more favourable than the one set by super-exchange mechanisms [78][79][80][81][82][83]. In particular, for larger values of gS, the order survives to larger temperatures, as we will also show below.
The generation of a dynamical mass described in this section can be implemented following the previous cold-atom scheme, and experimentally probed using standard detection techniques. In particular, the scalar condensate (i.e. chargedensity-wave ordering) can be readily probed by measuring the imbalance I = (n A − n B )/(n A + n B ) between the occupation of Rb atoms on A-and B-sublattices by using superlattices [84][85][86] or via noise correlations [16,87,88]. The antiferromagnetic Néel ordering can be further revealed by evaluating the imbalance observable, or the noise correlations, for the fermionic K atoms in a spin-resolved manner, which requires separating the two hyperfine states during the detection by means of a Stern-Gerlach sequence [79], or other similar techniques [80,89]. Since in the SBB phase, however, there are two energetically-degenerate configurations shifted by one lattice site, and the experiment will consist of many independent copies of the one-dimensional chains, the ensemble-averaged observables may fail to signal the phase transition without an additional term that weakly breaks the symmetry between the configurations. If a small symmetry-breaking field cannot be globally introduced in all these copies, one may resort to a combination with quantum gas microscopy [82,[90][91][92][93], which now also enables full spin and charge read-out.
Let us now explore the effect of finite S and non-zero h t . In Fig. 4(a), the circles represent the critical values of the SSB phase transition for different values of S, and are obtained with the MPS method based on the iDMRG scheme for an infinite chain [76]. These critical points are estimated by localising the divergence of the spin susceptibility, χ S = ∂ n z /∂ (gS/t), where we use bond dimension D = 200 and a four-sites repeating unit cell. As can be observed in the figure, for a vanishing transverse field, the critical points for different S are all arranged along the same critical line which, furthermore, delimits the Néel-ordered phase obtained with the HF method, and agrees with the large-S predictions (11) in the regime where we expect to recover the continuum QFT from the lattice regularisation. We note that, in our model, changing the value of the spin S for a vanishing h t = 0, does not modify the quantum fluctuations, such that the large-S prediction works equally well for any value of S. This contrasts the typical situation in models with O(3) symmetry, where the the  5. Scalar condensate and chiral SSB: We represent the scalar fermion condensate Σ 0 = ΨΨ in terms of gS/t for the ground state of a half-filled chain of different lengths N s , with h S = 0.3t and h t S = 0.05t. The results are obtained using DMRG for S = 1/2. Using the critical exponents of the 2D Ising universality class (ν = 1, β = 1/8), the lines cross at the critical point obtained for the infinite system with iDMRG, and collapse for an appropriate rescaling (inset). classical limit is associated with S → ∞, and quantum fluctuations appear as soon as S is finite. Moreover, as outlined in the appendix, we confirm that there is no qualitative distinction in the underlying physics for integer or half-integer spins, as occurs for models with a continuous O(3) rotational symmetry.
In Figs. 4(b) and (c), we represent the phase diagram as the transverse field is switched on, which introduces quantum fluctuations on the spins. In this limit, we observe how the long-range Néel phase shrinks as a result of the competing quantum fluctuations. We also observe that, as the value of S increases, a better agreement with the HF and QFT predictions is obtained, confirming the generic expectation. Note that this agreement is remarkable, given that the considered values of the spins S are still very far away from the large-S limit.
So far, our numerical benchmark has focused on the SSB captured by the Néel field. Let us note, however, that the dynamical mass generation refers to the fermionic sector, and the gap opening is associated to an underlying non-zero scalar condensate Σ 0 = ΨΨ . In order to extract the value of this condensate from the lattice simulations, we use where the expectation value is calculated with the MPS groundstate obtained using a DMRG algorithm with bond dimension D = 200 for finite chains of variable length L = N s a and unit lattice spacing a = 1. In Fig. 5, we represent the finite-size scaling for the scalar condensate of the S = 1/2 model. The crossing of the lines in the main panel serve to locate the critical point of the model gS c /t, which agrees with our previous iDMRG results based on the Néel order parameter. Hence, this shows that the chiral-SSB occurs via the simultaneous onset of a Néel antiferromagnet and a scalar fermion condensate, which corresponds to a charge density wave as seen from the lattice perspective. Moreover, these results allow us to identify the universality class of the corresponding chiral phase transition. As proved by the data collapse shown in the inset of Fig. 5, the critical exponents correspond to those of the (1+1) Ising universality class. Chiral symmetry restoration at finite temperatures: Phase diagram for different T /t and h t S = 0.05t, where we represent the Néel order parameter n z calculated in the HF approximation. The AF phase, where chiral symmetry is spontaneously broken, shrinks as the temperature increases. In the high-energy-physics lore, one says that chiral symmetry is restored at high temperatures.
Thermal fluctuations and chiral symmetry restoration.-Let us now move on to the study of the corrections due to thermal fluctuations. In Fig. 6, we represent sections of the phase diagram as a function of (gS/t, h S/t) for several values of the temperature T obtained with the HF method. One can observe how the area that encloses a dynamically-generated mass, characterised by the Néel order parameter n z , shrinks with increasing T . Therefore, for sufficiently high temperatures, the chiral SSB phase would eventually disappear in favour of the disordered paramagnet and the massless Dirac fermions, both of which respect the chiral symmetry.
As shown in Fig. 6, the required temperature scale for a robust observation of the different phases may lie above that of the tunnel coupling t, which is a very promising feature of our QS scheme. Using state-of-the-art cooling techniques, temperatures as low as T /t = 0.2 have been reported for twocomponent Fermi gases [82] and T /U = 0.05 for bosonic atoms [94]. As displayed in Fig. 6, the dynamical mass generation leading to the anti-ferromagnetic order and the scalar condensate can be accessed at larger temperatures T /t ∼ 0.4. The underlying reason is that the onset of antiferromagnetic correlations does not rest on the super-exchange mechanism, the scale of which is t 2 /U U. In our case, the mechanism is the dynamical breaking of chiral symmetry, the scale /gap of which is directly set by the interaction coupling g ∼ U b f , which can be on the order of the Rb-K interactions U b f .
Let us now interpret these results in light of thermal chiral symmetry restoration in the standard model. In particular, for a large region of its phase diagram, the vacuum of QCD is expected to break spontaneously a chiral symmetry associated to the quarks' flavour [95]. This mechanism is confirmed experimentally by the measured mass of light baryons, such as protons and neutrons, where chiral symmetry breaking yields the largest contribution to their masses [96], while only a small part comes from the masses of their constituent quarks. On the other hand, at very high temperatures or densities, corresponding to the first instants after the Big Bang or to the dense core of neutron stars, respectively, chiral symmetry is restored and quarks become massless, as shown in experiments involving heavy-ion collisions [97]. This phenomenon is captured by several effective theories of nucleons, such as the Nambu-Jona-Lasinio quantum field theory [98] and, as discussed above, also occurs in our model.
There are, however, several unsolved questions in QCD regarding the restoration of chiral symmetry. One is whether a phase transition at finite T c or a crossover exists for intermediate values of the baryon chemical potential µ B [99]. Another one and is the relation to the deconfinement of quarks where, instead of forming hadronic bound states, deconfinement gives rise to a so-called quark-gluon plasma [37]. For large values of µ B , it is not known if both transitions occur simultaneously or, alternatively, intermediate phases exist [100]. Many effective theories, however, fail to address these questions, since they do not include any confinement mechanism even if they correctly capture the essence of dynamical mass generation. As we show in the next section, our model presents such confinement-deconfinement phase transition between fractionally-charged quasi-particles, allowing for the investigation of the interplay between the latter and chiral symmetry restoration in a simple setup.
Confinement of fractionally-charged quasi-particles.-As mentioned before, the fundamental fields of the QCD sector of the standard model correspond to fractionally-charged fermion fields, the so-called quarks, coupled to bosonic Yang-Mills fields, the so-called gluons [1]. In contrast, the fundamental fields of our model (6) are fermion fields with integer charges coupled to the constrained Néel and canting fields. As noted in the introduction, however, the renormalised quasiparticles of a strongly-coupled QFT may sometimes differ completely from its fundamental constituents. As shown below, our QFT (6) displays a Jackiw-Rebbi-like mechanism of fractionalisation [29], whereby soliton configurations of the Néel field host localised fermions with a fractional charge q = ±e/2, which will play the role of quarks, allowing us to discuss various aspects of a confinement mechanism.
In Figs. 7(a)-(b), we present the MPS numerical results for the real-space configurations of the Néel lattice field n z ( j) = 1 2S ( S z 2 j − S z 2 j−1 ), and the integrated fermion charge above the half-filled vacuum N( j) = Q( j)/e = ∑ i< j ( c † i c i − 1 2 ) when one extra fermion is introduced above half filling. These figures show that the Néel field presents a kink-antikink pair that interpolates between the different SSB sectors, and that each of these topological solitons hosts a localised fermionic excitation with charge q = e/2, henceforth referred to as a 'quark' by the analogy with the fractionally-charged fundamental fermion fields of QCD. We note that similar quasiparticles with charges q = −e/2 would appear for dopings below half filling, playing the role of 'anti-quarks', and that quark-antiquark pairs could appear in the vacuum due to thermal fluctuations, as they correspond to higher-energy states of the theory.
Contrary to the general situation in (1+1) lattice gauge theories, which can only host confining phases [36], we can find regimes where quarks/anti-quarks can be confined/deconfined depending on the microscopic parameters. In order to explore this phenomenon, let us remark that the distance of the pair  of quarks of Fig. 7(b) is determined by the external pinning of the Néel solitons of Fig. 7(a). As described in more detail in Appendix F, we introduce an external potential that breaks explicitly the translational invariance and localises the solitons, which would otherwise travel freely through the chain, at the desired positions. This pinning makes our fractionallycharged quasi-particles static, such that we can discuss the analogue of the static quark potential [101]. We note that it is a standard practice in lattice gauge theories to use this terminology whether or not the charges actually correspond to dynamical quarks [36]. Therefore, the static quark potential quantifies the interaction energy between any pair of static charges as their distance is modified, giving information about confinement not only in QCD, but in other effective models.
In our model, the static quark potential V 1 (d) = E 1 (d) − E 0 can be obtained using the MPS numerics by calculating the energy E 1 (d) of the doped system with an extra fermion that fractionalises into a pair of quarks pinned at a distance d, measured in unit cells (2a), with respect to the energy E 0 of a pair pinned far apart (d 1). We note that the standard situation of (1+1) lattice gauge theories, such as the Schwinger model [24], is that the preservation of gauge symmetry requires that the static charges must be connected through an in- The potential energy between them decreases with the distance, and the minimum is reached when they are located at neighbouring sites.
termediate electric-field string, such that the energy increases with the separation d and leads to a linearly-increasing static quark potential [102,103]. In our case, the situation is completely different, as the energy of the deformed Néel field that connects the two quarks in Fig. 7(c) is independent of the pinning distance, and confinement is thus not enforced a priori. As argued below, there exists a competing mechanism that either favours confinement or deconfinement depending on the microscopic parameters of the model.
In Fig. 7(c), we depict the distance-dependence of the static quark potential for the S = 1/2 rotor Jackiw-Rebbi model for coupling gS/t = 0.9, and setting the other parameters such that we lie in the chiral-SSB phase. As can be observed in this figure, the potential decreases with the distance for small soliton separations, which means that the quarks repel each other. Hence, in the absence of the external pinning, the fractionally-charged quasiparticles would move freely at large distances from each other, and thus appear as asymptotic excitations in the spectrum of the rJR QFT. As we increase the coupling to gS/t = 1.1, Figs. 7(d)-(e) show that a completely-different quasi-particle emerges. In this case, the Néel field no longer interpolates between the two SSB sectors, but is instead suppressed within a small region of space where, as shown in Fig. 7(e), an integer-charged fermion is localised. This situation is reminiscent of the so-called quark bag models [104,105], where quarks and gluons are locked within a finite region of space, in which a phenomenological term that compresses the bag compensates the outward pressure of the quarks that are held inside, and confinement results from the competition of these two terms. The present situation is closer in spirit to the soliton quark model [106,107], where the fermions deplete a SSB condensate in a finite region of space, gaining kinetic energy at the expense of the cost of deforming the condensate. In our case, as the Néel field vanishes in the inner region, there scalar condensate will become zero Σ 0 (x) = 0, ∀x ∈ [x 0 −ξ , x 0 +ξ ], such that the bound fermions have a vanishing dynamically-generated mass, increasing their kinetic energy and the outward pressure. As mentioned above, this is compensated by the energy cost due to the inhomogeneous layout of the Néel field and the accompanying scalar condensate. Longrange order can be detected using the structure factor S(k), which shows a peak at a momentum commensurate with the fermionic density, k D /π = 2N f /N s . (c) In the confined phase, the attractive quark bags create an extensive region where the dynamical fermionic mass is screened to zero. (d) In this case, the peak at k D disappears, while S(k) is non-zero around k C = π, signalling a reduced but nonvanishing AF order. The insets show the finite-size scaling of the peak in S(k) for each case, where the density is kept fixed.
To substantiate this neat picture and connect it to the quark confinement, we should provide evidence that the integervalued charges shown in Fig. 7(e) are the result of an attractive force between the fractionally-charged quasiparticles of Figs. 7(a)-(b). This evidence of confinement is supported by the numerical results presented in Fig. 7(f), where we show that the static quark potential increases with the inter-quark distance in this case. Hence, as advanced in the introduction, it is possible to understand the microscopic confinement mechanism in our model. In the regime g < g c , the outward pressure of the quarks overcomes the inner force that tends to re-establish the homogeneity of the condensate, such that the quarks get deconfined and can move synchronous with the kin/antikink. This changes for g > g c , where the inward force to reestablish condensate homogeneity prevails, and the quarks get confined within the so-called quark bag. For h S/t = 0.4 and h t S/t = 0.01, we estimate a critical value of this confinement-deconfinement phase transition to be g c S/t ≈ 1.01 (see App. F).
We note that these integer-charged excitations also occur in other QFT with a non-classical scalar condensate due to chiral-SSB, such as the Gross-Neveu model [104]. However, to the best of our knowledge, there is no deconfinement transition where they become a pair of fractionally-charged fermions and, additionally, they require at least two fermion flavours to exist, which anyway masks the occurrence of fractionalization as happens for polyacetilene [108,109].
Quark crystals and chiral symmetry restoration.-Let us now move towards finite densities, and explore the rJR groundstate properties in the confined and deconfined regimes. In the deconfined regime, we have shown that the isolated quarks repel each other, such that they will maximise the inter-quark distance and broaden the density profiles if left unpinned. Accordingly, for finite fermion densities above the half-filled vacuum, one would expect the formation of a crystalline structure of equidistant kinks and anti-kinks with the corresponding periodic distribution of fractionallycharged fermions, namely a quark crystal.
On the other hand, in the confined regime, we have no a priori intuition of the possible groundstate ordering. To gain such intuition, let us first calculate the static potential between two distant quark-bag excitations, each of which contains a pair of confined quarks and thus an integer-charged fermion. The static bag potential can be obtained using the MPS numerics by considering in this case the energy of the system doped with a pair of fermions that are held inside the bags E 2 (d), and pinned at a distance d, where E 0 is again the energy of two quark bags pinned far apart. In order to fix the quark-bag distance as depicted in Fig. 8(a), we impose again an external potential that localises them at two given locations (see App. F). Our numerical results show that the static fermion-bag potential of Fig. 8(b) increases with the inter-bag distance, proving that the quarkbags attract each other. Therefore, if left unpinned, we expect that the two bags will merge yielding a wider depletion region of the condensate that can accommodate two integer-charged fermions, each of which can be thought of being composed of two confined quarks. This trend can be generalised to finite density regimes, where we expect the appearance of an extensive quark bag that is sufficiently wide to host all of the extra fermions. In the context of ultracold atoms, this phase can be understood as a phase separation phenomenon.
Let us now confirm this intuition by presenting the MPS numerics for the finite-density regime. In Fig. 9(a) and (c), we display the real-space dependence of the scalar condensate for the deconfined and confined regimes without the pinning potentials. As can be clearly observed, Fig. 9(a) presents a periodic sequence of kinks and antikinks, each of which hosts a single localised quark, giving rise to the aforementioned quark crystal. On the other hand, as we increase the coupling strength, Fig. 9(c) displays an extensive region where the Néel field vanishes, and the dynamical fermionic mass is screened to zero. This wide bag accommodates for all the extra fermions, leading to a phase separation with respect to the region where the vacuum displays a large dynamicallygenerated mass inhibiting the penetration of the massless confined quarks.
We note that the corresponding phases can be quantitatively distinguished by means of the static spin structure factor S(k) = 1 N 2 s ∑ i, j S z i S z j e ik(i− j) , which will peak at different momenta k D , k C for the deconfined/confined phases. For the deconfined quark crystal, Fig. 9(b) shows that the corresponding peak of the structure factor occurs for a momentum that is commensurate with the fermionic density modulation of the scalar condensate k D = 2πN f /N s . Conversely, for the confined phase-separated bag phase, Fig. 9(d) shows that the peak at k D vanishes, and one gets instead a non-zero structure factor FIG. 10. Deconfinement and chiral symmetry restoration: Phase diagram for different values of T /t obtained in the HF approximation for a periodic chain with N s = 64 sites and N f = 35 fermions, with h t S = 0.01t. The order parameter O (see main text) allows us to distinguish between the three different phases that appear for a finite doping. As the temperature increases, the phase transition line between the quark-bag and the quark-crystal phases gets modified, until both of them disappear at a sufficiently high value of T /t. around k C = π, signalling Néel order, which is partially broadened by the condensate deformation due to the quark bag. The inset of both figures displays the finite-size scaling of each peak, where we increase both the size N s and the number of fermions N f , such that the density ρ = N f /N s remains fixed. The extrapolated non-zero values of the corresponding peaks for 1/L → 0 show that the quark-crystal and phase-separated bag phases are both stable in the thermodynamic limit.
The vanishing value of of the structure factor at momentum k = π in the quark crystal suggests the possibility of a zero-temperature quantum chiral symmetry restoration for finite dopings-since this quantity is equivalent to the Néel order parameter n z used in our discussion of thermal chiral symmetry restoration for the rJR vacuum. This is confirmed by calculating the average value of the fermionic condensate where we get Σ 0 = 0.06 and Σ 0 = 0.48 for the parameters used in Fig. 9(a) and Fig. 9(b), respectively. In this case, it is not the thermal fluctuations, but instead the finite density of topological solitons, which reduces the average value of the fermionic condensate to zero (up to finitesize effects), indicating that chiral symmetry is restored in the quark-crystal phase. In this phase, therefore, chiral symmetry coexists with deconfinement. The situation is analogous in QCD, where both properties appear in the quark-gluon plasma. The presence of a single transition from this phase to a confined symmetry-broken phase, or the possibility of intermediate phases with only one of these properties, is still an open question in particles physics [100]. The investigation of such interplay in simple models could help to gain a better understanding of it in more complicated theories, specially in regimes where the chemical potential is large and Monte Carlo simulations suffer from the sign problem [99].
Let us thus explore this interplay in the presence of thermal fluctuations. Figure 10 depicts the HF phase diagram for a finite density of fermions over the half-filled vacuum for different values of T /t. To distinguish the two phases, we use as an order parameter the difference between the structure factor at the two different momenta of the confined and deconfined phases, O = S(k C ) − S(k D ). This quantity is zero in the disordered paramagnetic phase, which corresponds to a Luttinger liquid, as in the case of half filling. Positive and negative finite values of O corresponds, on the other hand, to the quark-bag and quark-crystal phases, respectively. For T /t = 0, we can clearly distinguished these three phases. Note that, it is only in the ordered phases, surrounded by a disorder LL, where the notion of confinement and deconfinement of fractionallycharged quasi-particles is well defined. Within this region, we observe that both the quark-bag and the quark-crystal phases have a finite extension, with a phase transition line separating them. The ordered region shrinks as the temperature increases (Fig. 10). It is interesting to notice that, in our model, the quark crystal disappears more rapidly than the quark-bag phase.
Confinement-deconfinement phase transition.-As we have shown in the previous sections, a characteristic feature of our (1+1)-dimensional QFT (6) is the possibility to understand the mechanism of confinement microscopically and, moreover, the existence of a deconfinement quantum phase transition as we vary the microscopic parameters. In order to study the latter in more detail, we make use again of the static structure factor peaks, which serve as order parameters to detect the corresponding phase transitions. In this section, we study two different types of phase transitions occurring at finite densities with DMRG, confirming the results obtained above using the HF method. The first one, corresponding to Fig. 11(a), describes the transition between the quark crystal and a longitudinal paramagnetic phase. As shown in this figure, following the the spin structure factor at the two characteristic momenta k D and k C , one can see that we start from a disordered phase at small interactions, where both S(k D ) and S(k C ) are zero. As gS/t increases, the quark crystal order parameter reaches a non-zero value S(k D ) > 0, while S(k C ) remains zero. This phase transition is continuous, similarly to the direct order-disorder transition we found for half filling. For larger values of gS/t, we find another phase transition, again in correspondence with the HF results. Fig. 11(b) shows the transition between the deconfined quark crystal S(k D ) > 0, S(k C ) = 0 towards the confined quark-bag phase S(k D ) = 0, S(k C ) > 0. In this case, this deconfinement transition is a first-order phase transition. This is believed to be the case also for the confinement-deconfinement transition in QCD at finite chemical potential. For h S/t = 0.4 and h t S/t = 0.01, this transition is located at g c S/t = 0.94, which roughly agrees with the prediction we obtained using the static quark potential (i.e. g c S/t = 1.01). The difference shows the presence of many-body effects in the case of a quark crystal, where the interaction between two quarks is influenced by the presence of a finite density of them. This agreement supports our claim that the mechanism behind the transition between a quark crystal and a quark-bag phase is the confinement of quark-like fractional quasi-particles.
Conclusions and outlook.-In this work, we studied several high-energy non-perturbative phenomena using a neat (1+1) quantum field theory, the rotor Jackiw-Rebbi model, and proposed a quantum simulation scheme using a Fermi- Bose mixture of ultracold atoms in an optical lattice. Dirac fermions, whose interactions are mediated by spin-S rotors in this model, acquire a dynamical mass through the spontaneous breaking of chiral symmetry. The generation of a mass term in the theory is accompanied here by antiferromagnetic order in the rotors, that we predict analytically in the large-S limit of the continuum model. Using a lattice model that regularises the theory, we studied the phase diagram at half filling in the presence of quantum and thermal fluctuation, showing how dynamical mass generation, in particular, survives in the ultimate quantum limit of S = 1/2. We also showed how, in this limit, the chiral symmetry breaking quantum phase transition lies in the Ising universality class and, moreover, we observed chiral symmetry restoration at sufficiently high temperatures.
We then focused on the regime of finite chemical potentials, where we find a confinement-deconfinement phase transition between quark bags and a crystal of fractional quarklike quasi-particles. This transition is characterised using the spin structure factor, and its microscopic origin is uncovered by means of the static quark potential. We also showed how deconfinement coexists with a restoration of chiral symmetry, even in the absence of thermal fluctuations. In this case, the latter occurs due to a proliferation of topological solitons at finite densities, drawing an interesting analogy to the quark-gluon plasma of QCD. Our results indicate how confinement between fractional charges could be further investigated in atomic experiments, shedding light into one of the long-standing questions of particle physics.
In the future, it would be interesting to extend the large-S techniques to a different family of lattice models that allow for symmetry-protected topological phases, and which have been studied so far for bosonic matter coupled to link spins in the S = 1/2 limit [110][111][112][113], or for fermionic matter coupled to auxiliary scalar [114][115][116] or gauge [53,117,118] fields, which can be substituted by the Néel and canting fields hereby introduced. These extensions would allow to explore the interplay of confinement-deconiment transitions and topological phases of matter in a completely new framework. It would also be interesting to use matrix-product-state sim-ulations to study real-time dynamics, serving as alternative benchmark for quantum simulations in addition to the finitedensity regime hereby studied, where Monte Carlo simulations for a single flavour of fermion fields are expected to suffer from the sign problem. The model can also be easily extended to higher dimensions, where the simulation proposal can be generalized in a straightforward manner by using higher-dimensional optical lattices, which would also reach the limits of efficient tensor-network numerical techniques. It is precisely in the cases where numerical simulations show limitations where cold atoms represent an efficient alternative to provide a full solution of the quantum many-body problem.

Appendix A: Implementation with cold atoms
We propose to realize the S = 1 2 limit of the lattice Hamiltonian (1) with a mixture of bosonic and spinfull fermionic atoms. Notice that the roles of bosonic and fermion degrees of freedom are interchanged here. In particular, we propose to use hard-core bosons to simulate the dynamical fermionic matter, while the fermionic ones will be used to implement the rotor fields. As will become clear below, this approach is motivated by the specific choice of bosonic and fermionic species, which presents a well-characterised Feshbach resonance that will allow an accurate experimental control of the inter-species scattering, the crucial ingredient in our scheme.
Bose-Fermi Hamiltonian.-We aim at realizing the following grand-canonical Hamiltonian for the Bose-Fermi mixture where H b describes the dynamics of the bosonic atoms of mass m b , H f the one of the fermionic atoms of mass m f , and H b f the interaction between the two species. The bosonic part of the Hamiltonian is defined as follows where b i and b † i are the bosonic creation and annihilation operators acting on lattice site i. This grand-canonical Hamiltonian describes the tunnelling of bosonic atoms in a 1D lattice with strength t, and the Hubbard interactions with energy U. Here, µ i = µ − V b,i is expressed in terms of the chemical potential µ and the on-site optical trapping potential V b,i , and controls the bosonic filling in the local density approximation. The fermionic contribution can be divided in two parts where the first term describes the external motional degrees of freedom Here, f iσ and f † iσ are the fermionic creation/annihilation operators acting on lattice site i with internal state σ = {↑, ↓}. Fermions are trapped in a very deep optical lattice, such that t f t, and we can neglect their tunnelling along the 1D lattice during the timescale of interest. They interact with Hubbard interaction U ↑↓ , and their filling is controlled by the local chemical potentials µ iσ = µ σ −V f ,i , where µ σ is the chemical potential for the fermionic atoms in each internal state, and V f ,i is an optical trapping potential. In addition, the internal degrees of freedom shall be described by where we have introduced the atomic energy levels ε ↑ , ε ↓ for the two states of the fermionic species, and a local driving of frequency ω d that induces local oscillations between these two states with a Rabi frequency Ω d , whereh = 1 henceforth. As will become clear below, this driving stems from radiofrequency radiation, which has a negligible momentum, and one can thus neglect recoil effects that would couple the internal and external degrees of freedom. Finally, the interaction between the two species iŝ which describes the on-site interaction between bosonic and fermionic atoms, and depends on the internal spin state of the fermions, i.e. in general U b↑ = U b↓ . As will be clear for the particular Bose-Fermi mixture discussed below, there are also boson-fermion scattering processes where the internal states of the fermions is changed by populating other bosonic states such that the total angular momentum along the quantisation axis is conserved. Nonetheless, these so-called spin-flipping collisions of strength U sf are negligible for a sufficiently-large difference of the on-site energies U sf |ε ↑ − ε ↓ |.
Let us now discuss the steps to arrive at the desired model (1), and find the relation between the model microscopic couplings and the cold-atom experimental parameters. Using two internal states |↑ and |↓ with different magnetic moments, the energy splitting corresponds to the Zeeman energy and can be tuned via an external magnetic. The driving frequency will be near-resonant and fulfils ω d = (ε ↑ − ε ↓ ) − ∆ d , where ∆ d is the so-called detuning. Let us note at this point that the optical trapping potential V f ,i of Eq. (A2) is assumed to be state-independent, which is generally the case, such that it does not modify the above resonance condition, and need not be included in the present discussion. Moving to an interaction picture with respect to the fermionic operators become f i,σ → f i,σ (t) = f i,σ e −iε σ t . Furthermore, assuming that |Ω d |, |∆ d | ω d , we can perform a rotating-wave approximation such that Finally, by moving to a frame that rotates with the drive frequency, this Hamiltonian can be written as the following timeindependent term where we define the Bloch sphere in a way that Ω d ∈ R, i.e. we use the phase of the driving as a reference for subsequent measurements. Here, we have introduced the following relation between the fermionic operators in the original Schrodinger picture f i,σ , and thosef i,σ in the rotating framê We can now define the following spin-1 2 operators in terms of these fermionic annihilation and creation operators Therefore, to measure the spin operators discussed in the main text, one would have to lock the phase evolution to the one set by the source that drives the transition.
After these derivations, we should enforce that only 1 fermion resides at each lattice site, which can be adjusted by the filling, and maintained by working in the regime where t f t, thus suppressing double occupancies. In this case, we obtain the desired S = S z i = 1 2 limit. Moreover, since the fermion tunneling t f is negligible, super-exchange processes stemming from virtual double occupancies occurring at order t 2 f /U ↑↓ are also negligible, and we can finally arrive at an effective description according to the following grand-canonical Hamiltonian where the boson operators are not altered by moving to the rotating frameb i = b i . In the Hamiltonian above, the coupling constant is defined as and the external field given by the driving term the specifics of which are discussed below. The local chemical potential µ i will be adjusted to vary the filling of the bosonic atoms, which is homogeneous in the central region of the trap. The variation of the filling will allow to explore the different quasi-particle regimes discussed in the main text. from this singular point, i.e. for finite values of U/t, provided that the Hubbard interactions are sufficiently strong. This is true in particular for the AF phase found in the fermionic case, which survives away from the hardcore limit, giving rise to the generation of a dynamical mass for strongly-correlated bosons. In Fig. 12, we represent the finite-size scaling for the Néel order parameter, and find that the disorder-order transitions occurs for interactions of the order U c ≈ 10t. Below this value the soft-core bosons are in a superfluid state, while the spins form a longitudinal paramagnet, similar to the disorder phase in the hardcore limit. Moreover, we find the the phase transition is also in the Ising universality class. Once all the steps for the derivation of the target model (A10) have been discussed for a generic cold-atom setting, let us estimate the specific parameters for a mixture of bosonic 87 Rb and spinfull fermionic 40 K atoms, and discuss the viability of the experimental realisation. The following main ingredients are relevant for the implementation: Optical lattice potential.-The dominant transitions of the two alkali atoms, bosonic 87 Rb and fermionic 40 K are λ Rb,D2 780 nm, λ Rb,D1 795 nm (A13) λ K,D2 767 nm, λ K,D1 770 nm which allows for widely tunable polarizabilities and opticallattice potentials. In order to achieve a large separation of timescales, it is desirable to achieve a large ratio of the tunnelling couplings t/t f , where t is the strength of the tunnel coupling for the bosonic species, and t f denotes the spinindependent tunnelling amplitude of the fermionic species. In essence, this means that the lattice potential experienced by the bosonic species should be much weaker compared to the one seen by the fermionic atoms. In order to reduce offresonant photon scattering, which would result in additional heating, at the same time the detuning from any internal transition has to be maximized. Due to the large fine-structure splitting of Rb, there is a convenient tuning range around the zero crossing of the polarizability at λ bz = 790 nm. This range offers both a wide tunability of the tunnelling ratio t/t f (Fig. 13) and a large detuning from all resonances to minimize heating. Tunable interspecies interactions.-In order to provide good control over the parameter g (A11) that appears in the effective Hamiltonian (A10), we propose to make use of the well-calibrated and easily-accessible interspecies Feshbach resonance between the absolute ground states |F = 1, m F = 1 of 87 Rb and |↑ ≡ |F = 9/2, m F = −9/2 of 40 K [119]. Although there does not seem to be a Feshbach resonance nearby to tune the scattering length between the ground state of Rubidium and |↓ ≡ |F = 9/2, m F = −7/2 of 40 K, the interaction parameter g can still be fully tuned over a wide range of values as shown in Fig. 14 where µ bσ = m Rb m K /(m Rb + m K ) is the reduced mass and m Rb = 86.9 u is the mass of one 87 Rb atom and m K = 39.96 u the mass of one 40 K atom; u is the atomic mass unit. The functions w b (r) and w f (r) denote the Wannier functions of Rb and K respectively, which are different because the two species see a lattice potential of different depth. Through these expressions, we can obtain the parameter g, as discussed below.
Interaction ratio gS/U.-Rubidium in its absolute ground state has a scattering length of a b 100.4 a 0 , and the on-site Hubbard interaction is defined as Coupling between spin states.-Ket us now discuss the specifics of the driving term discussed previously. The external fields h , h t can be realized with a radio-frequency or two-photon microwave transitions at frequency ω d almost resonant with the Zeeman energy difference ∆E Z ↑,↓ = ε i,↑ − ε i,↓ between |↑ and |↓ atoms. For the Feshbach resonance shown in Fig. 14, the resonance occurs at B 0 , which corresponds to ∆E Z ↑,↓ /h ≈ 80 MHz, where h denotes Planck's constant. The energy offset h is then realized by detuning the coupling frequency from resonance, i.e. h = ∆E Z ↑,↓ /h − ω d . With single-photon transitions, Rabi frequencies Ω d /2π of several 10 kHz can be easily achieved, which corresponds to the regime |h t | = |Ω d | t. Moreover, the pair of states |↑ and |↓ is well isolated from the other internal states, even in the presence of this coupling. The energetically closest transition is between the levels |↓ and |F = 9/2, m F = −5/2 , which is detuned by ∼ 5 MHz Ω d /h.
For some of the proposed phenomena, it is desirable to tune the coupling in order to enter the regime where h t is on the order of the tunnel coupling t. The level of control one can achieve is limited by the stability of the magnetic field that defines the Zeeman shift ∆E Z ↑,↓ . We have calculated the detuning from resonance δ h that occurs due to an imperfect control of the external magnetic field (Fig. 15). It has been demonstrated that the fluctuations in the external field can be suppressed below 3 × 10 −6 , as reported in Ref. [120]. The most sensitive regime occurs for h = 0, where the lower bound for h t would be on the order of a few 100 Hz, which coincides with typical experimental values for t. Note, that the typical timescale for these fluctuations is large compared to the duration of the experiment, hence, the detuning δ h can be assumed constant but will fluctuate between individual experimental realizations. one obtains the algebra of position n n n and angular momentum = n n n × (−i∇ ∇ ∇ n n n ) of a quantum-mechanical particle [58,59]. Moreover, in this limit, one also finds that n n n(x) · n n n(x) = 1 + 1/S (B7) such that the particle will be confined to a unit sphere in the large-spin limit S 1. Therefore, in the large-S limit, the Néel and canting fields represent a the orientation of a quantum rotor and its angular momentum, respectively. Combining the expressions for the Dirac (B1) and rotor (B5) fields, and neglecting again rapidly-oscillating terms, we find that the spin-fermion coupling can be expressed as H sf = dxΨ(x)g g g s · n n n(x) + (x)γ 0 /S Ψ(x), g g g s = gSe z .
(B8) The first term can be understood as a Yukawa-like term that couples the fermion bilinear Ψ(x)Ψ(x) to the rotor projection n z (x) instead of the standard scalar field in a Yukawa theory [1]. The second term couples the time-like component of the fermion current j 0 (x) = Ψ(x)γ 0 Ψ(x), i.e. the charge density, to the projection of the rotor angular angular z (x).
Finally, the continuum limit of the spin precession yields In analogy to the original situation (1), the rotor angular momentum is subjected to a magnetic field with both longitudinal and transverse components, being the latter responsible for introducing quantum fluctuations since [ x (x), n z (x)] = 0. Altogether, the continuum limit of the lattice model (1) corresponds to the quantum field theory of Eq. (6).

Appendix C: Path integral formulation of the continuum QFT
In this Appendix, we provide a path-integral derivation of the continuum-limit rotor-fermion QFT (6). This derivation serves for two purposes: one the one hand, it allows to clarify the absence of a topological θ term for any spin S, differing markedly from O(N) non-linear sigma models arising from Heisenberg models [58]; on the other hand, it sets the stage for a large-S limit approach to dynamical mass generation.
We are interested in the partition function Z = Tr{e −β H }, where H is the spin-fermion lattice Hamiltonian (1). In the β → ∞ limit of zero temperature, this partition function contains all the relevant information about quantum phase transitions related to the chiral SSB and dynamical mass generation we are seeking. Using fermionic and spin coherent states in Euclidean time τ = it ∈ (0, β ) [72], this partition function can be expressed as a functional integral in terms of anti-commuting Grassmann fields for the fermions ψ(τ, x i ), ψ (τ, x i ) and commuting vector fields for the spins Ω Ω Ω(τ, x i ) lying in a 2-sphere of radius S, which can be rescaled in terms of unit vector fields Ω Ω Ω(τ, x i ) = Sω ω ω i (τ). Using the resolution of the identity in terms of both types of fields [72], one finds that where the Euclidean action S E = S sf − iSA WZ can be expressed as a functional over the fermion and spin fields where H(ψ , ψ, Sω ω ω) is obtained by substituting the fermion and spin operators in the normal-ordered Hamiltonian (1) by the Grassmann and vector fields. One also finds the so-called Wess-Zumino term, which corresponds to the area enclosed by the trajectory of the spin field in the unit 2-sphere τ)).
In O(3)-symmetric situations, such as those arising in antiferromagnetic Heisenberg models, this Wess-Zumino term plays a crucial role as it is responsible for the mapping to a non-linear sigma model with an additional topological θ term that depends on the integer or half-integer nature of the spins [58]. In the present case, however, there is no rotational symmetry in the classical Hamiltonian H(ψ , ψ, Sω ω ω), and one can see from the particular form of the spin-fermion coupling and the external field that it will suffice to use coherent states pointing along the meridian at vanishing longitude which correspond to the great circle in the xz plane of Fig. 1. Accordingly, ω ω ω i · (∂ s ω ω ω i × ∂ τ ω ω ω i ) = 0 and there is no enclosed area by the precession of the spins, such that the Wess-Zumino term vanishes A WZ = 0. We can now introduce the equivalent of the slowly-varying quantum fields in Eqs. (B1) and (B5) in terms of the Grasmmann fermion fields and the spin vector fields Proceeding in analogy to Appendix B, we perform the gradient expansion in the continuum limit a → 0, neglect rapidlyoscillating terms, and find S E ≈ d 2 x L (x x x) with the following Lagrangian density with coupling g g g s and external field h h h defined in Eqs. (B8) and (B9), respectively. Here, x x x = (cτ, x) is a 2-dimensional Euclidean space with metric g E µ,ν = diag(1, 1), and the Euclidean gamma matrices areγ 0 = γ 0 ,γ 1 = −iγ 1 . The Dirac spinor is composed of the right-and left-moving continuum Grassmann fields Ψ(x x x) = (ψ + (τ, x), ψ − (τ, x)) t , such that the adjoint becomes Ψ(x x x) = Ψ † (x x x)γ 0 = (ψ − (τ, x), ψ + (τ, x)). Likewise, the Néel and canting fields are the vector fields n n n(τ, x) = 1 2 (ω ω ω 2i (τ) − ω ω ω 2i−1 (τ)) , (ω ω ω 2i (τ) + ω ω ω 2i−1 (τ)) .

(C8)
In the large-S limit, and in situations dominated by Néel correlations | (a/S) 2 | → 0, these fields are additionally subjected to the rotor constraints n n n(x x x) · n n n(x x x) = 1, n n n(x x x) · (x x x) = 0, which can be included in the partition function through the functional integral measure D[Ψ, Ψ, n n n, ] = (2S+1) Appendix D: Effective rotor action and large-S limit In this Appendix, we give a detailed derivation of the effective rotor action (9). As outlined in the main text, one can integrate out the Grassmann fields from the partition function since the corresponding functional integral reduces to a product of Gaussian integrals. These integrals are obtained after transforming the fields in terms of Matsubara frequencies ω n = π β (2n + 1), and quasi-momentum qa ∈ [−π, π), according to The corresponding integrals lead to det (−iω n + h D (q, g g g s · n n n)) , (D3) where we have introduced a renormalised rotor Lagrangian for a homogeneous canting field (x x x) = . This term (D4), which contains the original precession under the external field −h h h, gets a contribution from the back-action of the fermion current (C7). For a homogeneous canting field, and at halffilling conditions dx j 0 (x) = N s /2, this back-action renormalizes the external field to −(h h h − g g g/2). In Eq. (D3), we have introduced the single-particle Hamiltonian of a (1+1)-dimensional Dirac fermion with a Dirac mass proportional to the homogeneous Néel field m = g g g s · n n n. In the continuum and T = 0 limits, the product of determinants involving these Dirac Hamiltonians can be expressed in terms of momentum integrals ∏ n,q det (−iω n + h D (q, g g g s · n n n)) = e 1 4π 2 c d 2 k log(k k k 2 +(g g g s ·n n n) 2 ) (D6) where we have introduced Euclidean momentum k k k = (ω n , cq).
This expression, together with Eq. (D3), leads to an effective action which, up to an irrelevant term independent of the Néel and canting fields, reads as follows Despite the non-linearity, one can readily find an exact analytical solution to these saddle-point equations in the limit of vanishing transverse fields h t = 0, where the only possible solutions correspond to cos θ A = cos θ B = 0. There are four possible solutions within θ θ θ ∈ [−π, π) × [−π, π), namely (i) North-pole paramagnet: This solution is θ A = θ B = π/2, where the canting field points towards the north pole (x x x) = (S/a)e e e z , ∀x x x. The corresponding spin coherent state is |g + P = |S, S z,1 ⊗ |S, S z,2 ⊗ · · · ⊗ |S, S z,N s , where |S, m α,i is the common eigenstate of S S S 2 i , S α i with eigenvalues S(S + 1) and m ∈ {−S, −S + 1, · · · , S − 1, S}. Hence, all spin-S particles are aligned towards the north pole. We note that this ordering is not caused by collective effects, but induced by the external longitudinal field, and that is why we refer to this state as a spin paramagnet.
(iii) Néel anti-ferromagnets: These solutions are θ A = −θ B = π/2, or θ A = −θ B = −π/2. In this case, and it is the Néel field which points towards the north n n n(x x x) = e e e z , ∀x x x and south pole n n n(x x x) = −e e e z , ∀x x x, respectively. The corresponding spin coherent states are |g + N = |S, S z,1 ⊗ |S, −S z,2 ⊗ · · · ⊗ |S, S z,N s −1 ⊗ |S, −S z,N s , |g − N = |S, −S z,1 ⊗ |S, S z,2 ⊗ · · · ⊗ |S, −S z,N s −1 ⊗ |S, S z,N s , Hence, all neighbouring spin-S particles are aligned in an antiparallel manner. We note that, in this case, this ordering is caused by collective effects, and the north-or south-pole solutions occur via the spontaneous symmetry breaking of the Z 2 chiral symmetry n n n(x x x) → −n n n(x x x) mentioned in the main text. This is why we refer to these states as anti-ferromagnets. In order to decide which of the above orderings occurs in the system, we can compare the corresponding free energies per unit length f(θ A , θ B ) = 1 L F(θ A , θ B ) = − 1 Lβ log Z(θ A , θ B ), namely Comparing these energies, we find two critical lines h ± that separate the Néel anti-ferromagnet from the two longitudinal paramagnets. The first one is obtained from f (−π/2, −π/2) = f (+π/2, −π/2). Accordingly, for h < h − , the groundstate corresponds to the southpole longitudinal paramagnet. Conversely, for h > h − , we enter into the Néel anti-ferromagnet. The second line is and is obtained from f (+π/2, −π/2) = f (+π/2, +π/2). In this case, for h < h + , the groundstate corresponds to the Néel i h h h s,i · S S S i , h h h s,i = h h h − g g g n i (E3) where f (s) stands for fermions (spins). Accordingly, the fermion tunnel in a potential-energy landscape ε f,i set by the expectation value of the spins (see Fig. 1), while the spins precess in an effective external field, which becomes inhomogeneous depending on the average distribution of fermions. This mean-field approximation requires the obsrevables { S S S i , n i } N s i=1 to be determined self-consistently, and we must deal with a number of self-consistency equations that grows linearly with the number of sites. The self-consistent loop consists in the following steps: We start by setting an initial spin configuration S S S i , and compute the expectations values n i by solving the fermionic tight-biding model (E2) for a given temperature T = β −1 . These mean-field parameters are then used as input to determine the effective external field in the spin Hamiltonian (E3), which is subsequently diagonalised, such that we can calculate the corresponding spin expectation values S S S i for a given temperature T . This process must be iterated until reaching convergence for the free energy, which can be expressed as F MF = F f + F s +C, where F f (F s ) is the free energy of the fermions (spins) and C is the constant term appearing in the Hartree-Fock decoupling (E1). The free energy for the system of fermions can be expressed in terms of the grand partition function where we have introduced the chemical potential µ, and the total fermion number N f = ∑ i n i . For the mean-field decoupled model, this partition function can be readily expressed as where E f,n are the eigenvalues of the quadratic fermionic Hamiltonian (E2). The free energy, which can be expressed as F f = − 1 β log Z f + µ β ∂ ∂ µ log Z f , thus becomes where n FD (µ, β ) = µ∂ µ Z f,n β Z f,n = 1/(e β (E f,n −µ) + 1) is the socalled Fermi-Dirac distribution. For the spins, the number of which is conserved, the free energy can be written in terms of the canonical partition function Z s = Tr e −β H s , For the mean-field decoupled system, the Hamiltonian (E3) can be diagonalised for each spin independently where E s,n are the the (2S + 1) energies of the spins, which are identical for all sites of the chain. The corresponding free energy is Appendix F: Static quasi-particle potentials In this Appendix, we give further details on the calculation of the static potentials between fractional quark-like quasi-particles and between quark bags. In both cases, this is achieved by first pinning two of these quasi-particles to certain locations in the chain, separated by a distance d measured in unit cells (2a). In the first case, since each of them is associated to a soliton in the Néel ordered rotor field, we introduce a small parallel field with This field breaks translational invariance and pines two solitons, and the associated fractional fermionic charges, to positions i 0 and i 0 + d. In Figure 17(a), we choose i 0 = (N s + 1)/3 and ε = 0.02t, and represent the static quark potential V 1 (d) = E 1 (d) − E 0 , where E 1 (d) and E 0 are the energies of a chain with N s sites and N f = (N s + 1)/2 + 1 fermions under the pinning potential (F2) for a given d and i 0 = 10, respectively. The potential changes from repulsive to attractive as we increased gS/t. For the parameters h S = 0.4t, h t S = 0.01t, the transition is located at g c S/t ≈ 1.01.
The static potential between two quark-bag quasi-particles is obtained by first pinning them at a distance d, in a chain with N s sites and N f = (N s + 1)/2 + 2 fermions, using now a parallel field (F1) with The potential V 2 (d) = E 2 (d) − E 0 , where E 2 (d) is now the corresponding energy when two quark bags are pinned at a distance d, is shown in Figure 17(b) for different values of gS/t. In this case, the potential is always attractive.