String tension and robustness of confinement properties of in the Schwinger-Thirring model

Confinement properties of the $1+1$ Schwinger model can bestudied by computing the string tension between two charges. It is finite (vanishing) if the fermions are massive (massless) corresponding to the occurrence of confinement (screening). Motivated by the possibility of experimentally simulate the Schwinger model, we investigate here the robustness of its screened and confined phases. Firstly, we analyze the effect of nearest-neighbour density-density interaction terms, which -- in the absence of the gauge fields -- give rise to the Thirring model. The resulting Schwinger-Thirring model is studied, also in presence of a topological $\theta$ term, showing that the massless (massive) model remains screened (confined) and that there is deconfinement only for $\theta=\pm\pi$ in the massive case. Estimates of the parameters of the Schwinger-Thirring model are provided with a discussion of a possible experimental setup for its realization with ultracold atoms. The possibility that the gauge fields live in higher dimensions while the fermions remain in $1+1$ is also considered. One may refer to this model as the Pseudo-Schwinger-Thirring model. It is shown that the screening of external charges occurs for $2+1$ and $3+1$ gauge fields, exactly as it occurs in $1+1$ dimensions, with a radical change of the long distance interaction induced by the gauge fields. The massive (massless) model continues to exhibit confinement (screening), signalling that it is the dimensionality of the matter fields, and not of the gauge fields to determine confinement properties. A computation for the string tension is presented in perturbation theory. Our conclusion is that $1+1$ models exhibiting confinement or screening -- massless or massive, in presence of a topological term or not -- retain their main properties when the Thirring interaction is added or the gauge fields live in higher dimension.


