Implementing quantum electrodynamics with ultracold atomic systems

We discuss the experimental engineering of model systems for the description of QED in one spatial dimension via a mixture of bosonic $^{23}$Na and fermionic $^6$Li atoms. The local gauge symmetry is realized in an optical superlattice, using heteronuclear boson-fermion spin-changing interactions which preserve the total spin in every local collision. We consider a large number of bosons residing in the coherent state of a Bose-Einstein condensate on each link between the fermion lattice sites, such that the behavior of lattice QED in the continuum limit can be recovered. The discussion about the range of possible experimental parameters builds, in particular, upon experiences with related setups of fermions interacting with coherent samples of bosonic atoms. We determine the atomic system's parameters required for the description of fundamental QED processes, such as Schwinger pair production and string breaking. This is achieved by benchmark calculations of the atomic system and of QED itself using functional integral techniques. Our results demonstrate that the dynamics of one-dimensional QED may be realized with ultracold atoms using state-of-the-art experimental resources. The experimental setup proposed may provide a unique access to longstanding open questions for which classical computational methods are no longer applicable.


I. INTRODUCTION
The experimental engineering of atomic model systems for the description of dynamical gauge fields represents a major challenge with most important applications. Fundamental gauge fields mediate the strong and electroweak forces between matter in the standard model of particle physics, where the photon of quantum electrodynamics (QED) is a most prominent representative [1]. Gauge fields can also arise as emerging degrees of freedom in strongly correlated condensed matter systems such as related to the quantum Hall effect [2] or effective theories of spin liquids [3]. Most pressing questions concern the realtime evolution of strongly interacting gauge fields coupled to fermionic matter (QCD), such as realized during the early stages of our universe and explored in collisions of ultrarelativistic nuclei at giant laboratory facilities [4].
The complex many-body dynamics of gauge fields is often very difficult to study in the original systems both experimentally and theoretically. For instance, in a heavyion collision most experimental observables give information only about the integrated space-time evolution of the system. Its theoretical description is complicated by the fact that ab initio computer simulations of the realtime dynamics can only be achieved in certain limiting cases because of the so-called sign-problem [5]. Here, the experimental engineering of atomic quantum simulators appears as an attractive alternative, which may provide a unique access to longstanding open questions [6][7][8].
Compact table-top experiments with ultracold atoms * kasper@physics.harvard.edu are rather easily accessible and provide a very flexible testbed, with tunable interactions or reduced dimensionality by shaping the confining optical potential [9,10]. Since setups employing ultracold quantum gases can be largely isolated from the environment, they offer the possibility to study fundamental aspects such as the unitary real-time evolution of systems with engineered Hamiltonians to high accuracy [11]. The realization of external static gauge fields has been achieved in many impressive experiments with ultracold atoms, such as in Refs. [12,13]. In contrast, no experimental implementation of gauge fields as dynamical degrees of freedom, as in QED or QCD, in an atomic setup has been achieved yet, although it has been theoretically discussed in [6][7][8]. The presence of dynamical gauge fields in the description of a physical system reflects an underlying local symmetry, whose space-time dependent transformation properties significantly constrains the quantum dynamics allowed. Implementing such a gauge symmetry for a system of bosonic and fermionic atoms can lead to involved constructions, which often rely on higherorder processes that are challenging to control experimentally.
Further experimental progress can be facilitated with the identification and implementation of simple gauge field theory examples. Abelian gauge symmetries, such as the local U (1) symmetry of QED, clearly stand out in this respect, in particular, if they are implemented in one spatial dimension. Abelian gauge theories in the continuum are simpler than non-Abelian theories such as QCD because of the absence of self-interactions of the gauge bosons. In QED the photon interacts only via processes involving electrons and positrons, while the gluons in QCD directly interact with each other. In a lattice-discretized theory in more than one spatial dimension even QED magnetic fields appear as ring exchange interactions which makes a possible experimental implementation via effective interactions [14,15], higher-order perturbative processes [16] or ancillary degrees of freedom [17,18] involved. In one spatial dimension, however, no QED magnetic field exists. This dramatically simplifies possible descriptions of the interaction of the remaining electric field with the fermions. In this case, the interaction terms respecting the local gauge symmetry can be realized using heteronuclear boson-fermion spin-changing interactions which preserve the total spin in every local collision [16,19]. While a Jordan-Wigner transformation can be used to express the fermions as quantum spins, our construction does not rely on this mapping but directly simulates the fermionic degrees of freedom. This is particularly important in view of experimental implementations of gauge theories in dimensions larger than one, for which the Jordan-Wigner mapping is less useful. Despite the reduced complexity, the Abelian gauge theory setup in (1 + 1) space-time dimensions still offers a rich phenomenology, including important dynamical strong-field phenomena such as Schwinger pair production [19,20] and string breaking [21], which are highly relevant for many systems also in more than one spatial dimension.
Very interesting and detailed suggestions have been made to realize gauge field dynamics in atomic systems, where many proposals concentrate on quantum link models [6]. In these models the gauge fields are regularized using quantum link variables which have a finitedimensional link Hilbert space, and whose dimension is determined by the number of bosonic atoms residing on a given link between fermionic atoms in an optical lattice. Since the Hilbert space of a quantum link model is finite, the mapping to atomic systems is in general facilitated. Many ground-breaking investigations have been performed using a small number of bosons per link [8,15,16,21,22]. A low-dimensional Hilbert space also allows one to achieve theoretical estimates based on diagonalization or tensor network techniques [23,24]. Since the Hilbert space of QED or QCD itself is infinitedimensional, the Hamiltonian formulation of the original gauge field theory on a spatial lattice [25] can be recovered for a sufficiently large number of bosons 1 .
In this work we follow Ref. [19] and consider a mixture of bosonic 23 Na and fermionic 6 Li atoms in a onedimensional optical superlattice. We concentrate on the 1 The mapping becomes more involved if only a small number of bosons per link is employed. In this case, the quantum fields of the original gauge theory can arise as low-energy effective degrees of freedom of the theory of quantum links after dimensional reduction. More precisely, the quantum fields of a gauge theory in D dimensions are obtained from a (D + 1)-dimensional theory of quantum links. To recover one-dimensional QED or QCD would, therefore, require a two-dimensional quantum link setup where the extra dimension could also be implemented with the help of internal degrees of freedom [6].
regime with a large number of bosons residing on each link, such that the results of the original lattice gauge theory in the continuum limit are recovered. Our discussion about the range of possible atomic system's parameters builds, in particular, upon experiences with related experimental setups of fermions interacting with coherent samples of bosonic atoms [26][27][28][29][30][31]. The other important ingredient of our investigation concerns benchmark calculations of the atomic system and of the original gauge theory. Since exact diagonalization techniques are no longer applicable in this case, we employ powerful functional integral (FI) techniques [19,[32][33][34][35]. They allow us to do ab initio calculations in an important range of strong-field phenomena. Reproducing these benchmark results with future experimental realizations will be a crucial milestone, before new regimes can be explored that are no longer accessible with classical computational methods. This publication is organized as follows. In section II we briefly review QED in 1+1 space-time dimensions, i.e. the massive Schwinger model. We employ a lattice discretization with staggered fermions to connect the gauge theory to an atomic model in an optical superlattice with angular momentum conserving scattering processes. In section III we discuss in detail a possible experimental implementation of the Schwinger model in a mixture of bosonic and fermionic atoms, where gauge invariance requires a correlated hopping of the staggered fermions with the Schwinger bosons residing on the links. While the basic discussion follows to a large extent Ref. [16], we concentrate on the implementation in specific systems with given experimental parameters to ensure the relevant separation of scales that is required to suppress contributions from gauge-symmetry violating states. Moreover, we employ species-dependent lattices to separate the bosonic and fermionic degrees of freedom in order to simplify the experimental realization. In that section we also translate the basic quantities of the cold atom system to the fundamental parameters of the corresponding lattice gauge theory. A set of viable parameters to be employed in an upcoming experiment is presented in section IV. The later sections are devoted to benchmark calculations in order to demonstrate that relevant QED processes can indeed be described using the available experimental parameters of the atomic setup. In section V we review the functional integral approach and derive equations of motion for the cold atom system. This method allows us to accurately describe the nonequilibrium dynamics of coherent bosonic fields coupled to fermions from first principles. In section VI we study the real-time dynamics of Schwinger pair production and string breaking in the cold atom system. This section is an extension of the pair-production results of Ref. [19] to the new parameter sets established in this work, and to the phenomenon of string breaking that has not been considered in the large boson number regime of the atomic setup before. We present the dynamics of various experimentally accessible observables and discuss the accuracy with which QED can be represented in practice by a finite atomic system. We conclude and give an outlook in section VII.

II. THE SCHWINGER MODEL REVISITED
Quantum electrodynamics for massless fermions in one spatial dimension (Schwinger model) is an exactly solvable field theory [36]. On the other hand, no analytic solution is known for massive fermions (massive Schwinger model) [37]. From a phenomenological point of view, a particular interest in this model stems from the fact that it shares several characteristic aspects with the theory of strong interactions (QCD) such as spontaneous chiral symmetry breaking or dynamical string breaking (see e.g. Refs. [38][39][40][41]).
Nonperturbative studies of the massive Schwinger model are typically based on a lattice discretization of the continuum theory. The construction of a hermitean, local and translation-invariant lattice theory of fermions necessarily entails the appearence of unphysical degrees of freedom, the so-called fermion doublers [42]. One possibility to resolve this problem is by making the spurious doubler modes heavy via a Wilson term [43] and we refer to Refs. [32,33,44] for numerical studies using this approach. In this work, we employ the alternative staggered fermion discretization where the Dirac spinor is decomposed in such a way that the doubler modes can be disregarded as they decouple. As a consequence, the particle and antiparticle components reside on neighbouring lattice sites [25]. The Hamiltonian of the theory is given by where a S denotes the lattice spacing and M the fermion mass [25,36] and we choose to work in natural units ( = c = 1). Here, the staggered fermion field operator ψ n , which resides on lattice sites n, fulfills the canonical anticommutation relation {ψ n , ψ † m } = δ nm . The fermionic charge operator is defined as Accordingly, the presence of a fermion at an even site is interpreted as particle (Q n = +1) whereas the absence of a fermion at an odd site is interpreted as antiparticle (Q n = −1). The unitary link operator U n and the electric field operator E n act between neighboring lattice sites n and n + 1 and obey the commutation relations where g denotes the gauge coupling. This algebra entails an infinite-dimensional local Hilbert space. For the U (1) gauge theory, the Gauss's law operator generates local gauge transformations and commutes with the Hamiltonian [H KS , G n ] = 0. Quantum link models have been proposed as an alternative formulation of gauge theories for which finitedimensional representations of the link algebra exist [45][46][47]. Recently, the prospect of constructing quantum simulators for gauge theories has boosted the interest in these models as their implementation in atomic systems could be greatly facilitated. In this approach, the electric field operator is identified with the z-component of the quantum spin operator L z,n , whereas the link operators are regarded as raising and lowering operators with L ±,n = L x,n ± iL y,n . The quantum spin operators fulfill the angular momentum algebra [L i,n , L j,m ] = iδ nm ijk L k,n and denotes the spin magnitude. Consequently, the commutation relation (3b) is identically fulfilled, whereas the commutation relation (3a) is no longer valid if is kept finite, but replaced by [U n , U † m ] = 2δ nm E m /[g ( + 1)] which goes to zero only as → ∞. Only in this limit the unitarity of the link operator U n is restored again. However, it is central for the whole construction that a finite does not affect gauge invariance as generated by the Gauss's law operator with the Hamiltonian of the quantum link Schwinger model The finite-dimensional representation of the angular momentum algebra makes its implementation in systems of ultracold atoms feasible. For representations with small , numerical methods based on diagonalization or tensor network states provide valuable information about static and dynamic properties [23,24,[48][49][50]. Of course, it is an important question to understand the connection between the finite-dimensional representation of cold atom gauge theories and the infinite-dimensional representation corresponding to QED. In Ref. [19] the large-regime and the convergence to → ∞ results, i.e. QED itself, has been established using powerful functional integral (FI) techniques for strong-field phenomena. In section V we describe the FI approach and apply it to obtain benchmark results for the Hamiltonian (8) using parameter sets motivated by possible experimental realizations.

III. EXPERIMENTAL REALIZATION
Our starting point for the realization of the U (1) gauge theory coupled to fermionic matter in an ultracold atom experiment is a genuine interacting gas of fermionic and bosonic atoms [7]. To facilitate the connection with the experiments we reintroduce where appropriate. The bosons φ α (x) and fermions ψ α (x) fulfill the canonical commutator and anti-commutator relations, respectively, Here, the greek labels α, β denote magnetic hyperfine states of the atoms. The particles are confined by external potentials and interact via inter-and intra-species scattering processes. The corresponding Hamiltonian consists of three parts, H = H T + H V + H I . The kinetic part, H T , describes the movement of the atoms, with masses M b and M f . The potential energy contribution, H V , is determined by the external potentials according to whereas the atomic scattering processes are described by The coupling constants are determined by the scattering lengths of the inter-and intra-species scattering processes. Throughout this work, we indicate purely bosonic terms by the superscript b, fermionic terms by the superscript f and boson-fermion interactions by the superscript bf .
In the following, we describe in detail all the steps that are required so that the gas of fermionic and bosonic atoms in three spatial dimensions (10) -(12) behaves according to the Hamiltonian (8) in one spatial dimension. To this end, we show how to reduce the dimensionality from three to one dimensions and how to construct the staggered lattice for fermions. Afterwards, we describe how to select only those interactions from (12) that correspond to the gauge-invariant interactions in (8).

A. One-dimensional staggered geometry
The basic ingredient for realizing a one-dimensional lattice structure with lattice constant a is an optical lattice with tight radial confinement. Employing a laser frequency which is blue detuned for fermions and red detuned for bosons allows us to place a mesoscopic bosonic gas on the links between the fermionic lattice sites, as indicated in the left graph of Fig. 2. In fact, the potential energy contributions (11) can be split into an axial and radial part, Here, l s,α denotes the radial confinement length scale, where we distinguish bosons and fermions by the superscript s ∈ {b, f }. Owing to tight radial confinement, the three-dimensional system is effectively rendered onedimensional and we employ the product form where ϕ s (y) and ϕ s (z) are the ground state wave functions in the y and z directions, respectively. We assume that these states are independent of the magnetic quantum number.
To generate a staggered structure for the fermions, the original optical lattice with period a is superimposed by an optical superlattice with period 2a, as indicated in Realization of the staggered lattice structure. The black arrows indicate the dipole lasers which are used in the experimental setup.
(1) Blue/red detuning for fermions/bosons generates phase-shifted optical potentials for bosons and fermions with a lattice period a. (2) The superposition of the lattice with period a by a superlattice of period 2a generates the staggered structure for the fermions.
the right graph of Fig. 2. Owing to the fact that the frequency of the second laser is tuned closer to resonance with respect to the fermions than to the bosons, the second lattice does practically not affect the bosonic degrees of freedom. Disregarding the effect of overall confinement in the axial direction, the axial part of the potential is then given by The amplitudes V s i,α , i = 1, 2 are determined by the AC Stark shift of the corresponding magnetic substates α. The tunneling of bosons between adjacent sites is suppressed by tuning of the laser amplitude. As already noted in Ref. [16], such a construction is spin-independent and, hence, distinct from e.g. Ref. [21].
In the following, it is useful to switch to a representation in terms of localized Wannier functions. To this end, we first consider the bosonic degrees of freedom and focus on the two lowest energy bands. By tuning of the laser amplitude, we may choose V f i,α such that two Wannier functions w f α,n p (x) are obtained which are sufficiently localized in the left (p = L) and right (p = R) minimum of the elementary cell n ∈ {0, . . . , N − 1} with positive integer N of the optical lattice, respectively. The corresponding expansion of the fermionic field operator reads We note that the total number of bosonic/fermionic lattice sites is 2N and we will label them by n ∈ {0, . . . , 2N − 1}. Similarly, we may expand the bosonic field operators, φ α (x) = n ,p w b α,n p (x)φ α,n p , where the Wannier functions w b α,n p (x) are again localized in the two minima of the elementary cell. In fact, the structure of the superlattice suggests the definition We note that the kinetic energy contributions (10) are suppressed since the Wannier functions that correspond to the different minima in the optical lattice do not have a sizable overlap (cf. also section IV).

B. Angular momentum conservation
In the previous section, we reviewed how the potential energy (11) can be used to generate the staggered lattice structure. Moreover, is was noted that the kinetic energy (10) is suppressed due to the localization of the Wannier functions at the potential minima. To ensure local gauge invariance and create dynamics in the cold atom system, we have to tune the interaction Hamiltonian (12) such that only a selection of terms contributes. In the following, we discuss all interaction terms in more detail and explain the connection between the various scattering lengths and coupling constants g s αβγδ for s ∈ {b, f, bf }. For pedagogical reasons, we discuss the construction in free space first and take into account the lattice later.
We suppose that the inter-and intra-species interactions of bosons and fermions are local and conserve angular momentum. Specifically, we consider bosonic degrees of freedom with spin f b = 1 and fermionic degrees of freedom with spin f f = 1/2. Therefore, the two-particle potentials are given by where the total spin can take the values F b ∈ {0, 2}, F f ∈ {0, 1} and F bf ∈ {1/2, 3/2} [51]. The interaction strengths g s,Fs are related to the s-wave scattering lengths a s,Fs via Here M r,s denotes the reduced mass of the two scattering partners. In general, the projector P F for two particles with individual spins f 1 and f 2 on the subspace with total spin F can be written as where M ∈ {−F, −F + 1, . . . , F − 1, F} are the possible magnetic quantum numbers. We may relate the interaction strengths g s,Fs to the constants appearing in the interaction part of the Hamiltonian (12) according to Here, f 1 ; α; f 2 ; β|f 1 , f 2 ; F s , M are the Clebsch-Gordan coefficients for coupling the individual spins f 1 and f 2 to the total spin F s . Specifically, we have Following along the lines of the previous section, we reduce the three-dimensional system to a setup with effectively one spatial dimension and expand the field operators in terms of Wannier functions. Using a compact notation, where n = (n 1 , n 2 , n 3 , n 4 ) denotes the site indices and µ = (α, β, γ, δ) the magnetic quantum numbers, we can write for the interaction part of the Hamiltonian: The coupling constants U s n are determined by the dimensional reduction and the explicit form of the overlap integrals of the Wannier functions, as given in Appendix A.
Based on this interaction Hamiltonian, we see that a plethora of possible interaction terms are generated in general. In order to realize the Hamiltonian (8), however, we have to guarantee that only specific terms contribute that respect the gauge symmetry. To this end, we use the fact that the application of an appropriate magnetic field and radio-frequency (rf) dressing allows for a selection of a small number of relevant interaction terms, whereas all other contributions become suppressed. We emphasize that this selection is achieved by the unequal shift of the bosonic and fermionic energy levels, as depicted in Fig. 3. Most notably, this procedure results in the bosonic spin exchange with the simultaneous fermion hopping that corresponds to the gauge-invariant interaction term in (8). We note that this selection process does not exclude elastic scattering terms, i.e. scattering processes without changing the individual spins of the atoms.
All bosonic states are prepared in α b = {−1, 0} states whereas the fermionic degrees of freedom are generated in the staggered configuration with α f = 1/2 on even sites and α f = −1/2 on odd sites. As a consequence, interactions including the α b = 1 sector, which are allowed in principle, are suppressed at all times if initialized accordingly. We further elaborate on this issue in the following sections.

C. Bosonic intra-species interactions
In this section, we discuss the intra-species interaction terms of bosons in more detail. Owing to localization in the optical lattice, only on-site interactions of bosons contribute, i.e. U b n = 0 for n = (n, n, n, n) and all others effectively vanish. Accordingly, the relevant part of the purely bosonic term in the interaction Hamiltonian (22) is given by where all entries of µ = (α, β, γ, δ) may take values α b ∈ {−1, 0}. Again, we note that we disregard terms including the magnetic substates α b = 1 which are excluded by the spin conservation if initialized accordingly. The interaction term (23) reads where the coupling constants result from (21) and we assumed that the overlap integrals U b n are the same for all terms. For later convenience, we denote the bosonic degrees of freedom on even sites as φ 2n,−1 ≡ d 2n and φ 2n,0 ≡ b 2n , whereas we interchange their role on odd sites such that φ 2n+1,−1 ≡ b 2n+1 and φ 2n+1,0 ≡ d 2n+1 . In fact, the bosons b n and d n can be understood as Schwinger bosons [52] with the identification which constitute a representation of the angular momentum algebra [L i,n , L j,m ] = iδ nm ijk L k,n with L ±,n = L x,n ± iL y,n . The constraint 2 = b † n b n + d † n d n is fulfilled because the hopping of bosons between neighboring sites n → n ± 1 is suppressed. The bosonic intra-species interaction Hamiltonian is then given by where we disregarded an irrelevant constant and introduced the abbreviation ∆ b,0 ≡ (2 − 1)

D. Fermionic intra-species interaction term
In this section, we discuss the intra-species interaction terms of fermions in more detail. Again, only on-site interaction terms contribute owing to localization such that U f n = 0 only for n = (n, n, n, n). Taking into account the Clebsch-Gordon coefficients, the purely fermionic term in the interaction Hamiltonian (22) can be reduced according to In general, this four-fermion interaction term influences the dynamics. However, the contribution can be written as a density-density interaction between α f = −1/2 and α f = 1/2 particles with density operators ρ n,±1/2 = ψ † n,±1/2 ψ n,±1/2 . Restricting ourselves to an initial state |Ψ with only α f = 1/2 particles on even sites and only α f = −1/2 particles on odd sites, one immediately finds that H f I |Ψ = 0. Consequently, this four-fermion interaction does not contribute to the time evolution due to an appropriate initial-state preparation.

E. Inter-species interaction term
Regarding the fermion-boson scattering contributions to the Hamiltonian (22), we have to consider both the spin exchange process as well as elastic scattering processes. According to the interaction selection process described above, the spin exchange term that involves the correlated hopping of fermions and bosons is given by with µ = (α, β, γ, δ) = (−1/2, 0, 1/2, −1). The first term corresponds to Fig. 4a whereas the second term is shown in Fig. 4b. According to (21), the coupling constant for this specific scattering process is given by We emphasize that the spin exchange term does not change the staggered occupation of fermions such that the four-fermion term (27) still does not contribute. We anticipate that this applies to the elastic scattering terms as well. Accordingly, the fermions are completely determined by their parity (even/odd sites) and we may therefore drop the spin label completely, such that Employing the Schwinger boson representation and taking into account that the overlap integral U bf n does not depend on the specific lattice site n, the spin exchange Hamiltonian can be written as The elastic scattering processes, on the other hand, are given by where β ∈ {−1, 0}. We note that the coupling constants g bf µ still depend on the magnetic substates and are, therefore, not identical for the different terms. Moreover, we find that each U bf n is independent of n in (33), and identical in the first and second line (further denoted by U bf n1 ) as well as in the third and fourth line (further denoted by U bf n3 ), cf. Appendix A. For the first term in (33) with where we used (21) again. The second term in (33) is the same as the first one upon replacing b 2n → d 2n−1 and d 2n → b 2n−1 . The third term in (33), however, is different owing to µ = (1/2, β, 1/2, β) corresponding to the different fermionic parity and reads The fourth term in (33) is the same as the third one upon replacing b 2n → d 2n+1 and d 2n → b 2n+1 . Employing the Schwinger-boson representation b † n b n = + L z,n and d † n d n = − L z,n , the first term (34) can be written as Similarly, the third term (35) is given by and similar expressions are also obtained for the second and fourth term by replacing L z,2n → −L z,2n∓1 . Accordingly, we obtain a contribution withg bf = (g bf,1/2 − g bf,3/2 )/6 that does not appear in the target Hamiltonian (8). The difference of the L z operators can be expressed in terms of the Gauss's law operator, L z,n −L z,n−1 = G n +ψ † n ψ n +[(−1) n −1]/2. Any term proportional to G n exactly vanishes upon acting on physical states so that we can disregard them. In summary, the elastic scattering terms result in bilinear, parity-dependent fermionic contributions that account to the potential energy of the fermions: Summing up all contributions that originate from the kinetic, potential and interaction terms, the cold-atom QED Hamiltonian takes the form The terms in the last line include the potential energy contributions, in particular the energy due to the trapping of the atoms as well as the bilinear term (39). Again, we emphasize that the fermionic contribution V f n depends on the parity of n. Introducing the en- , the fermionic potential contribution is given by Owing to total particle number conservation, the first term does not contribute to the dynamics and can thus be disregarded. The bosonic potential term is treated in a similar fashion. In fact, defining ∆ b,1 according to where we disregarded an irrelevant constant proportional to V b 0 which only depends on . Adding the second contribution in (40) and defining ∆ b ≡ ∆ b,0 +∆ b,1 , we obtain Comparison of the cold-atom QED Hamiltonian with the target Hamiltonian (8) then shows that we have a term linear in L z,n . However, this term does not contribute. To see this, we transform the Hamiltonian into the interaction picture. To this end, we split the cold-atom QED Hamiltonian into two parts, H = H 0 + H 1 , with Here, we introduced the detuning 2∆ ≡ ∆ b − ∆ f , the effective boson self-interaction χ BB and the boson-fermion interaction χ BF , The index n is the same as in (23) and (29) and the overlaps U n are determined in (A11) of the appendix. Upon acting with the unitary transformation Performing the canonical transformation ψ n → (−i) n ψ n , we can finally identify H 1 with the quantum link Hamiltonian (8) where we introduced the abbreviations and time is measured in units of M instead of ∆. We note that we have to take the limit → ∞ in order to recover the Hamiltonian formulation of lattice QED corresponding to (1).

IV. MICROSCOPIC PARAMETERS
At this point, we are now able to determine the accessible parameters for an experimental implementation of the Schwinger model via a mixture of bosonic 23 Na and fermionic 6 Li atoms [30], which is determined by the parameters χ BB , χ BF , ∆ and the occupation numbers of the links.
The wave functions that appear in the overlap integrals for χ BF and χ BB are determined by the lattice spacing a, the lattice depths V b , V f 1 and V f 2 appearing in (15), the difference of scattering lengths that drives the spin-changing collisions as well as the atomic density on each site. For Na + Na in the f b = 1 manifold, one obtains a b,0 − a b,2 = 5 a 0 [53], where a 0 is the Bohr radius. For the Na + Li collisions, we have a bf,3/2 − a bf,1/2 = 0.9 a 0 [54]. Further, we choose a = 4 µm, Here the number of lattice sites can range from a few sites up to N 100, where the limit arises typically due to different experimental constraints. A judicious choice of the orthogonal confinement allows us to use the same wave function for bosons and fermions in the orthogonal direction ϕ s (y) = ϕ s (z) with ω ⊥ /2π = 5 kHz. This potential results in χ BF /h = 0.05 Hz, χ BB /h = 0.58 Hz. We choose the occupation of the bosonic links to be N B ≈ 100 atoms per well, such that the estimated time scale for three-body loss is of the order of several seconds. This means that the Bose-Fermi coupling will be typically of the order of χ BF N B /h ≈ 5 Hz. It will lead to a hopping of fermions between neighboring lattice sites if the detuning is not too large and we therefore choose a detuning of ∆/h ≈ 10 Hz. Further, the chosen experimental parameters ensure that the bosonic and fermionic tunneling J B /h ≈ 0.25 Hz, J F /h ≈ 0.25 Hz (see appendix (A12)) is ten times slower than the expected gauge field dynamics. So we can indeed neglect the direct tunneling for a description of the cold-atom system as done in the previous derivation.
Given the microscopic parameters of the model, we have in mind applications to strong-field phenomena such as Schwinger pair production or string breaking, for which we aim to perform benchmark simulations of the cold-atom Hamiltonian (40). To study Schwinger pair production, the electric field E is supposed to exceed the critical field strength E c = M 2 /g such that the amplitude E/E c 1. In the atomic system, this is given by The achievable electric field exceeds the critical field strength with E/E c ≈ 1.5 for the given experimental parameters.
This puts the exciting prospect of studying Schwinger pair production using ultracold atoms within reach of current experimental techniques. For a successful quantum simulation of Schwinger pair production, however, it does not suffice to implement only the appropriate Hamiltonian. Additionally, we also need to verify (i) that a proper initial state can be prepared, (ii) that the QED dynamics is indeed accessible to the experimental protocol, and (iii) that we can read-out the relevant quantities experimentally. In fact, all these requirement can be fulfilled with current experimental techniques for a range of important strong-field phenomena, as further described in section VI B.

V. FUNCTIONAL INTEGRAL APPROACH
In this section we outline the functional integral (FI) method to investigate the real-time dynamics of fermions coupled to bosonic fields and for convenience we perform the following theoretical calculations in natural units. Since identical fermions cannot occupy the same state, their quantum nature is highly relevant and a consistent quantum theory of nonequilibrium Schwinger pair production including backreactions is envisaged. Based on the defining functional integral of the quantum theory, a systematic expansion can be achieved where the corrections are given in terms of a small (dimensionless coupling) parameter for strong-field phenomena as reviewed in Ref. [34]. At lowest order in this expansion one recovers a classical-statistical field theory for bosonic fields or the so-called Truncated Wigner approximation. The inclusion of fermions requires to go beyond lowest order, which we describe below and apply to the cold-atom Hamiltonian (46) in section VI. To study the strong-field regime of QED, the field strength needs to be of the order of the critical field E c = M 2 /g. The corresponding cold atom setup is characterized by To make contact with the Truncated Wigner approach, we start from a purely bosonic theory [57]. The expectation value of an observable O is given by the trace where the time-dependent density operator ρ(t) in the Schrödinger picture obeys the von-Neumann equation Its discretized version may be used as a starting point for deriving a functional integral expression of the time evolution for the density matrix. We perform a Wigner transformation of the discretized von-Neumann equation (51) where the transformation only acts on bosonic degrees of freedom in the density matrix and the Hamiltonian. The Wigner transform, which depends on the center field ϕ 1 , is denoted by ( · ) W , and we refer to appendix B for further details. Here, ϕ 1 may carry additional indices that will be suppressed in the following. We emphasize that there are no operators present anymore as they are replaced by c-number variables. Employing the product formula for the Wigner transform (B7), the last equation can be written as The product formula of Wigner transforms introduces an additional integration with respect to the difference field η 1 . The naming 'center field' and 'difference field' is motivated by the Schwinger-Keldysh formalism as reviewed in Ref. [34], in which the fields on the forward and backward branch of the closed-time path correspond to ϕ ± = ϕ ± η/2. Iterating this expression, we obtain the functional integral representation with t N = t 0 + N ∆t. The exponent in the last expression can be written as which is a discretized version of the Schwinger-Keldysh action. Expanding the Hamiltonian up to linear order in the difference field η k , we obtain Within this approximation, the difference field η k can be integrated out as it appears only linearly in the exponent. This gives The arguments of the complex Dirac-delta functions impose the discrete time evolution equation [57] i In fact, this equation can be solved implicitly, such that ϕ k = g(ϕ k−1 ). To actually integrate out ϕ k for k ∈ {0, . . . , N − 1}, we have to change the argument of the delta function, such that A similar derivation can be performed for real scalar field theories [58]. The Jacobi determinant takes the form where we used and being a small number. This shows that the Jacobi determinant does not contribute [59], since it affects the dynamics only at order O(∆t 2 ). The Wigner function ρ W (ϕ N ; t N ) for arbitrary times t N is obtained by sampling over initial conditions ϕ 0 , which are then evolved by (60). Expectation values of bosonic observables O are then given by

A. Interacting boson-fermion theory
We now include the quantum corrections induced by the coupling to the fermions. Here, Wigner transforms are only calculated with respect to the bosonic variables whereas the fermionic operators and their appropriate time-ordering still has to be taken into account. For the sake of simplicity, we denote the fermionic fields by ψ and emphasize that they carry additional indices which will be suppressed in the following. Furthermore, we assume that the Hamiltonian H = H B + H F can be separated into a purely bosonic part H B and a quadratic fermionic part, H F = ψ † h F ψ, where h F is a matrix and may contain bosonic degrees of freedom.
The Wigner transform of the discretized time evolution equation (51) for a single time step is now given by where we emphasize again that the Wigner transform H W (ψ, ϕ) depends on both fermionic operators ψ and c-number variables ϕ. Iterating this expression, we obtain again a functional integral representation of the time evolution of the Wigner function where we introduced the time ordering operator T and the anti-time ordering operatorT . Based on the assumption that H = H B +H F , we may also separate its Wigner transform into a purely bosonic and a fermionic contribution Furthermore, we assume that the initial density matrix factorizes into a purely bosonic part ρ B as well as a purely fermionic contribution ρ F . Accordingly, the Wigner transform at initial times only affects the bosonic part, such that also factorizes. However, the factorization property of the density matrix may be lost during the time evolution. We integrate out the fermions to get the time evolution of bosonic observables O B . Denoting the trace over fermionic operators by Tr F , we obtain Owing to the fact that the purely bosonic part is the same as in (54), we focus on the fermionic contributions in the following. Since H F is taken to be quadratic in the fermionic operators, it is convenient to introduce the abbreviation Here, h F W (ϕ) is a matrix that only depends on bosonic variables but not on fermionic operators. Further, we introduce ϕ ± k as the linear combination of center field ϕ k and difference field η k , and we assume that the initial fermionic density matrix can be written as where Z is an appropriate normalization. By utilizing the identity [60] Tr F e ψ † M1ψ . . . e ψ † Mnψ = det(1 + e M1 . . . e Mn ) (74) with matrices M k for k = 1, . . . , n, we may explicitly perform the fermionic trace in (70). Introducing the evolution matrix where S k,m (ϕ) depends on the string of fields ϕ k , ..., ϕ m for k ≥ m, we obtain for the last two lines of (70). For later convenience, we also define the fermionic propagator whose time evolution at order O(∆t) is governed by The components of D k can be identified with the equaltime correlation function ψ † m ψ n as will be shown below. This discrete equation can be solved via a mode expansion of the operator ψ m as illustrated in Ref. [34]. Knowing the mode functions then allows one to compute fermionic correlation functions.
The expansion of the Tr log in (76) up to linear order in the difference field η k then yields Since the difference field η k appears only linearly in the exponent, we may integrate it out again. As compared to the purely bosonic case, cf. (60), the resulting complex Dirac-delta function now impose the discrete time evolution equation [61] i which is implicitly solved via ϕ k =g(ϕ k−1 ). As in the purely bosonic case, the Jacobi determinant does not contribute at leading order. Accordingly, the Wigner function ρ W (ϕ N ; t N ) is obtained by sampling over initial conditions and subsequent time evolution of ϕ and D according to (80) and (78), respectively. Finally, to calculate the expectation value of bilinear fermionic observables O F = ψ † Aψ with A being a matrix, we consider While the bosonic contributions are identical to (70), the fermionic trace now takes the form we may again perform the fermionic trace explicitly. Expanding the exponent up to linear order in the response field η k while setting it to zero otherwise, we finally obtain Accordingly, the discrete time evolution of O F is determined by By explicitly introducing spatial indices and choosing A mn = δ mm1 δ nn1 , we recover the evolution equation of D k as given in (78). This shows that D k can indeed be identified with the equal-time correlation function ψ † m ψ n .

B. Equations of motion for the cold atom system
To derive the equations of motion for the cold atom system, we apply the method from the previous section to the Hamiltonian (46). To this end, we rewrite H QL in terms of the Schwinger boson operators b n and d n . Upon employing symmetric operator ordering and skipping an irrelevant constant, the Wigner transform of the Hamiltonian is given by [62] where b n and d n are now c-numbers, but the ψ n and ψ † n are still fermionic operators. In the following, time is treated as a continuous variable. Accordingly, all equations correspond to their discrete versions up to order O(∆t 2 ). The equations of motion for the fermionic correlator D(t) are given by where the fermionic matrix h F W (t) reads The fermionic two-point function D mn = ψ † m ψ n = 1 2 (δ mn − F nm ) can be expressed in terms of the statistical two-point function F mn = [ψ m , ψ † n ] . Accordingly, the equations of motion for the bosonic degrees of freedom are given by We specify initial conditions to solve the system of time evolution equations (87) and (89). The method outlined above allows us to use Gaussian initial states for the fermions. In the following, we focus on the vacuum of the Hamiltonian Since H f,0 is quadratic, the ground state and the dispersion relation can be determined analytically, cf. Appendix C. Given an optical lattice with N elementary cells, i.e. 2N lattice sites, the dispersion relation is given by two bands ±ω q with with q ∈ {0, . . . , N − 1}. The corresponding mode function expansion of the fermionic field operator reads and the momentum space creation/annihilation operators are defined with respect to the fully filled lower band |vac according to a q |vac = c q |vac = 0. The bosonic samples are prepared in an excited eigenstate of determined by the number of bosonic atoms on each site 2 = b † n b n + d † n d n and the eigenvalue of the operator The initial conditions for the bosons need to be chosen such that Gauss's law is fulfilled. Since the bosonic degrees of freedom are highly occupied, we approximate the initial Wigner function by In fact, the corrections to the exact Wigner function are O(1/ ) [57]. The explicit values of B 0,n and D 0,n are specified in the next section. To initiate the dynamic evolution, the system is then quenched to an interacting field theory governed by the Hamiltonian (46).

VI. SCHWINGER PAIR PRODUCTION AND STRING BREAKING
In the following we discuss two fundamental phenomena of high-energy physics described by Schwinger model and whose dynamics might be addressed in the cold-atom framework presented. First we discuss how the cold atom system approaches the continuum results of Schwinger pair production and then we concentrate on parameter sets that characterize given experimental systems.

A. Theoretical results
Schwinger Pair Production. In quantum electrodynamics the presence of a sufficiently strong electric field results in the spontaneous breakdown of the vacuum by the emission of charged particle-antiparticle pairs (Schwinger effect) [63][64][65]. This fundamental process has not been experimentally observed yet due to the large required field strength of the order of the critical electric field E c = M 2 /g. The observation of this effect in the cold-atom framework, however, seems to be feasible with current technology as discussed in section IV.
To study the Schwinger effect in the cold-atom framework, we consider a one-dimensional lattice with 2N lattice sites, periodic boundary conditions and finite spin magnitude = 1 2 (N b + N d ). An initially constant electric field E/E c = 1 in QED then corresponds to an initial configuration with a bosonic species imbalance |I(0)| = |N b − N d | = 2M 2 /g 2 . We first solve the equations of motion in the limit → ∞ for g/M = 0.1, a S M = 0.005 and N = 512. We checked that our results are insensitive to changes of both the infrared and the ultraviolet cutoff.
In fact, the information of the fermionic sector is encoded in the correlation function F nm . Even though the concept of a particle number is not uniquely defined in an interacting theory [66], it is useful to define a quasiparticle distribution n(k) from F nm [32,34,67], cf. also Appendix C. We display the time-evolution of the total particle number N = k n(k) in Fig. 5. The production rate at early times, when the backreaction of produced particles does not yet substantially influence the dynamics, coincides with the analytically known result [32,67] At later times, the backreaction of particles becomes important and leads to the expected deviations from the analytic curve. We find phases in which particle production terminates and plateaus are formed [32,34].
In the cold-atom setup no fundamental particles are produced since the number of atoms is fixed. However, the physics of pair production is still encoded in the correlation function F nm since the staggered structure of the cold-atom setup results in the representation of fermions on even/odd sites as particles/antiparticles. Accordingly, the hopping of fermionic atoms between neighboring sites which generates correlations F nm can be interpreted as pair production. As the cold-atom setup shows a truncation error of O(δρ/ ) compared to QED it is interesting to investigate this error as a function of . In Fig. 5, we demonstrate the convergence of the cold-atom behavior towards the QED result upon increasing . For the chosen parameters and = 2500, we still observe sizable quantitative deviations from the QED result whereas this discrepancy further decreases for increasing values of .
Besides studying fermionic observables based on F nm , it is also instructive to evaluate the bosonic species imbalance I(t) = N b (t)−N d (t) which is related to the QED electric field according to E(t) = gI(t)/2. In Fig. 5, we display the time-evolution of the electric field. Starting from QED ( → ∞), we observe the expected behavior according to which the production and subsequent acceleration of particle-antiparticle pairs results in plasma oscillations [32,34,68]. Accordingly, the electric field decreases as the particle number increases and particle creation effectively terminates once the field drops below a certain level, corresponding to the plateaus in the particle number in Fig. 5.
In the cold-atom setup, the physics of plasma oscillations is observed as well. The fermionic hopping reduces the initial bosonic species imbalance I(0) = 2M 2 /g 2 > 0 until it changes sign and reaches a local minimum I(t min ) < 0. Subsequently, the species imbalance increases again, changes sign and reaches a local maximum and so forth. As for the particle number, we observe that the cold-atom behavior converges towards the QED results upon increasing the value of .
The results in Fig. 5 are all based on system sizes of N = 512 for which infrared artifacts are suppressed. In the following we study the influence of the system size and display the time-evolution of the bosonic species imbalance I(t) for different values of N and fixed = 10 5 in Fig. 6. Whereas the behavior remains the same qualitatively, the actual quantitative behavior might substantially change upon decreasing the value of N . For N being too small, we observe oscillations on top of the plasma oscillations which can be attributed to the finite momentum resolution. Moreover, even though a reasonable agreement between QED and the coldatom setup is found for the first oscillation period for N = 512, we still observe sizable deviation at later times.

String Breaking.
The physics of confinement in the theory of quantum chromodynamics (QCD) manifests itself by the formation of a string between two external, static quarks. This confining string can break in theories with dynamical fermions by the production of charged particle-antiparticle pairs which result in a screening of the static sources [69][70][71][72][73][74]. QED in one spatial dimension shares important aspect of dynamical string breaking and, therefore, serves as a toy model for addressing related questions [33,38,40,41,44]. To study dynamical string breaking in QED in one spatial dimension we prepare two static charges ±Q located at ±d/2 on the spatial lattice with 2N lattice sites. The corresponding electric field between the charges is given by E 0 = Q while it vanishes outside. In the cold-atom setup, this corresponds to a bosonic species imbalance of I(0) = 2Q/g inside the string |x| < d/2 whereas it vanishes outside of it.
We first make contact to the corresponding QED literature [33,44] by considering the limit → ∞ and choosing g/M = 1, a S M = 0.1 and N = 1024. To this end, we study the time-evolution of the electric field E n for d/a S = 287 and display different instances of time in Fig. 7. Starting from the initial field configuration, the field energy is transferred to the fermionic sector by particle-antiparticle production such that the amplitude decreases. The dynamics is such that the opposite charges are produced locally on top of each other and are then accelerated by the electric field. Depending on the value of d, the initial string may or may not con- tain enough energy to produce the required charges ±Q to screen the external charges. In the first scenario the string does not break completely and in the former case the string starts oscillating.
In Fig. 8 we display the electric field in the center of the string, and we choose d such that the produced amount of charge exactly screens the external charges, which is attributed to the phenomenon of string breaking. Considering the cold-atom setup, the finite value of then again introduces deviations from the QED behavior. Most notably, we observe that the breaking of the string happens already for smaller distances d CA < d QED for the same parameters g/M = 1, a S M = 0.1 and N = 1024. As expected, we observe convergence towards the QED results upon increasing the value of .
Finally, we consider fermionic observables which are defined in terms of the correlation function F nm . Unlike in the Schwinger mechanism, however, we may observe charge separation directly owing to the spatially inhomogeneous configuration. Accordingly, we focus on the expectation value of the charge density. More specifically, we consider the average charge density on two neighboring lattice sites in order to coarse grain the staggered structure, which is an artifact of the chosen fermion discretization In Fig. 9 we display the time evolution of the charge densityq n . As described previously, the dynamical charges are produced on top of each other such that the total charge density vanishes initially. The dynamical charges are then separated by the existing field such that positive charges are accelerated towards −Q and negative charge towards +Q. As the dynamical charges cannot be considered as hardcore particles, the charge density spreads beyond the static charges resulting in the outwards directed parts of the charge density. Accordingly, the external charges are gradually screened and finally result in the breaking of the string. At asymptotic times, the external charges are then supposed to become screened by an exponential cloud of dynamical fermions [24,33,44,75].

B. Experimental protocol
After discussing strong field QED in the continuum limit, we now present an experimental protocol of initialization, evolution and detection that allows us to observe QED with cold atoms. It is important to note that it does not suffice to only provide the required symmetry of the Hamiltonian in order to quantum simulate a gauge theory but, equally important, also the initial conditions have to fulfill Gauss's law as accurately as possible. The second requirement can only be guaranteed to a certain degree, however, the presented preparation scheme is consistent with the initial conditions chosen for the theoretical description. According to section III D, the fermions can be prepared in the lowest band via a subsequential adiabatic ramp of V f 2,α and then V f 1,α . If the wavelength of this lattice is well chosen, the bosonic links with N B atoms are also directly prepared at the intersection between the fermionic sites. In particular, the chosen wavelength is supposed to be red detuned for lithium and blue detuned for sodium. To apply an initial electric field according to (94), the bosons need to be prepared in a staggered structure of alternating imbalance. This can be achieved by controlling the imbalance with a linear coupling between the two bosonic states, e.g., via radio-frequency coupling or two photon microwave coupling. The detuning of the linear coupling for every second site is controlled by utilizing a species selective standing light wave with twice the lattice period for the bosons. Properly chosen, one can implement a π pulse for every second site that leads to the required 'staggered' imbalance of the bosons. The subsequent dynamics of the system is initiated with a quench of the mass term from being far off-resonant.
We start our benchmarking with the parameters presented in section IV, which correspond to g/M = 2.6, a S M = 0.05 and N = 100. For these initial conditions we can benchmark a possible experimental realization via the functional integral approach. This allows us to investigate the role of the experimentally relevant parameters on the physics of the Schwinger effect, especially the spin magnitude = N B /2 and the coupling strength g. To check the convergence towards the lattice QED result for given parameters we also vary the spin magnitude, = 10, 20, ∞. Since the experimental parameter allow for the exploration of the strong coupling regime g/M > 1, the range of validity of the theoretical treatment has to be further investigated [76]. We expect, however, that the experiment shows the same qualitative behavior as shown in Fig. 10.
In Fig. 10a we display the time evolution of particle FIG. 10: Particle production in 23 Na-6 Li setup: The initial production of particles is driven by χBF NB with an initial imbalance of I(0) = 0.3, which subsequently levels off so that the particle number reaches a plateau. We observe quick convergence towards QED as function of , where the parameters for this simulation are g/M = 2.6, aSM = 0.05 and N = 100.
number per lattice site, where we use again the adiabatic definition from appendix C. First, we observe that particle production happens initially on time scales of the order of χ BF N B , which is short compared to the limiting particles losses, and thus accessible in the experiment. As to be expected from χ BF N B > ∆ we also observe initial oscillations, which are smeared on times t > h/∆. The particle number per lattice site then reaches a plateau with N (t)/N 10 −3 . For the given parameters, we find convergence towards the QED results already for = 20 which, as compared to the ideal-typical parameters from the previous section, can be traced back to the larger value of the coupling g/M = 2.6. In general it turns out that the required value of is inversely proportional to g/M to obtain convergence towards the QED result, cf. also the simulations of dynamical string breaking in Fig. 8. Owing to the fact that realistic values of the spin magnitude are of the order of = O(100), we can expect genuine QED behavior for the proposed experimental setup. The fluctuations in the number of bosons, which is expected to be of the order of a few percent, are not expected to substantially alter the reported behavior on the time scales considered. We also display the electric field E(t) = gI(t)/2 as determined by the species imbalance Fig. 10b. The production of particles results in a decrease of the species imbalance, which quickly drops below the critical value so that particle production terminates. Conceptually the momentum distribution of the produced particles can be read out precisely via band mapping [30], however, this is very challenging for the given parameters as one can deduce from Fig. 10. Nevertheless, the underlying physics of particle production can be already accessed from the integrated number of produced particles. Nevertheless, the underlying physics of particle production can be already accessed from the integrated number of produced particles. According to Fig. 10, around 0.2 particles have to be detected on average for a lattice with N = 100 sites. Current experimental se-tups allow realizing ≈ 30 copies of the cold atom QED system by employing an array of one dimensional traps. Thus one expects integrated O(10) particles which can be detected with fluorescence in a subsequent magnetooptical trap for which single particle resolution is well established [28]. While the detection of the produced particles is very challenging, the bosonic species imbalance, i.e. the electric field, changes significantly. Thus, the Schwinger effect in a cold atomic setup can also be observed by measuring the integrated boson imbalance via standard absorption imaging techniques.
The coupling strength g/M is another experimentally relevant parameter of importance. Accordingly, we study the dependence of the particle number on choosing slightly different couplings g/M = 0.4, 1.0, 2.6 for fixed parameters a S M = 0.05, N = 100 and → ∞ in Fig. 11. We note that we have to adapt the initial species imbal- ance I(0) = 2M 2 /g 2 in order to provide an initial critical field E cr in each case. Most notably, the simulations indicate that the first plateau in the particle number is a stable signature for the non-trivial interplay of the electric field and the produced fermions. We find, however, that the number of produced particles at the first plateau is inversely proportional to the ratio g/M . As expected, smaller couplings g/M result in a slowdown of dynamics such that the first plateau is reached at later times.

VII. CONCLUSION
This work is meant as a guide towards a first largescale quantum simulation of a lattice gauge theory with dynamical gauge fields. For that purpose, we concentrate on one of the simplest, yet highly nontrivial, example of QED in 1+1 space-time dimensions and provide strong evidence that present-day experimental resources and protocols are sufficient to observe the dynamical phenomena of Schwinger pair production and string breaking in the laboratory using ultracold atoms.
Our results point out that experimental realizations using coherent many-body states residing on the links of an optical lattice can be highly efficient for quantum simulations of such high-energy particle physics phenomena. This represents a paradigmatic change in view of the large number of studies in the literature focusing on a small number of atoms per link. To substantiate our findings, we exploit the long-term experience that has been gained with the engineering and manipulation of related setups of fermions interacting with coherent samples of bosonic atoms.
For the example of a Bose-Fermi mixture of 23 Na and 6 Li atoms, characterized by a plethora of potentially gauge-symmetry breaking interactions, we apply external potentials and fields such that one ends up with a lattice Hamiltonian with local gauge invariance. In this way, the microscopic parameters describing the bosonic and fermionic atom degrees of freedom are connected with the parameters describing the gauge field theory.
We use the experimentally available parameter range of the cold-atom system in benchmark calculations and convergence against the full QED results is observed. The very detailed comparisons of the real-time dynamics of the atomic system in the large boson number regime are possible using powerful functional integral techniques, which complement exact diagonalization or tensor network methods that are applicable in the small boson number regime.
A future experimental implementation of gauge symmetries in a flexible cold-atom setup can explore new parameter ranges and phenomena, even beyond what is realized by nature so far. While the gauge coupling of QED is weak, with α em 1/137 in nature's three spatial dimensions, studies at stronger couplings in various dimensions would be extremely interesting. No conventional computational technique has so far been able to predict the real-time dynamics of QED or QCD for couplings of order one -despite the fact that this is a crucial missing link in our understanding of the thermalization process of QCD as explored in relativistic heavy-ion collision experiments. The experimental setup discussed in this work may provide already access to answering important aspects of longstanding open questions.
In the vicinity of the potential minimum x b,n , we may expand the potential in a Taylor series to quadratic order The minima in the fermionic potential are determined by For the fermionic atoms, there are two different types of local minima and we need to expand the optical potential in a Taylor series to quadratic order around both of them.

(B7)
Appendix C: Particle number distribution In this appendix, we define and determine the particle number distribution n(q) in terms of the correlation function F nm = vac| [ψ n , ψ † m ] |vac . To this end, we consider the fermionic part of the Kogut-Susskind Hamiltonian (1) We may diagonalize the Hamiltonian by treating the link variables U n as c-number background for the fermions, as it is also done in the functional integral approach. To this end, we define the Fourier transformation according to We note that the system is still translation invariant over two lattice sites if we study Schwinger pair production, i.e. U n = U n+2l with l ∈ {0, . . . , N − 1}. Denoting the even links by U even = U 2n and the odd links by U odd = U 2n+1 , we have with U A = (U even + U odd )/2 and U B = (U even − U odd )/2. The Hamiltonian can then be written in matrix notation according to Fermions without a gauge field correspond to U n = 1 leading to U A = 1 and U B = 0. The eigenvalues of this Hamiltonian are given by ω ±,q = ±ω q with The corresponding normalized eigenvectors are In fact, every q-mode is diagonalized by a unitary transformation matrix U q = (u + q , u − q ). This defines quasiparticle creation/annihilation operators with respect to the instantaneous vacuum state |Ω , fulfilling a q |Ω = c q |Ω = 0. The Hamiltonian is then given by We define the quasi-particle distribution function n(q) as the expectation value of the instantaneous number operator n(q) ≡ vac| a † q a q + c † q c q |vac , where the asymptotic ground state |vac is given by the Hamiltonian with U n = 1. The expectation value of the Hamiltonian is The two contributions of n(q) are found by employing the Bogoliubov transformation (C8) such that vac| a † q a q |vac = |M q | 2 2ω q (ω q − π q ) vac|ψ † q ψ q |vac + M q 2ω q vac|ψ † q ψ q+N |vac + M * q 2ω q vac|ψ † q+N ψ q |vac + ω q − π q 2ω q vac|ψ † q+N ψ q+N |vac , (C12) and vac| c † q c q |vac = |M q | 2 2ω q (ω q + π q ) vac|ψ q ψ † q |vac The Fourier transformation of the correlation function F mn is determined bỹ with q, q ∈ {0, . . . , 2N − 1}. Accordingly, n(q) can be written as where the energy density in Fourier space is given by The total particle number is then found by summing over all Fourier modes N = q n(q) .
In the cold atom system, we proceed analogously but replace the link operators by the Schwinger bosons U n → [ ( + 1)] −1/2 b * n d n . As the bosonic degrees of freedom are again considered as c-numbers, we obtain very similar expressions for the particle number.