I. INTRODUCTION
The study of confinement properties in gauge theories is a long-lasting subject of research, with applications in a variety of physical systems ranging from Quantum Chromodynamics (QCD) [1] to effective gauge theories emerging in strongly correlated systems [2]. The origin of the concept of confinement is associated to the fact that no quarks appear as asymptotic particle states. Instead, they seem to be confined inside hadrons and mesons. At date, there is not yet final and detailed description of confinement in QCD. Such difficulties are due to the complexity of QCD, a strongly coupled non-Abelian gauge theory with symmetry group SU (3). Due to this reason, historically, an important role in the understanding of confinement was played by solvable theories in 1 + 1 dimensions -an archetypical example being the Schwinger model [3]. This is a well studied field theory [4], where relativistic fermions are coupled to a U (1) gauge field. It exhibits confinement of the fermions and from this point of view can be seen as a toy model for QCD [5].
The Schwinger model and its multi-flavour generalization can be mapped by bosonization [6] onto massive sine-Gordon models, having a mass proportional to the fermion charge, but frequency β fixed to √ 4π [7,8] (see more references in [9]). In the regime of vanishing charge, and upon addition of an interaction term between the fermions, one obtains a (massless) sine-Gordon model with variable frequency [10]. This model is the Thirring model [11]. In the massless limit its correlation functions are known [12,13]. In the massive case it is solvable by Bethe ansatz [14].
Both Schwinger and Thirring models separately have been heavily investigated for their relevance as a toy model for phenomena occurring in higher dimensions and for the appealing possibility to study solvable/integrable interacting models which may exhibit confinement. In comparison, the model in which 1 + 1-d fermions are simultaneously charged and interact via local interactions -a model to which one may refer as the Schwinger-Thirring model -has been relatively less investigated [7,[15][16][17][18][19][20]. In spite of the fact that, to the best of our knowledge, the Schwinger-Thirring model is not integrable, rigorous results on the mapping on the massive sine-Gordon model have been established [16].
An important way to characterize the confinement properties of models such as the Schwinger and the Thirring one is provided by the study of the string tension via the determination of the energy of the configuration of two probe charges [4]. More precisely, one has to study how the energy E is modified when the two static charges are taken further and further apart. If the energy is divergent the theory is said to be confining. The prefactor of the linear term of the energy with respect to the distance L, with E ∝ L for large L, is called string tension. Due to the phenomenon known as string breaking, some models which have no free charged states (un-der the gauge group) may exhibit finite energy for arbitrarily separated charges. The final result, is therefore, a screening of the charges. For this reason sometimes is preferred to define confinement as the absence of charged asymptotic states. In this paper this position will not be adopted and the phases will be referred to be either in a confined or screened phase. For a general introduction on the confinement problem see, e.g., [21].
In the 1+1-d Schwinger model, one can compute the string tension between two charges. It is found to be finite if the fermions are massive and vanishing if they are massless [4,5]. It is therefore concluded that for massive fermions one has confinement, and for massless fermions there is the occurrence of screening [5].
This paper focuses on the confinement properties of the Schwinger-Thirring model in order to investigate their robustness and the role played by the one-dimensionality of the fermions. By studying the string tension, the effect of the Thirring interactions is addressed first. Secondly, the system where the fermions live in one spatial dimension and the gauge fields are defined in D = 2, 3 dimensions is explored.
This work is directly motivated by the following two reasons. From one side it is intended to determine whether a 1 + 1 model exhibiting confinement or screening -massless or massive, in presence of a topological term or not -maintains its properties when interaction is added and especially when the gauge fields are allowed to live in higher dimension. The latter question is especially relevant since the confinement property of the Schwinger model could be intuitively explained by the fact that, at the classical level, the energy between two-point particles grows linearly with the distance for a gauge field in 1 + 1 dimensions. When naïvely applied, this argument would lead to conclude that there should be deconfinement when the gauge fields are living in three dimensions. Therefore, to test the robustness of the confined phase, the case where gauge fields are living in higher dimensions (2 + 1 and 3 + 1) while the fermions remain in 1 + 1 is studied. This is the equivalent for the Schwinger model of the Pseudo-QED, in which the fermions are confined in a 2D plane while interacting with the electromagnetic field living in the 3D space [22,23]. The Schwinger-Thirring model in which the gauge fields live in higher dimension is referred as the Pseudo-Schwinger-Thirring model.
The second motivation for our study comes from recent theoretical [24][25][26][27][28][29][30][31][32][33][34][35] and experimental [36] progress on the emulation of gauge theories with ultracold atoms and trapped ions. In general, quantum simulators [37,38] are built up to emulate quantum mechanical systems by properly shaping the system dynamics using external fields. Due to the general complexity of the many problem of quantum mechanics, a quantum simulator could provide answers to long-standing problems in physics. The realization of quantum simulations of interacting theories could help the understanding of target models, the validation of analytical or numerical techniques, and the exploration of physical phenomena which are not currently reachable by other approaches. Ultracold atoms and trapped ion setups provide a variety of reliable different tools to perform simulations of different many-body phenomena [39,40].
One of the most challenging goals certainly concerns the implementation of gauge theories where non-perturbative phenomena, such as confinement [41], occur. This motivates the study of simpler models which may exhibit such phenomena. As mentioned earlier, the Schwinger model is probably the simplest, non-trivial gauge theory involving fermions one can think of, and, more relevantly, it also exhibits confinement. It was also the target of the first quantum simulation of a lattice gauge theory addressing the real-time dynamics on a few-qubit trapped-ion quantum computer [36].
Departing from the original model Hamiltonian, two variations are experimentally accessible (if not unavoidable) according to recent implementation schemes: the addition of tunable Thirring interactions, and the engineering of a gauge field living in D = 2, 3. As mentioned above, these are exactly the two scenarios focused here. A further interesting ingredient, the topological θ-term, can be added as well. In the presence of this extra term, deconfinement is possible for θ = ±π, while the system retains its confining character for any other angle in between [8,42].
As previously reminded, it is well known that the massless Schwinger model is in the screened phase, while the massive one exhibits confinement [5]. Here it is shown that for massless fermions the screening phase survives when the four-point local interaction term is turned on. The confined phase remains present, as well, when a finite mass is turned on in the presence of interactions, as in the Schwinger model. In both cases (massless and massive) the string tension does not depend on the Thirring interaction coefficient g. Consistency of the theory requires, for both cases, that g > −π. When the gauge fields are allowed to live in higher dimensions (2 + 1 or 3 + 1), giving rise to the Pseudo-Schwinger-Thirring model, the massless (massive) model will be shown to remain screened (confined). In particular, it is shown that in leading order on the mass, the string tension for gauge fields in 1 + 1 and 2+1 dimensions is the same. The same scenario occurs in presence of a topological θ-term. Our results shows that both Thirring interactions and gauge fields living in higher dimensions do not alter the confinement properties of the Schwinger model with the θ term, and for θ = π the massive (massless) model remains deconfined (screened).
The paper is organized as follows. In Sec. II, the continuum and lattice formulations of the Schwinger model are introduced, briefly discussing the quantum link formulation. The generation of a Thirring-like interaction based on nearest-neighbour density-density interactions between nonrelativistic tight binding is presented. In Sec. II A the definition of the string tension is reminded and the criterion that will be used to determine confinement properties is stated. In Sec. III, it is shown that the confinement properties of the Schwinger model are retained under the presence of a Thirring-like interaction. In Sec. IV an estimate of the parameters of the simulated quantum system is provided. This is done using, as reference, the proposal Ref. 27, where a Thirring-like interaction appears naturally as a byproduct of the proposal. The difficulties arising in such implementation are discussed and the achievable range of parameters as a function of the amplitude of the external superlattice poten-tial is explored. Finally, in Sec. V, the Pseudo-Schwinger-Thirring model featuring the gauge fields in higher dimensions is addressed analysing both the massless and massive theories. The conclusions are drawn in Sec. VI, while additional, more technical material is presented in the Appendices.

II. IMPLEMENTATION AND DEVIATIONS FROM THE SCHWINGER MODEL
The Hamiltonian for the Schwinger model in the axial gauge A 0 = 0 is given by: In lattice field theories it is possible to provide a regularization of their continuum counterpart by introducing a lattice spacing a s . As a consequence, different lattice formulations may correspond to the same model in the continuum limit a s → 0. An example is given by the Kogut-Susskind (KS) Hamiltonian [43] which, in the continuum limit for the U (1) gauge theories, gives the QED Hamiltonian in the axial gauge. For the case of one spatial dimension, the KS Hamiltonian is given by where c n annihilates a fermion in the lattice site n. In the KS model the spinor degrees of freedom are encoded in the spatial coordinates (staggered fermions) [44] and the continuum limit is taken by sending the lattice spacing to zero a s → 0.
More precisely the continuum variables are identified through c n / √ a s → ψ up (x) for n even and c n / √ a s → ψ down (x) for n odd, to form the two dimensional spinors. E (x) denotes the electric field. Furthermore the representation of the gamma matrices is fixed to be γ 0 = σ z and γ 1 = iσ y where the σ's are the Pauli matrices. The gauge operators obey the commutation relations [U m , E n ] = eδ mn and U † m , L n = −eδ mn , with the remaining commutation relations vanishing. The parameter e is the gauge coupling. The lattice model, besides regularizing the continuum theory, provides also a bridge towards the experimental implementation.
One of the first and well documented difficulties in implementing such model consists in reproducing the infinite dimensional space spanned by the gauge fields present on the links (the electric field in each link varies from −∞ to +∞). An alternative approach consists on truncating the Hilbert space making it finite and more suitable for quantum simulation. Such models are known as quantum link models [45][46][47]. In a quantum link model, the link algebra described above is replaced by the algebra of angular momentum [L i,m , L j,n ] = iδ nm ε ijk L k,n . Each representation of this algebra, identified by the spin magnitude S, provides a finite Hilbert space. In the limit of S → +∞ the KS Hamiltonian should be recovered. The commutation relations between the U m 's and E m 's are exactly realized by identifying E n → eL n and iU n → L +,n /S (S + 1) (where L +,n = L x,n + iL y,n ). However there are now extra non-zero commutation relations: U m , U † n = 2δ mn E m / [eS (S + 1)], reflecting the request that H KS is recovered when S → +∞. This could constitute a possible drawback on the experimental implementation of these models. This subject was extensively studied in [35] where it was shown how the finite dimension of the corresponding Hilbert space may deviate from the infinite dimensional case. Qualitatively, the finiteness of the quantum links plays a small role and, as expected, the results converge to the QED result as the value of the total spin is increased.
Another problem that may emerge is the possibility of nearest-neighbour density-density interaction. In fact, as in [27], such term may be present as a byproduct of the implementation scheme. In more general cases, terms of these type may appear quite naturally in implementation schemes where gauge invariance is imposed via energetic constraints, due to the fact that such terms are indeed gauge invariant.
From a different perspective, the Schwinger-Thirring model can also host interesting physics which may motivate direct implementations of it in its own. In the scheme of Ref. [27], the four-Fermi interaction is repulsive and it is not possible to reverse the sign of interaction. The reason of such contribution directly comes from the Fermi statistics. For this reason a possible path towards an implementation where such interactions are tunable may include, for example, an extra species of bosons with correlated hopping with the fermions through all the lattice, producing an analogous term with opposite sign. This term could give extra control on the sign and strength of this interaction.
Building on this intuition, the consequences of the presence of such term in the gauge theory are investigated. This is done in a general setting where it is admitted a θ-term on the Lagrangian. The lattice Hamiltonian is then modified to be (3) where n n = c † n c n . The presence of l 0 corresponds to a background electric field which in the continuum limit gives rise to the θ term with l 0 = θ/2π. The other extra term is a nearest-neighbour density-density interaction. This term, which is of the form λ n n n n+1 , scales with a s to get a finite contribution in the continuum limit, as it can be seen by taking into consideration that n → 1 2as dx and the four operators c x bring a factor of a 2 s forcing a pre-factor∼ a −1 s . In this way the possible terms of the form ψ∂ψ are of higher order and the remaining part in terms of spinor com- When comparing it to the Hamiltonian of the Thirring model g 2 dx ψ γ µ ψ ψ γ µ ψ , which results in a term 2g dxψ † up (x) ψ up ψ † down (x) ψ down (x) one concludes that the correspondent lattice parameter is given by 2g/a s .
The inclusion of the l 0 term, and equivalently of the θ term, corresponds to a shift of the background electric field of the vacuum. This is a consequence of the presence of a linear term ∝ L n . In the scheme [27], that will be discussed in Section IV for the estimate of the achievable range of parameters, this corresponds to create an extra imbalance between the chemical potentials of the boson species.
A. Confinement and screening in 1 + 1 gauge theories Here the existence of confinement is characterized by computing the string tension σ between two external charges added to the system. Since such characterization of confinement properties in Schwinger-Thirring models is the main subject of the present paper, a pause is made here to relate confinement, deconfinement and screening to the string tension σ. This quantity is defined as the constant of proportionality between the energy T of two added external charges and their distance L: when, for large L, σ is positive one has confinement, while σ < 0 one has a deconfined phase. In the case in which σ = 0, one then studies the limit for large L of the energy T : if it is infinite and positive/negative then one has respectively confinement/deconfinement. When σ tends to 0 and |T | is not diverging, then the ratio T with /T without can be considered, where T with (T without ) is the energy with (without) the fermionic fields. If this ratio is vanishing, one has a screened phase, while in the other cases it is not possible to conclude about confinement, deconfinement or screening just looking at the energy T . Instead one should look at the behaviour (and the poles) of correlation functions in order to determine the presence/absence of charged asymptotic states [4]. This way of characterizing confinement properties covers the cases to be considered here.

III. ROBUSTNESS UNDER THIRRING INTERACTIONS
The problem of confinement in the Schwinger and Thirring models in 1 + 1 dimensions can be addressed through bosonization [6]. The general procedure adopted here follows closely [17] and it will be be the basis for the results presented in Section V where the gauge fields are living in higher dimensions.
The continuum Lagrangian in Euclidean time for the Schwinger-Thirring model reads: For g = 0, the model (5) is the Schwinger model with the θterm, which is known to exhibit (partial) deconfinement only for θ = ±π [8,42]. In turn, for e = 0 one has the the Thirring model. Such theory makes sense only for g > −π, as discussed below, and it can mapped order by order in the mass m to a sine-Gordon model [10]. In the following it is showed that these results remain valid when both parameters e and g are finite in the Schwinger-Thirring model. In particular the Thirring term does not play any role on confined/deconfined phases of the model.
The quartic fermionic interaction in (5) can be re-expressed via an Hubbard-Stratonovich transformation. This amounts to replace in the Lagrangian g , with now the integration being performed over B µ as well. Even though the parameter e is entering the transformation, the integration of the fields B µ produce a result independent of e and such introduction of the matter-gauge coupling is present for convenience only. In order to isolate the coupling with the fermion fields the variable transformation A µ → C µ = A µ + B µ is performed. The resulting theory reads: where the indices b, c refers respectively to B and C fields. The result is then a standard Schwinger model plus an extra boson field coupled to the gauge field. It is worth noting that this is a peculiar type of interaction which is independent of the parameters e and g. These parameters control the interaction between the gauge field and the fermions, and the mass of the field B 2 µ . The gauge transformations are encoded in the usual way in the fields ψ and C µ . In turn the B µ field does not transform under gauge transformation. Exploiting gauge freedom, one can pick the Lorentz gauge where ∂ µ C µ = 0 and therefore parametrize the field C as as C µ ≡ −iε µν ∂ ν ϕ which is a general form for divergenceless fields in two dimensions. The field B has be parametrized with a gradient part too: It turns out that the field χ decouples from the remaining ones and therefore can be left out for this analysis.
In turn, the massive term is mapped to −m ψ ψ cos 2eϕ + iψγ 5 ψ sin 2eϕ . Furthermore, induced by the chiral anomaly, an extra term arises with the form − e 2 2π (∂ µ ϕ) 2 . The Lagrangian on ψ is now decoupled from the rest of the system and it can be mapped to the bosonic action 1 π ϕ the full Lagrangian is then quadratic on ϕ and ϕ and it can be integrated out.
The full Lagrangian after this procedure is given by: The field χ only appears in the kinetic term and is not coupled to the other fields as announced before. The θterm of ϕ can be re-cast in cosine form by transforming Crucially, the coupling between ϕ and ϕ has just a prefactor 1 in front of it, making the terms ∂ 2 ϕ 2 cancel out. This result is due to the fact that one of the gauge fields is actually a fictitious field derived from a Thirring interaction which finally controls its mass and not the interaction with the other gauge field. The next step consists in integrating out the field ϕ . This field is also related to a θ-term with the opposite sign, inducing a new transformation ϑ → ϑ − θ/ √ 4π (opposite to the one made above). However, the ϑ field acquired a mass, so as expected the dependence on θ is not erased but instead is explicit in the term e 2 2π ϑ − θ/ √ 4π 2 . Nonetheless this transformation is useful to perform the integration on the ϕ field, which appears in the form e 2 2g (∂ µ ϕ ) 2 − e √ π ∂ µ ϑ∂ µ ϕ . The integration brings the Thirring contribution to the bosonic action g 2π (∂ µ ϑ) 2 . Finally, the usual mass term is restored by The Lagrangian finally reads: The Lagrangian (8) shows the contribution of both the Schwinger and the Thirring models. No mixing between the g and e couplings occur. It is worth noting that the above Lagrangian corresponds exactly to the Lagrangian obtained from (5) if one applies the mapping obtained from bosonization of the free massless Dirac fermions integrating out the gauge field. More precisely, this amounts to replace the pure fermionic action by the respective sine-Gordon model and replace the current byψγ µ ψ = 1 √ π ε µν ∂ ν ϑ both in the coupling to A µ and in the Thirring term contribution. After translating the bosonic field ϑ by θ/ √ 4π and integrating the gauge field, the Lagrangian (8) is obtained. The fact that such substitution holds is not obvious and it guarantees that the final result is correct.
From the bosonized Lagrangian (8) it follows that, as in the Thirring model, it is required g > −π in order the model itself makes sense. The existence of a threshold for a minimum of g is expected from the lattice theory. For g → −∞ the dominant term on the Hamiltonian is a nearest-neighbour strong repulsion which will induce a phase separation in the system. For that limit the field theory, or, in other words, the continuum limit, cannot be taken.
Regarding the screening and confinement properties of the model, the simplest case corresponds to the massless theory where µ = 0. The propagator for the bosonic theory is given by From (9) one can compute, for instance, the two-point function of the scalarψψ and pseudoscalarψγ S ψ fields, showing that there are no charged fermions in the spectrum. The derivation of this result can be found in the Appendix A.
For the massive case a perturbative analysis can be found on Appendix B where the string tension is computed. The computation of the string tension for the massless case of the Schwinger model is treated in Section V as a particular case of the general construction with gauge fields in D + 1 dimensions. The screening of the external charges becomes then explicit with an exponential decay of the energy with the distance between charges. As a consequence, the string tension is vanishing.
The results of this Section, and in particular the form of the Lagrangian (8), clearly show that the confinement and screening phases of the Schwinger model, respectively for the massive and massless cases, are not altered by the Thirring interactions terms. In the next Section an estimate of the lattice parameters of the Schwinger-Thirring model in setups of ultracold atoms is provided, while in Section V the robustness of confinement when the gauge fields live in higher dimensions is discussed.

IV. ESTIMATES OF LATTICE PARAMETERS FOR ULTRACOLD ATOM SETUPS
A number of proposals have been put forward to simulate the Schwinger model. These include both proposals in which the gauge symmetry emerges as a symmetry of the low energy effective theory [27,32,48] or as an exact symmetry of the system [25,35]. An important point to be stressed is that fine tuning of parameters should be neither required or crucial since small inaccuracies in the experimental implementation could be spoil the validity of the quantum simulation. In this Section the issue of the deviations from the desired Hamiltonian due to the presence of a nearest-neighbour density-density interaction is discussed. It is argued that the presence of such four fermion interaction term corresponds to a Thirring interaction. A scheme for simulating the Thirring model was also put forward in [49].
This Section also provides an estimate of the values of a different parameters in a possible implementation with optical lattices [50]. The main goal here is to provide a quantitative estimate of the energy scales in the system, and in particular, of the role of Thirring terms. To this end, and in order to have a specific example at hand, the model described in [27] was chosen. There, the density-density interaction appears explicitly. This model makes use of one species of fermions and two species of bosons and it builds the quantum links using the Schwinger representation. The fermions are hopping between all lattice sites while odd links are associated with one species of bosons, and even links with the other. Each boson is then only allowed to hop between its designated link. In concrete, at any lattice site n, one can have bosons of both species. If n is odd then the species 1 can only hop to n + 1 and the species 2 to site n − 1. The opposite happens for n even. This is illustrated in Figure 1 where both species of bosons and fermions are represented around an even site.
The quantum links are realized through the Schwinger representation where with σ indicating the bosonic species (1 or 2) that has to be consistent with the parity of the link in question. In this language the generator of gauge transformations can be written as where n F n indicates the number of fermions and n σ n the number of bosons of species σ in the site n. Therefore if writing an Hamiltonian of the form H = H 0 +U G 2 n where the energy scale U is much larger than the energy scales of H 0 , one always obtains, in perturbation theory, a gauge invariant low energy Hamiltonian. When considering H 0 to be the sum of single particle Hamiltonians of the three considered species, the effective Hamiltonian is of the form (3).
Adopting the typical notation used for ultracold atoms [50], one has to evaluate the hopping parameters t α of the (nonrelativistic) ultracold atom mixture and the parameters of the lattice Hamiltonian (3). The index α is here used to denote the fermions or one of the two bosonic species: α ∈ {F, 1, 2}. To establish the connections between the parameters t α and the lattice Hamiltonian (3), one has to restore and c in the Hamiltonian, which corresponds to add c to all terms except the mass term which gets a c 2 . The parameters of the KS Hamiltonian are a s (which is not the lattice spacing of the cold atomic optical lattices), the electric charge e and the Thirring term g. The kinetic term is characterized by t B t F /U , the pure bosonic term by t 2 B /U and the nearest-neighbour density-density term by t 2 F /U . The connection with the parameters of the KS Hamiltonian is given by from which one conclude that the parameters cannot be varied independently. A key and challenging aspect of this implementation scheme is that the interactions between the different atomic species shall satisfy the conditions U 11 = U 22 ≡ U and 2U 12 = 2U 1F = 2U 2F = U [27], where U αβ is the onsite interaction parameter between species α and β [50]. The related problem is that the proposal requires the bosons sitting in asymmetric minima. This asymmetry may change the Wannier functions entering the evaluation of the U parame-ters [50]. This may lead to interactions of the form U +/− αβ where the labels +/− indicate the relative (+) and absolute (−) minima of the optical lattices involved in the experimental implementation. However, the existence of different parameters U +/− is not crucial. Gauge invariance at low energies is obtained in perturbation theory for large U so as long as both U +/− remain larger than the other parameters of the model the condition G x |ψ = 0 still holds. Since the goal of this Section is to provide an estimate of the parameters of the underlying lattice Hamiltonian, this asymmetry is disregarded, since it can be made small. Details about this point are given in Appendix C. Furthermore, to satisfy the conditions U 11 = U 22 ≡ U and 2U 12 = 2U 1F = 2U 2F = U one still has to match the different interactions between atomic species. Again, the system can be shown to be robust under small enough deviations from these conditions. Possible deviations lead to U G 2 x → U G 2 x +U ∆U αβ n α x n β x , where ∆U αβ is the deviation of the interaction between species α and β from the desired value. The fundamental requirement to obtain a gauge invariant theory is still valid as long as ∆U αβ /U 1. To estimate the lattice parameter, a mixture of 52 Cr is used as reference. This particular choice is related to the fact, as explained below, that their scattering lengths approximately have the required interaction strength between the two bosonic species foreseen by the proposal. The scattering length between pairs with total angular momentum equal to 0 is a b0 30 − 50 and for pairs with total angular momentum equal to 4 is a b4 58 ± 6, with a 0 the Bohr radius [51,52]. With a b0 = 30a 0 (different species) and a b4 = 60a 0 (same species) one has the required relations between U 11 , U 22 and U 12 . For the fermionic species it is assumed that a tuning of the interaction is possible in such a way that the remaining conditions can be obeyed. For the present calculation it is assumed, as well, that a 1F = a 2F = 30a 0 . Notice that the presence of bosonic-bosonic interactions does not have the fermionic-fermionic counterpart and that fermionic and bosonic densities are different (see the difference between the two equations in C10). However, as long as S is not very large, it is expected that the effect is not significant. Here the calculations will be performed for S = 1 and this assumption proves to be enough. In order to increase S one should carefully adjust a 1F and a 2F . Mixtures involving alkalineearth-like atoms may prove to be even more convenient, due to the possibility of easily engineering species dependent lattices there.
The Wannier functions are assumed to be Gaussian with a width σ i α , depending, also, on the direction i and determined variationally by minimizing the Gross-Pitaevskii energy (see more details in [53,54]). The width of the Gaussians in the directions perpendicular to the optical lattices is considered to be σ ⊥ , assumed here to be the same for all the species (fermions and bosons). The potential felt by the particles is characterized by an amplitude V 0α , which controls the height of the barrier, and by V 0α ∆ α , the off-set controlling the difference of energy between minima (with lattice spacing a).
It is assumed S = 1, V 0F = V 0B = V 0 and ∆ F = ∆ B = 0. For the lattice spacing it is taken a ∼ 0.5µm. The process is detailed in Appendix C. In the following, we specif-ically address the case of a 1D optical lattice, where transverse confinement in the direction perpendicular to the wire is shallow. Similar ideas are applicable in the 3D optical lattice case. Two different choices for σ ⊥ were considered: σ ⊥ = 2a and σ ⊥ = 5a. For smaller values, the variational approach is not expected to be accurate while, for larger values, very large potential depths are required to enter in the necessary perturbative regime. The main result is illustrated in Fig. 2, where the effective parameters t B t F /U , t B t F /U and t B t F /U are plotted in a single curve (one for each σ ⊥ ). This is because, at this scale, they are indistinguishable. Furthermore, with these parameters, t α /U becomes approximately 10 −1 at V 0 = 130E ref and V 0 = 180E ref , respectively. Note that, when compared to the usual recoiled energy, E ref is π 2 times larger. By adjusting some parameters, other choices are of course possible. Such values were chosen in order to give an illustrative example. For these range of parameters, one finds that the relation 2U 1F = 2U 2F = U is satisfied very accurately with an error inferior to 10 −2 %. At the same time, the ratio t F /t B remains always very close to 1.
In conclusion, the range of parameters one may reasonably have access to corresponds to −g ∼ 0.25 and e 2 a 2 s ∼ 0.5 bounded by the condition −8/g = e 2 a 2 s . The typical energy of the processes of the effective cold atomic system is of the order t α t β /U, V ∆ ∼ 10 −3 E ref . Moreover, in this regime the value of the three terms t α t β /U become indistinguishable. Finally it is observed that the strength of the mass parameter is less constrained and it can be changed via the quantity ∆. This parameter was set to zero during the parameter estimation but, as long as it is not too large, the approximation remains valid.

V. SCREENING WITH GAUGE FIELDS IN D + 1 DIMENSIONS
In this Section the Schwinger-Thirring model in presence of a gauge field living in higher dimension is considered. As briefly mentioned in the Introduction I, naïvely one would expect to find deconfinement is in D = 3 spatial dimensions, one may anticipate to find features of normal QED with the particularity of the "electrons" being restricted to one spatial dimension. Here it is showed that the situation is more subtle.
Models exhibiting dimensional mismatch between relativistic fermions and gauge fields have been studied and applied in a variety of situations: in the context of graphene [23,55,56] and transition-metal dichalcogenides [57,58], for topological insulators both for 2 + 1 [59] and 1 + 1 fermions [60], for 1 + 1 fermionic models [60][61][62] and as a source to generate effective short-range [63] and long-range [62] interactions via the dimensional mismatch. Besides the direct relevance of these models to experimental systems, such as graphene in a 3D electro-magnetic field [22], these models are also promising from the points of view of theoretical study and of experimental implementation perspectives. The most immediate application is given by the realization of an intermediate step towards the simulation of increasingly complicated gauge theories. For this, the necessary correlated hopping of bosons and fermions is restricted to the line, while the gauge fields still live in higher dimensions. In terms of implementation complexity, these models are expected to be more complicated to implement than the Schwinger model but still simpler than QED. The same applies if the gauge fields are non-Abelian.
In the following the basic construction of theories in which the gauge fields live in D + 1 dimensions is revised. The general form of the Lagrangian is given, in Euclidean time, by: The fermions are considered to live in the lower dimensionality d + 1, which is made explicit in the matter Lagrangian L d+1 M . The gauge field lives in the higher dimensionality D + 1. The D + 1 current is taken to be: In (10) the term L GF corresponds to the Fadeev-Popov gauge fixing term which is given by where different choices of ξ correspond to different gauges. The Feynman gauge, where ξ = 1 and one has a diagonal propagator G µν = 1 −∂ 2 δ µν , is adopted. The theory (10) can be suitably formulated only in the lower dimension d without explicitly invoking higher dimension. The effect of the higher dimensionality is encoded in a modified kinetic term for the gauge fields. For the case of d = 2 and D = 3 it goes under the name of Pseudo-QED [22]. Details on the general construction can be found in [62]. The general form of the Lagrangian is given by: where the operatorM D is given in terms of the propaga-torĜ D of the gauge fields according to the relationM D = In the previous expression ∂ 2 is the Laplacian in d + 1 dimensions. The explicit expression forM D (orĜ D ) depends on the higher dimension D but all the fields are now exclusively in d + 1 dimensions. The calculation of the effect of external charges on the system can be done as for the Schwinger case. This amounts to introduce an extra contribution −iQA µ j µ ext where Q is the absolute value of the two opposite external charges. This external current can be written in the form j µ ext ≡ ε µν ∂ ν K and it can be eliminated by a chiral transformation. The variable change corresponds to ψ = e iQKϕγ5 ψ where, again, one should take into account the chiral anomaly. Notice that the modified kinetic term of the gauge field has no effect on the procedure. Once the gauge fields are integrated out, the resulting bosonic theory is given by: With the field transformation the coupling between the field K and the bosonic field is translated to a sine-Gordon form. The resulting Lagrangian reads: where M D and α D can be seen as operators acting on φ and K D is a simple space-time function. After some algebra one can write α D = 2eF D K and The unperturbed theory, i.e. the theory with no external charges, can be easily recovered by setting Q = 0.
Despite the non-locality, the above Lagrange is still translational invariant in space and in time. The latter gives rise to the conservation of energy. The total energy can be computed through the energy-momentum tensor. This will be, in general, rather complicated. Nonetheless it is still possible to compute the difference of energies, since the more complicated terms cancel out (for the massive case one can also do it in first order in perturbation theory on the mass).
In Appendix D the construction of the energy-momentum tensor for theories with higher derivatives is revised. In the present case there is an arbitrarily high number of derivatives. The energy is given by the integral in space of the T 00 component of the energy-momentum tensor, Eq. (D7). When quantized the fields are promoted to operators and take T 00 in normal order so the total energy is given by E = dx : T 00 (x) : . The difference of energy as a result of introducing external charges can be written as: where it was denoted by T 00 Q (x) the energy-momentum tensor when external charges are also present according to Eq. (14).
The analysis start in the massless case moving then to the small mass limit.

A. Massless fermions
For massless fermions µ = 0. This case is much simpler to analyse since the effect of external charges is only in the term Q 2 K D , completely decoupled from the fields. However complicated the energy momentum tensor obtained from the system without external charges is, the same tensor results from the presence of external charges with an additional space-time function independent of the fields. In other words, the first term of the tensor in (D7) is not affected by the external charges while the second is just "translated" with no operator content [see Q 2 K D in (14)]. Therefore for the massless case one obtains ∆E m=0 = Q 2 dxK D (t = 0, x). Recall that K D = 1 2 ∂ µ KF D ∂ µ K. Each ∂ µ K encodes two Dirac deltas corresponding to two different external charges as described above. Therefore in this expression there are included interactions between the charges, corresponding to pick the Dirac deltas at different points, and "self-interactions", corresponding to pick the same Dirac delta in both K's. The latter ones are independent of L and therefore do not account for actual interaction between different charges. In the following they are then neglected. By performing the implicit integrals on the definition of K D and making use of the fact that ∂ µ K is time-independent, the energy can be then written as where F D (k 0 , k 1 ) are the Fourier components of F D . In the following these integrals are computed explicitly for D = 1, 2, 3.

Massless D = 1
This case corresponds to haveĜ 1 = 1/ − ∂ 2 and F 1 (0, k) = 1/ e 2 /π + k 2 . The integral results in: Without the presence of the fermion field (or turning off the coupling e = 0), the resulting energy would exhibit a linear growth with the distance. The exponential decay present here is a result of pair production that screens the external charges. This shows explicitly the charge screening known for the massless Schwinger model. As discussed in Section III, the same happens in presence of the Thirring interactions between fermions.

Massless D = 2
The function F 2 is given by F 2 (0, k) = 1/ e 2 /π + 2 |k| . Again the integral can be performed explicitly: The functions Ci (cosine integral) and Si (sine integral) are respectively given by Ci(x) = − +∞ −x dt cos t/t and Si(x) = x 0 dt sin t/t. In the limit of L → ∞ the cosine integral goes to zero and the sine integral converges to π/2. As a result also here the energy goes to zero (as 1/L 2 ) as the distance L increases despite the pure gauge theory is exhibiting a logarithm increase of the energy with the distance.

Massless D = 3
For the three-dimensional case one has to introduce an UV cut-off Λ in order to regularize the integral over the extra dimensions where the gauge field lives. The resulting function F will be dependent on this cut-off: it is found The integral (16) requires a careful study. Within a change of variables it can be written as: The dependence on the distance is now isolated in the prefactor 1/L, since the remaining was absorbed into the cut-off Λ = LΛ. In this expression the screening due to pair creation is evident: setting e = 0 the integral (19) simply givesΛ −1 (in the large cut-off limit) and what it remains is the expected Coulomb energy: Q 2 /4πL. When e acquires a finite value, one is coupling the gauge fields to the fermion fields and the pair production starts. This is made explicit in the integral since this extra positive term in the denominator will decrease the absolute value of the integrand.
It is now shown that for any finite charge e total screening occurs and actually ∆E D=3 = 0 in the large cut-off limit. The integral can be broken into pieces: +∞ 0 = n 2π(n+1)/Λ 2πn/Λ . Let us denote the non-oscillatory part by f (q) = log 1 + q −2 / 1 + e 2 /4π log 1 + q −2 . Since this function is not singular except for q = 0, most of the integrals to be computed vanish in the large cut-off limit. This can be seen by integrating by parts which will bring powers of Λ to the denominator. At lowest order, if the functions f has finite derivatives, the leading term goes likeΛ −3 . After the above procedure, the only part that remains is the case n = 0. This is treated by observing that f is strictly decreasing in the interval of integration. Thereforẽ as one can see by replacing the value of the function f by its maximum value in the interval (4π/e 2 ) when the cosine is positive and by its minimum value when the cosine is negative. Since the function is continuous one can make f 2π/Λ as close as desired to 4π/e 2 by increasingΛ, therefore the bound goes to zero. Since the integral is positive this shows that the energy goes to zero. One can finally write: When the coupling between the gauge fields and the fermions is turned on, the fermionic fields react to the presence of external charges by starting pair production. Remarkably they are able to completely screen the external charges. This suggests that when the gauge field is in 3 + 1 dimensions the fermions become more effective at screening external charges than at 2 + 1 or even 1 + 1. For the latter case the energy decreases exponentially with the distance while here it is zero for any distance.

B. The massive case
The massive case is more complicated and the effect of the external charges is no longer decoupled from the fields. Without the external charges, the interaction is still local and all the non-locality is on the kinetic term. When the external charges are introduced, the non-locality is carried over to the interaction via α D . This means that the insertion of external charges will modify both terms present in (D7).
Let us start when the system has no external charges. By performing first order perturbation theory in the mass of the fermionic theory, the ground-state has the structure: |Ω 0 = |0 + µ |1 . The state |0 is the vacuum of the massless theory, which is still a quadratic theory. The normal ordering is taken with respect to this state. When going to the system with external charges, even though the quadratic term is modified, the modification is of order µ so that the ground state of such theory is, at the lowest order in perturbation theory, given by |Ω Q = |0 + µ |1 , where |0 is the same vacuum state of the massless theory and the first order correction was modified due to the presence of external charges. The normal ordering is then taken with respect to the same state in both theories.
One can now analyse the energy-momentum tensor (D7). With no external charges T 00 = T 0 − µ cos √ 4πφ , where all the contributions independent of µ were inserted in T 0 . In the presence of external charges, this is modified to T 00 = T 0 + µT 0 − µ cos √ 4πφ + Qα D + Q 2 K D where µT 0 is the order-µ term obtained from the first part of (D7). Explicitly one has: The presence of this term is due to the fact that the non-locality is carried over to the interacting part, proportional to µ, by the presence of external charges. As a result one has at first order of perturbation theory: Due to the normal ordering, the only non-vanishing term from the Taylor expansion of the first cosine is 1. All the others average to zero in the ground state. The same kind of argument holds for cos √ 4πφ + Qα D , where only cos (Qα D ) is non vanishing. Finally, it is observed that, by construction, any term ofT 0 has always at least one φ, as it is clear from 22. Therefore it averages to zero in the ground state in presence of normal ordering. The result is then: Intuitively, from the above expression one expects a finite string tension when Qα D is "mostly" non multiple of 2π between −L/2 and L/2 and "mostly" multiple of 2π outside this interval. This, as it will be seen explicitly in the following, is what happens for the derived α D in the different dimensions. From the definition of α D , and following the same path used for K D when deriving (16), one can write: Again the different dimensions are considered separately in order to compute the increment to the energy due to the presence of the mass. From this point on the notation α D (x) ≡ α D (0, x) is adopted.

Massive D = 1
This is the well known case studied in detail in the literature [4]. In this Section the computation to retrieve the expected results in the present formalism is performed, serving also to set the scene and the notation for the subsequent computations in D = 2 and D = 3. One has This integral can be calculated explicitly giving: (27) It turns out that, in order to compute the string tension, it is enough to work out the limits |x| L/2 and |x| L/2, as it is shown below. By inspecting directly the function (27) one finds: This is enough to compute the string tension from (24) even without computing exactly the integral (24) itself. Indeed, making use of the fact that the integrand is symmetric over x → −x, the integral can be broken in three parts: . The value of x 0 is fixed to guarantee that exp − e √ π (L/2 − x 0 ) 1. Within this limit one can compute the first and the third integrals using the asymptotic expressions in (28) obtaining: where one can explicitly see the linear growth in L. Note that x 0 can be chosen independent of L for large enough values of L, and the corresponding term ∝ x 0 in ∆E m actually does not grow with L. This reflects the contribution arising for the two regions close to the charges which is independent of their distance (if large enough). Furthermore, the remaining integral is bounded by values independent of L. Explicitly, substituting the cosine by −1 one has an upper bound of 4x 0 and substituting the cosine by 1 one has a lower bound of 0.
One can therefore conclude that the linear behaviour in L is exclusive of the term ∝ L and one can finally write: where the dots indicate some bounded dependence on L. The string tension reads explicitly: This is the well known result [4], reviewed also in Appendix B, obtained here by a careful analysis of the energy excess due to the presence of external charges. For higher dimensions the same procedure shall be followed in the next Sections.
(34) The same procedure from before is followed, studying the limits of |x| L/2 and |x| L/2. The result obtained is actually the same here. For |x| L/2 this is seen by noticing that inside the cosine and sine integrals one can replace X ± by L and take the large L limit. Then one can replace the cosine integral by zero and the sine integral by π/2. For the second case X ± is replaced by ±2x inside the cosine and sine integrals [note also that Si (−y) = −Si (y)]. Then one finds the same result as for (28) and the results for the 1 + 1 case translates directly to 2 + 1. In particular the string tension is the same at this order in perturbation theory in the mass: Even though te presence of confinement in itself is not a surprise for this case, it is interesting to observe that the the resulting string tension, at this order in perturbation theory, is independent of the gauge fields living in 1 + 1 or 2 + 1 dimensions.

Massive D = 3
Finally the case in which the gauge field lives in 3 + 1 dimensions is considered. The function α reads: The same kind of analysis, that was done for the other two cases, can be followed here. This consists in looking at the two limits where x is much smaller or much larger than L. The three main steps consist on writing 2 sin (kL/2) cos (kx) = sin (k (L/2 + x)) + sin (k (L/2 − x)), breaking the integral in two contributions and performing the substitution q = k/Λ. A factor is absorbed in the cut-offΛ = Λ |L/2 ± x| choosing ± depending on the argument of the sign in each piece.
Since the limit of interest is where x is far away from L/2, this rescaling is well defined and the limitΛ → +∞ still makes sense. The integral reads: Each of the "sign" terms comes respectively from sin (k (L/2 + x)) and sin (k (L/2 − x)) factors in order to take care of the correct sign. One immediately sees that if the signs of L/2 ± x are different, as they are in one of the relevant cases |x| L/2, α 3 is zero. For the other case of |x| L/2 the integral bears similarities to (19) and part of the approach can be followed. Namely, one can divide the integral in small pieces +∞ 0 = n 2π(n+1)/Λ 2πn/Λ and observe that several of them converge to zero as the limit of large cut-off is taken. This is due to the rapid oscillation of the sine (or cosine as in Eq. (19)]. Then the only remaining part is: , |x| < L/2.
(38) In (19) the part analogous to this last piece was also zero as long as e was finite. Now this is no longer correct due to the 1/q factor which picks a large contribution near q = 0. To see this explicitly, one takes the leading order of log 1 + q −2 / 1 + e 2 4π 2 log 1 + q −2 for small q which is simply 4π 2 e 2 . Then the result is independent of the cut-offΛ: Summing up then the results: which again corresponds to the expected behaviour for a confined phase. The string tension is given by: which is finite in general. It is interesting to observe that even though σ 1 = σ 2 , nevertheless they are still different from σ 3 . Furthermore, in the two previous cases if the external charge Q is a multiple integer of e the string tension vanishes. For D = 3 it is no longer the case. The string tension remains finite when Q is a multiple integer of e with the factor of 2π in the argument of cosine replaced in (41) by 4Si (2π) 5.67 < 2π. Total screening is obtained instead for Q = π 2Si(2π) e.

VI. CONCLUSIONS
In this work, the robustness of the confined phase for 1 + 1 fermions was studied by determining the string tension between two probe static charges. In the first part of the paper, the robustness of the confinement properties when a Thirring interaction term is added to the Schwinger model was investigates, as well as the effect of a topological θ-term. Through bosonization, it was shown that the known results for both models (i.e., when different types of interactions are separately considered) still hold. The theory only makes sense when the Thirring coupling is g > −π (as in the Thirring model) and most importantly, the system only deconfines for θ = ±π (as in the Schwinger model). Through an Hubbard-Stratonovich transformation, one observes that this model can be regarded as a fermionic field interacting with a massless gauge field which in turn interacts with a 'massive gauge field.
It is possible that general interactions of this form may break confinement. This can be an interesting problem to address in the future since, in principle, such interactions may be experimentally available in the context of quantum simulations. However, the terms obtained from the Thirring interaction induce cancellations which would not appear under a general coupling between the gauge fields. The Thirring parameter does not allow to vary interaction between the bosonic fields but only the mass of one of them. This result shows that, with respect to confinement, a possible nearest-neighbour densitydensity interaction plays no role and therefore the phase is stable. In the case considered here the static charges coupled directly to the gauge fields (through −iQA µ j µ ext ), but not directly to the fermion current, even in the presence of the Thirring term. Such term would take the form 1 2 g extψ γ µ ψj µ ext and, therefore, has the appearance of an external field. The study of such systems is beyond the scope of the present work, but it may also be theoretically and experimentally relevant. Estimates of the parameters for a possible experimental implementation with ultracold atoms in optical lattices have been also presented. The possibility of higher dimension of the gauge fields while the fermions remain in 1 + 1 dimensions was also considered. It was found that in the massless case, as for the Schwinger model, there is a strong screening when static charges are introduced on the system. In the Schwinger case the linear growth of the energy with the distance is replaced by an exponential decay. In the case of a 2 + 1 dimensional gauge field the logarithm is replaced by oscillatory functions (which go to zero as power laws). Finally in 3 + 1 dimensions the 1/L decay is replaced by zero: external charges are completely screened. When a small mass is considered, a linear growth of the energy with the distance is observed and therefore a finite string tension is obtained for all 1 + 1, 2 + 1 or 3 + 1 dimensions. Furthermore, at this order in perturbation theory, the string tension is the same for the first two cases and smaller for the latter one: σ 1 = σ 2 > σ 3 . The last inequality only holds for small enough external charges (notice anyway that the string tensions are periodic functions of the external charges).
This result is somewhat counter-intuitive since the confinement in the Schwinger model is usually attributed to the fact that the Gauss law in 1+1 dimensions impose a constant electric field (rather than 1/r 2 of the 3 + 1 system). Our results suggest that this feature is not necessary to obtain confinement and, instead, it is the dimensionality of the space-time available to the fermion fields that is rather dictating confinement. In order to better test this hypothesis it would be interesting to study how far can one extend the space-time allowed for the fermion fields before leaving the confined phase (for gauge fields in 3 + 1 dimensions for example). Given that it is known that when the fermion fields span the full 3 + 1 dimensions the theory is deconfined (corresponding to regular QED), this point of transition does exist. Furthermore, with the advent of quantum simulation of gauge theories, one can hope that an experiment with tunable fermion dimensionality [64] could probe directly interesting phenomena like this transition.
We finally observe that the case of the gauge fields living in 3 + 1 dimensions exhibits a quantitative difference with respect to the other two in which the gauge fields are defined in 1 + 1 and 2 + 1. It would be interesting to understand in detail how the models differ. A quantity of interest would be the expected fermion distribution in the presence of external charges. In the 3 + 1 case the matter-gauge system creates, dynamically, a linear growth of the energy from a static 1/r 2 energy interaction. T he way this happens is expected to be quantitatively different from the case of 1 + 1 dimensions, where the static energy interaction is already linear.
if even or sinh if odd. By doing the Fourier transform of (9), exponentiating it and then Fourier transforming again, one can compute the Fourier transform of (A1) in terms of the momentum p. The result is given by: The integration of one of the variables, say q n , can be performed using the Dirac delta. The n − 1 integrations of the zero-th component of q i can be then carried out putting them on-shell. This results in: where it was abbreviated E i = q 2 i + m 2 . Using the notation Q = q 1 + . . . + q n−1 , one can write the denominator of (A3) in the form λ (p − q) 2 + m 2 . The momenta part can be where it was used that it is possible to eliminate the dependence on the spatial component of p by a suitable translation of the spatial variable of integration. The poles then obey Because the particles q i are on-shell, the maximum value of the total for momenta is λQ 2 = − (n − 1) 2 m 2 corresponding to the situation where all the n − 1 particles are at rest in a given frame and therefore Q 1 = 0. For this case one finds a pole at λp 2 0 = −n 2 m 2 . By increasing the total momentum of Q one finds a branch cut along the axis starting at −n 2 m 2 corresponding to multiparticle states. This holds for any n > 1. For the special case n = 1 there is an isolated pole at λp 2 0 = −m 2 and therefore the theory does not contain further states. may depend on the direction and that are fixed by energy minimization. For the present case the shape of the function in the directions y and z is fixed through a parameter σ ⊥α while the value of the longitudinal component, which is called simply σ α , can be fixed variationally. The relevant potentials, that can be realized experimentally, take a form that can be written as: To simplify the subsequent computation a quartic polynomial potential is considered. It is constructed by expanding the expression above in powers of kx and having two minima, one at −a/2 and other at a/2, with an offset between them of ∆V 0 . In order to derive the potential one can require that its first derivative is of the form (A/a) (2x/a + 1) (x/a − x 0 /a) (2x/a − 1) where x 0 corresponds to the position of the maxima between the two minima. The potential depth will be proportional to A. By integrating one obtains the form of the potential and an extra parameter c as a constant of integration. This parameter is fixed by requiring that the absolute minima, chosen arbitrarily to be the one at x = −a/2, corresponds to zero energy. The potential depth, V 0 , is the height of the barrier at x 0 (position of the maximum). With this definition one can replace A by V 0 according to 3 . The potential considered for, say, the boson species 1 is then: while for the boson species 2 is V B2 (x) = V B1 (−x). For the fermions the potential V F (x) has the same structure with an amplitude A F 0 : It is important to observe, however, that only bosons -unlike the fermions -should feel a double well potential according to the proposal. For this reason the polynomial double well potential approximation for the fermions is not as good as an approximation as it is for the bosons. Nonetheless, as one is interested in the strong coupling regime of the model, this provide a reasonable approximation, as we verified. The difference between the energies of the minima, V 0 ∆, should be small when compared to V 0 or, equivalently, x 0 /a should be small. Consequently, their influence on the parameter determination is small, which was checked explicitly. In what follows it will then be taken x 0 /a = ∆ = 0 avoiding unnecessary complicated formulas. This results in and For the considered scheme it is required that U = U 11 = U 22 and U 12 = U 1F = U 2F = 2U . The fine-tuning of this condition is not crucial as discussed in the main text. With the parameters σ ⊥ and σ fixed, one has to rely on the control of the scattering length in order to fulfil this condition. Within the variational approach, one computes the average energy per site and requires that σ minimizes it. The problem of the different shape of the minima, also referred in the main text, can be addressed as follows. The total energy is given by: ε = d 3 r r ,α n α 2 2m α |∇φ α r | 2 + n α V ext |φ α r | 2 + + β>α n α n β g αβ 2 |φ α r | 2 |φ β r | 2 (C8) (n α is the number of atoms per well of the species α). Due to the asymmetry, the total energy per site is different depending on which minima one is referring V ext |φ α r | 2 . For the minima at x = ±a/2 the result is: where it was already included the approximation that the spreading in the perpendicular direction is the same for all species and characterized by σ ⊥ . Assuming that all masses are the same and the offsets ∆ α = 0, the only parameters species-dependent are the densities n α and the amplitudes V 0α . The asymmetry of the minima is present whenever ∆ α = 0. Therefore the problem of the asymmetry of the different Wannier function is not present here. It was checked that the obtained estimates do not depend very much on the ∆ parameters if they are not too large. The density for fermions is n F = 1 while for the other two species of bosons is n 1,2 = S. With this there are in total two parameters to fix, σ B and σ F , given two coupled equations ∂ε/∂σ F = 0 and ∂ε/∂σ B = 0. The reference energy is denoted by E ref = 2 /2ma 2 with an assumed equal mass for all the species m. The following dimensionless parameters are now introduced:Ṽ 0α = V 0α /E ref ,σ α = σ α /a and a αβ = a αβ /a. Regarding the scattering lengths, one is working on the assumption that a 1F = a 2F = a 12 ≡ 2a scatt and a 11 = a 22 = a scatt . The two equations are then: As discussed in the main text and above one takes a 12 = 2a scat and a 11 = a 22 = a scat . Due to the difference on the interaction terms of the two equations of (C10) the assumption a 1F = a 2F = 2a scat does not automatically satisfy the requirement of the Hamiltonian parameters of the proposal. For S small, at least, the result is approximately valid so these values are also taken as reference for the scattering between bosons and fermions. Equations (C10) are used to obtain the data reported in Figure 2.
An important check concerns whether the values of the parameters validate the perturbative approximation obtained for large values of U . This amounts to check that t α /U and V 0α ∆ α are actually small (here it is considered that they should be ∼ 0.1 or smaller). For illustrative purposes it was fixed S = 1, V 0F = V 0B = V 0 and ∆ F = ∆ B = ∆. Direct analysis of the above equations shows that in order to guarantee that the perturbative regime is valid in the interval V 0 ∼ 3 − 10, then one should haveσ ⊥ ∼ 0.2 and ∆ 10 −3 . If one takesσ ⊥ to be, say, two or three times higher than this, larger potential amplitudes are required. Alternatively, larger scattering lengths could also be used. From the other side, there is some freedom in choosing the values of ∆ in order to remain in the perturbative regime. However this choice should respect the fact the two minima should still be present at x = ±a/2 which is translated into |∆| < 1/3. Finally, the mass parameter of the target model will scale as V 0 ∆ and the choice was taken such that the energy scale of this term matches the order of magnitude of the other terms in the Hamiltonian t α t β /U ∼ V 0 ∆, which results to be ∆ 10 −3 , as referred in the main text.
Appendix D: Equations of motion and energy-momentum tensor for theories with higher derivatives Here a classical field theory with higher derivatives is considered. The well known Euler-Lagrange equations are derived by the extremization of the action. The inclusion of higher derivatives on the Lagrangian lead to a reformulation of the equations. In fact by calculating explicitly δS = 0, integrating by parts whenever necessary one obtains: where N is the highest number of derivatives appearing in a term of the Lagrangian. For N = 1 one recovers the usual Euler-Lagrange equations. Consider now a general translation x µ → x µ + ε µ . The total change of the Lagrangian is which results in The derivatives can be written as acting on ∂L ∂(∂µ 1 ...∂µ n φ) with a minus sign plus a total derivative term: By continuing this process with every ∂ µi acting on φ in the terms that are not an exact derivative (last term), one obtains: where, in order to simplify the notation, it was written ∂L ∂ ∂φ ≡ ∂L ∂ (∂ µ1 . . . ∂ µn φ) .
The special case of n = 0 gives ∂L ∂φ ∂ ν φ. Summing over all n, joining all terms coming from the part of the above expression one can recognize the equations of motion (D1) and therefore they are put to zero. What remains is a series of total derivatives. Furthermore in itself the Lagrange density changes as δL = ε ν ∂ ν φ. By rearranging the variables one obtains: and therefore one can identify the energy-momentum tensor. This is just the conserved current that follows from Noether's theorem for the special case of space-time translations: