Implementing quantum electrodynamics with ultracold atomic systems

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 figure 2. In fact, the potential energy contributions (11) can be split into an axial and radial part,

$$\begin{eqnarray}&&{V}_{\alpha }^{s}({\bf{x}})={V}_{\parallel ,\alpha }^{s}(x)+{V}_{\perp ,\alpha }^{s}{{\rm{e}}}^{-({y}^{2}+{z}^{2})/{l}_{s,\alpha }^{2}}.\end{eqnarray} \tag{ 13 }$$

Here, ${l}_{s,\alpha }$ denotes the radial confinement length scale, where we distinguish bosons and fermions by the superscript $s\in \{b,f\}$. Owing to tight radial confinement, the three-dimensional system is effectively rendered one-dimensional and we employ the product form

$$\begin{eqnarray}&&{\phi }_{\alpha }({\bf{x}})={\varphi }_{b}(y){\varphi }_{b}(z){\phi }_{\alpha }(x),\end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray}&&{\psi }_{\alpha }({\bf{x}})={\varphi }_{f}(y){\varphi }_{f}(z){\psi }_{\alpha }(x),\end{eqnarray} \tag{ 14b }$$

where ${\varphi }_{s}(y)$ and ${\varphi }_{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.

Zoom In Zoom Out Reset image size

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 the right graph of figure 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

$$\begin{eqnarray}&&{V}_{\parallel ,\alpha }^{b}(x)={V}_{1,\alpha }^{b}{\cos }^{2}\left(\displaystyle \frac{2\pi x}{a}\right),\end{eqnarray} \tag{ 15a }$$

$$\begin{eqnarray}&&{V}_{\parallel ,\alpha }^{f}(x)={V}_{1,\alpha }^{f}{\sin }^{2}\left(\displaystyle \frac{2\pi x}{a}\right)+{V}_{2,\alpha }^{f}{\cos }^{2}\left(\displaystyle \frac{\pi x}{a}\right).\end{eqnarray} \tag{ 15b }$$

The amplitudes ${V}_{i,\alpha }^{s}$, 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 [16], such a construction is spin-independent and, hence, distinct from e.g. [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 fermionic degrees of freedom and focus on the two lowest energy bands. By tuning of the laser amplitude, we may choose ${V}_{i,\alpha }^{f}$ such that two Wannier functions ${w}_{\alpha ,n^{\prime} p}^{f}(x)$ are obtained which are sufficiently localized in the left (p = L) and right (p = R) minimum of the elementary cell $n^{\prime} \in \{0,\ldots ,N-1\}$ with positive integer N of the optical lattice, respectively. The corresponding expansion of the fermionic field operator reads

$$\begin{eqnarray}&&{\psi }_{\alpha }(x)=\displaystyle \sum _{n^{\prime} ,p}{w}_{\alpha ,n^{\prime} p}^{f}(x){\psi }_{\alpha ,n^{\prime} p}.\end{eqnarray} \tag{ 16 }$$

We note that the total number of bosonic/fermionic lattice sites is $2N$ and we will label them by $n\in \{0,\ldots ,2N-1\}$. Similarly, we may expand the bosonic field operators, ${\phi }_{\alpha }(x)={\sum }_{n^{\prime} ,p}{w}_{\alpha ,n^{\prime} p}^{b}(x){\phi }_{\alpha ,n^{\prime} p}$, where the Wannier functions ${w}_{\alpha ,n^{\prime} p}^{b}(x)$ are again localized in the two minima of the elementary cell. In fact, the structure of the superlattice suggests the definition

$$\begin{eqnarray}&&{\psi }_{2n^{\prime} ,\alpha }\equiv {\psi }_{\alpha ,n^{\prime} L},\quad {\psi }_{2n^{\prime} +1,\alpha }\equiv {\psi }_{\alpha ,n^{\prime} R},\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&{\phi }_{2n^{\prime} ,\alpha }\equiv {\phi }_{\alpha ,n^{\prime} L},\quad {\phi }_{2n^{\prime} +1,\alpha }\equiv {\phi }_{\alpha ,n^{\prime} R}.\end{eqnarray} \tag{ 17b }$$

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 (see also section 4).

In the previous section, we reviewed how the potential energy (11) can be used to generate the staggered lattice structure. Moreover, it 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}_{\alpha \beta \gamma \delta }^{s}$ for $s\in \{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

$$\begin{eqnarray}&&{V}^{s}({{\bf{x}}}_{1},{{\bf{x}}}_{2})=\delta ({{\bf{x}}}_{1}-{{\bf{x}}}_{2})\displaystyle \sum _{{{ \mathcal F }}_{s}}\,{g}_{s,{{ \mathcal F }}_{s}}{{\bf{P}}}_{{{ \mathcal F }}_{s}},\end{eqnarray} \tag{ 18 }$$

where the total spin can take the values ${{ \mathcal F }}_{b}\in \{0,2\}$, ${{ \mathcal F }}_{f}\in \{0,1\}$ and ${{ \mathcal F }}_{{bf}}\in \{1/2,3/2\}$ [51]. The interaction strengths ${g}_{s,{{ \mathcal F }}_{s}}$ are related to the s-wave scattering lengths ${a}_{s,{{ \mathcal F }}_{s}}$ via

$$\begin{eqnarray}&&{g}_{s,{{ \mathcal F }}_{s}}=\displaystyle \frac{2\pi {{\hslash }}^{2}{a}_{s,{{ \mathcal F }}_{s}}}{{M}_{r,s}}.\end{eqnarray} \tag{ 19 }$$

Here ${M}_{r,s}$ denotes the reduced mass of the two scattering partners. In general, the projector ${{\bf{P}}}_{{ \mathcal F }}$ for two particles with individual spins f₁ and f₂ on the subspace with total spin ${ \mathcal F }$ can be written as

$$\begin{eqnarray}&&{{\bf{P}}}_{{ \mathcal F }}=\displaystyle \sum _{M}| {f}_{1},{f}_{2};{ \mathcal F },M\rangle \langle {f}_{1},{f}_{2};{ \mathcal F },M| ,\end{eqnarray} \tag{ 20 }$$

where $M\in \{-{ \mathcal F },-{ \mathcal F }+1,\ldots ,{ \mathcal F }-1,{ \mathcal F }\}$ are the possible magnetic quantum numbers. We may relate the interaction strengths ${g}_{s,{{ \mathcal F }}_{s}}$ to the constants appearing in the interaction part of the Hamiltonian (12) according to

$$\begin{eqnarray}&&{g}_{\alpha \beta \gamma \delta }^{s}=\displaystyle \sum _{{{ \mathcal F }}_{s}}\displaystyle \sum _{M}\,{g}_{s,{{ \mathcal F }}_{s}}\langle {f}_{1};\alpha ;{f}_{2};\beta | {f}_{1},{f}_{2};{{ \mathcal F }}_{s},M\rangle \langle {f}_{1},{f}_{2};{{ \mathcal F }}_{s},M| {f}_{1};\gamma ;{f}_{2};\delta \rangle .\end{eqnarray} \tag{ 21 }$$

Here, $\langle {f}_{1};\alpha ;{f}_{2};\beta | {f}_{1},{f}_{2};{{ \mathcal F }}_{s},M\rangle$ are the Clebsch–Gordan coefficients for coupling the individual spins f₁ and f₂ to the total spin ${{ \mathcal F }}_{s}$. Specifically, we have ${f}_{1}={f}_{2}=1$ for boson–boson interactions (s = b), ${f}_{1}={f}_{2}=1/2$ for fermion–fermion interaction (s = f) and ${f}_{1}=1/2$, ${f}_{2}=1$ for boson–fermion interaction (s = bf).

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 ${\bf{n}}=({n}_{1},{n}_{2},{n}_{3},{n}_{4})$ denotes the site indices and ${\boldsymbol{\mu }}=(\alpha ,\beta ,\gamma ,\delta )$ the magnetic quantum numbers, we can write for the interaction part of the Hamiltonian:

$$\begin{eqnarray}&&{H}_{I}=\displaystyle \frac{1}{2}\displaystyle \sum _{{\bf{n}},{\boldsymbol{\mu }}}{U}_{{\bf{n}}}^{b}{g}_{{\boldsymbol{\mu }}}^{b}{\phi }_{{n}_{1},\alpha }^{\dagger }{\phi }_{{n}_{2},\beta }^{\dagger }{\phi }_{{n}_{4},\delta }{\phi }_{{n}_{3},\gamma }+\displaystyle \frac{1}{2}\displaystyle \sum _{{\bf{n}},{\boldsymbol{\mu }}}{U}_{{\bf{n}}}^{f}{g}_{{\boldsymbol{\mu }}}^{f}{\psi }_{{n}_{1},\alpha }^{\dagger }{\psi }_{{n}_{2},\beta }^{\dagger }{\psi }_{{n}_{4},\delta }{\psi }_{{n}_{3},\gamma }+\displaystyle \frac{1}{2}\displaystyle \sum _{{\bf{n}},{\boldsymbol{\mu }}}{U}_{{\bf{n}}}^{{bf}}{g}_{{\boldsymbol{\mu }}}^{{bf}}{\psi }_{{n}_{1},\alpha }^{\dagger }{\phi }_{{n}_{2},\beta }^{\dagger }{\phi }_{{n}_{4},\delta }{\psi }_{{n}_{3},\gamma }.\end{eqnarray} \tag{ 22 }$$

The coupling constants ${U}_{{\bf{n}}}^{s}$ 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 figure 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.

Zoom In Zoom Out Reset image size

All bosonic states are prepared in ${\alpha }_{b}=\{-1,0\}$ states whereas the fermionic degrees of freedom are generated in the staggered configuration with ${\alpha }_{f}=1/2$ on even sites and ${\alpha }_{f}=-1/2$ on odd sites. As a consequence, interactions including the ${\alpha }_{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.

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}_{{\bf{n}}}^{b}\ne 0$ for ${\bf{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

$$\begin{eqnarray}&&{H}_{I}^{b}=\displaystyle \frac{1}{2}\displaystyle \sum _{n,{\boldsymbol{\mu }}}{U}_{{\bf{n}}}^{b}{g}_{{\boldsymbol{\mu }}}^{b}{\phi }_{n,\alpha }^{\dagger }{\phi }_{n,\beta }^{\dagger }{\phi }_{n,\delta }{\phi }_{n,\gamma },\end{eqnarray} \tag{ 23 }$$

where all entries of ${\boldsymbol{\mu }}=(\alpha ,\beta ,\gamma ,\delta )$ may take values ${\alpha }_{b}\in \{-1,0\}$. Again, we note that we disregard terms including the magnetic substates ${\alpha }_{b}=1$ which are excluded by the spin conservation if initialized accordingly. The interaction term (23) reads

$$\begin{eqnarray}&&{H}_{I}^{b}=\displaystyle \frac{1}{2}\displaystyle \sum _{n}{U}_{{\bf{n}}}^{b}\left({g}_{b,2}\,{\phi }_{n,-1}^{\dagger }{\phi }_{n,-1}^{\dagger }{\phi }_{n,-1}{\phi }_{n,-1}+\displaystyle \frac{{g}_{b,0}+2{g}_{b,2}}{3}{\phi }_{n,0}^{\dagger }{\phi }_{n,0}^{\dagger }{\phi }_{n,0}{\phi }_{n,0}+2{g}_{b,2}\,{\phi }_{n,0}^{\dagger }{\phi }_{n,-1}^{\dagger }{\phi }_{n,-1}{\phi }_{n,0}\right),\end{eqnarray} \tag{ 24 }$$

where the coupling constants result from (21) and we assumed that the overlap integrals ${U}_{{\bf{n}}}^{b}$ are the same for all terms. For later convenience, we denote the bosonic degrees of freedom on even sites as ${\phi }_{2n,-1}\equiv {d}_{2n}$ and ${\phi }_{2n,0}\equiv {b}_{2n}$, whereas we interchange their role on odd sites such that ${\phi }_{2n+1,-1}\equiv {b}_{2n+1}$ and ${\phi }_{2n+\mathrm{1,0}}\equiv {d}_{2n+1}$. In fact, the bosons b_n and d_n can be understood as Schwinger bosons [52] with the identification

$$\begin{eqnarray}&&{L}_{+,n}={b}_{n}^{\dagger }{d}_{n},\quad {L}_{-,n}={d}_{n}^{\dagger }{b}_{n},\end{eqnarray} \tag{ 25a }$$

$$\begin{eqnarray}&&{L}_{z,n}=\displaystyle \frac{1}{2}({b}_{n}^{\dagger }{b}_{n}-{d}_{n}^{\dagger }{d}_{n}),\end{eqnarray} \tag{ 25b }$$

which constitute a representation of the angular momentum algebra $[{L}_{i,n},{L}_{j,m}]={\rm{i}}{\delta }_{{nm}}{\epsilon }_{{ijk}}{L}_{k,n}$ with ${L}_{\pm ,n}={L}_{x,n}\pm {{\rm{i}}{L}}_{y,n}$. The constraint $2{\ell }={b}_{n}^{\dagger }{b}_{n}+{d}_{n}^{\dagger }{d}_{n}$ is fulfilled because the hopping of bosons between neighboring sites $n\to n\pm 1$ is suppressed. The bosonic intra-species interaction Hamiltonian is then given by

$$\begin{eqnarray}&&{H}_{I}^{b}=\displaystyle \frac{{g}_{b,0}-{g}_{b,2}}{6}{U}_{{\bf{n}}}^{b}\displaystyle \sum _{n}\,{L}_{z,n}^{2}+{{\rm{\Delta }}}_{b,0}\displaystyle \sum _{n}{(-1)}^{n}{L}_{z,n},\end{eqnarray} \tag{ 26 }$$

where we disregarded an irrelevant constant and introduced the abbreviation ${{\rm{\Delta }}}_{b,0}\equiv (2{\ell }-1)({g}_{b,0}-{g}_{b,2}){U}_{{\bf{n}}}^{b}/6$.

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}_{{\bf{n}}}^{f}\ne 0$ only for ${\bf{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

$$\begin{eqnarray}&&{H}_{I}^{f}=\displaystyle \sum _{n}{U}_{{\bf{n}}}^{f}{g}_{f,0}{\psi }_{n,1/2}^{\dagger }{\psi }_{n,-1/2}^{\dagger }{\psi }_{n,-1/2}{\psi }_{n,1/2}.\end{eqnarray} \tag{ 27 }$$

In general, this four-fermion interaction term influences the dynamics. However, the contribution can be written as a density–density interaction

$$\begin{eqnarray}&&{H}_{I}^{f}=\displaystyle \sum _{n}{U}_{{\bf{n}}}^{f}{g}_{f,0}{\rho }_{n,1/2}{\rho }_{n,-1/2}\,\end{eqnarray} \tag{ 28 }$$

between ${\alpha }_{f}=-1/2$ and ${\alpha }_{f}=1/2$ particles with density operators ${\rho }_{n,\pm 1/2}={\psi }_{n,\pm 1/2}^{\dagger }{\psi }_{n,\pm 1/2}$. Restricting ourselves to an initial state $| {\rm{\Psi }}\rangle$ with only ${\alpha }_{f}=1/2$ particles on even sites and only ${\alpha }_{f}=-1/2$ particles on odd sites, one immediately finds that ${H}_{I}^{f}| {\rm{\Psi }}\rangle =0$. Consequently, this four-fermion interaction does not contribute to the time evolution due to an appropriate initial-state preparation.

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

$$\begin{eqnarray}&&{H}_{{I}_{{se}}}^{{bf}}=\displaystyle \frac{1}{2}\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}}^{{bf}}{g}_{{\boldsymbol{\mu }}}^{{bf}}({\psi }_{2n,\alpha }^{\dagger }{\phi }_{2n,\beta }^{\dagger }{\phi }_{2n,\delta }{\psi }_{2n+1,\gamma }+{\psi }_{2n,\alpha }^{\dagger }{\phi }_{2n-1,\beta }^{\dagger }{\phi }_{2n-1,\delta }{\psi }_{2n-1,\gamma }+{\rm{h}}{\rm{.}}{\rm{c}}.),\end{eqnarray} \tag{ 29 }$$

with ${\boldsymbol{\mu }}=(\alpha ,\beta ,\gamma ,\delta )=(-1/2,0,1/2,-1)$. The first term corresponds to figure 4(a) whereas the second term is shown in figure 4(b). According to (21), the coupling constant for this specific scattering process is given by

$$\begin{eqnarray}&&{g}_{{\boldsymbol{\mu }}}^{{bf}}=\displaystyle \frac{\sqrt{2}}{3}({g}_{{bf},3/2}-{g}_{{bf},1/2}).\end{eqnarray} \tag{ 30 }$$

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

$$\begin{eqnarray}&&{\psi }_{2n}\equiv {\psi }_{2n,-1/2},\quad {\psi }_{2n+1}\equiv {\psi }_{2n+\mathrm{1,1}/2}.\end{eqnarray} \tag{ 31 }$$

Employing the Schwinger boson representation and taking into account that the overlap integral ${U}_{{\bf{n}}}^{{bf}}$ does not depend on the specific lattice site n, the spin exchange Hamiltonian can be written as

$$\begin{eqnarray}&&{H}_{{I}_{{se}}}^{{bf}}=\displaystyle \frac{1}{2}{U}_{{\bf{n}}}^{{bf}}{g}_{{\boldsymbol{\mu }}}^{{bf}}\displaystyle \sum _{n}({\psi }_{n}^{\dagger }{L}_{+,n}{\psi }_{n+1}+{\rm{h}}.{\rm{c}}.).\end{eqnarray} \tag{ 32 }$$

The elastic scattering processes, on the other hand, are given by

$$\begin{eqnarray}{H}_{{I}_{{el}}}^{{bf}} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{n^{\prime} ,\beta }{U}_{{\bf{n}}}^{{bf}}{g}_{{\boldsymbol{\mu }}}^{{bf}}{\psi }_{2n}^{\dagger }{\phi }_{2n,\beta }^{\dagger }{\phi }_{2n,\beta }{\psi }_{2n}\\ & & +\displaystyle \frac{1}{2}\displaystyle \sum _{n^{\prime} ,\beta }{U}_{{\bf{n}}}^{{bf}}{g}_{{\boldsymbol{\mu }}}^{{bf}}{\psi }_{2n}^{\dagger }{\phi }_{2n-1,\beta }^{\dagger }{\phi }_{2n-1,\beta }{\psi }_{2n}\\ & & +\displaystyle \frac{1}{2}\displaystyle \sum _{n^{\prime} ,\beta }{U}_{{\bf{n}}}^{{bf}}{g}_{{\boldsymbol{\mu }}}^{{bf}}{\psi }_{2n+1}^{\dagger }{\phi }_{2n,\beta }^{\dagger }{\phi }_{2n,\beta }{\psi }_{2n+1}\\ & & +\displaystyle \frac{1}{2}\displaystyle \sum _{n^{\prime} ,\beta }{U}_{{\bf{n}}}^{{bf}}{g}_{{\boldsymbol{\mu }}}^{{bf}}{\psi }_{2n+1}^{\dagger }{\phi }_{2n+1,\beta }^{\dagger }{\phi }_{2n+1,\beta }{\psi }_{2n+1},\end{eqnarray} \tag{ 33 }$$

where $\beta \in \{-1,0\}$. We note that the coupling constants ${g}_{{\boldsymbol{\mu }}}^{{bf}}$ still depend on the magnetic substates and are therefore not identical for the different terms. Moreover, we find that each ${U}_{{\bf{n}}}^{{bf}}$ is independent of n in (33), and identical in the first and second line (further denoted by ${U}_{{\bf{n}}1}^{{bf}}$) as well as in the third and fourth line (further denoted by ${U}_{{\bf{n}}3}^{{bf}}$), see appendix A. For the first term in (33) with ${\boldsymbol{\mu }}=(-1/2,\beta ,-1/2,\beta )$, we obtain

$$\begin{eqnarray}&&{H}_{{I}_{{el},1}}^{{bf}}=\displaystyle \frac{{g}_{{bf},1/2}+2{g}_{{bf},3/2}}{6}\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}1}^{{bf}}{\psi }_{2n}^{\dagger }{b}_{2n}^{\dagger }{b}_{2n}{\psi }_{2n}+\displaystyle \frac{{g}_{{bf},3/2}}{2}\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}1}^{{bf}}{\psi }_{2n}^{\dagger }{d}_{2n}^{\dagger }{d}_{2n}{\psi }_{2n},\end{eqnarray} \tag{ 34 }$$

where we used (21) again. The second term in (33) is the same as the first one upon replacing ${b}_{2n}\to {d}_{2n-1}$ and ${d}_{2n}\to {b}_{2n-1}$. The third term in (33), however, is different owing to ${\boldsymbol{\mu }}=(1/2,\beta ,1/2,\beta )$ corresponding to the different fermionic parity and reads

$$\begin{eqnarray}&&{H}_{{I}_{{el},3}}^{{bf}}=\displaystyle \frac{{g}_{{bf},1/2}+2{g}_{{bf},3/2}}{6}\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}3}^{{bf}}{\psi }_{2n+1}^{\dagger }{b}_{2n}^{\dagger }{b}_{2n}{\psi }_{2n+1}+\displaystyle \frac{2{g}_{{bf},1/2}+{g}_{{bf},3/2}}{6}\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}3}^{{bf}}{\psi }_{2n+1}^{\dagger }{d}_{2n}^{\dagger }{d}_{2n}{\psi }_{2n+1}.\end{eqnarray} \tag{ 35 }$$

The fourth term in (33) is the same as the third one upon replacing ${b}_{2n}\to {d}_{2n+1}$ and ${d}_{2n}\to {b}_{2n+1}$. Employing the Schwinger boson representation ${b}_{n}^{\dagger }{b}_{n}={\ell }+{L}_{z,n}$ and ${d}_{n}^{\dagger }{d}_{n}={\ell }-{L}_{z,n}$, the first term (34) can be written as

$$\begin{eqnarray}&&{H}_{{I}_{{el},1}}^{{bf}}=\displaystyle \frac{{g}_{{bf},1/2}-{g}_{{bf},3/2}}{6}\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}1}^{{bf}}{\psi }_{2n}^{\dagger }{\psi }_{2n}{L}_{z,2n}+\displaystyle \frac{{g}_{{bf},1/2}+5{g}_{{bf},3/2}}{6}{\ell }\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}1}^{{bf}}{\psi }_{2n}^{\dagger }{\psi }_{2n}.\end{eqnarray} \tag{ 36 }$$

Similarly, the third term (35) is given by

$$\begin{eqnarray}&&{H}_{{I}_{{el},3}}^{{bf}}=\displaystyle \frac{{g}_{{bf},3/2}-{g}_{{bf},1/2}}{6}\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}3}^{{bf}}{\psi }_{2n+1}^{\dagger }{\psi }_{2n+1}{L}_{z,2n}+\displaystyle \frac{{g}_{{bf},1/2}+{g}_{{bf},3/2}}{2}{\ell }\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}3}^{{bf}}{\psi }_{2n+1}^{\dagger }{\psi }_{2n+1},\end{eqnarray} \tag{ 37 }$$

and similar expressions are also obtained for the second and fourth term by replacing ${L}_{z,2n}\to -{L}_{z,2n\mp 1}$. Accordingly, we obtain a contribution

$$\begin{eqnarray}&&{H}_{{I}_{{el}},{L}_{z}}^{{bf}}={\tilde{g}}_{{bf}}\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}3}^{{bf}}{\psi }_{2n+1}^{\dagger }{\psi }_{2n+1}({L}_{z,2n+1}-{L}_{z,2n})+{\tilde{g}}_{{bf}}\displaystyle \sum _{n^{\prime} }{U}_{{\bf{n}}1}^{{bf}}{\psi }_{2n}^{\dagger }{\psi }_{2n}({L}_{z,2n}-{L}_{z,2n-1})\end{eqnarray} \tag{ 38 }$$

with ${\tilde{g}}_{{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}+{\psi }_{n}^{\dagger }{\psi }_{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:

$$\begin{eqnarray}&&{H}_{{I}_{{el}}}^{{bf}}={U}_{{\bf{n}}3}^{{bf}}[({g}_{{bf},1/2}+{g}_{{bf},3/2}){\ell }+{\tilde{g}}_{{bf}}]\displaystyle \sum _{n^{\prime} }{\psi }_{2n+1}^{\dagger }{\psi }_{2n+1}+{U}_{{\bf{n}}1}^{{bf}}\left(\displaystyle \frac{{g}_{{bf},1/2}+5{g}_{{bf},3/2}}{3}{\ell }+{\tilde{g}}_{{bf}}\right)\displaystyle \sum _{n^{\prime} }{\psi }_{2n}^{\dagger }{\psi }_{2n}.\end{eqnarray} \tag{ 39 }$$

Zoom In Zoom Out Reset image size

Summing up all contributions that originate from the kinetic, potential and interaction terms, the cold atom QED Hamiltonian takes the form

$$\begin{eqnarray}H & = & \displaystyle \frac{{g}_{b,0}-{g}_{b,2}}{6}{U}_{{\bf{n}}}^{b}\displaystyle \sum _{n}\,{L}_{z,n}^{2}+{{\rm{\Delta }}}_{b,0}\displaystyle \sum _{n}{(-1)}^{n}{L}_{z,n}+\displaystyle \frac{1}{2}{U}_{{\bf{n}}}^{{bf}}{g}_{{\boldsymbol{\mu }}}^{{bf}}\displaystyle \sum _{n}({\psi }_{n}^{\dagger }{L}_{+,n}{\psi }_{n+1}+{\rm{h}}.{\rm{c}}.)\\ & & +\displaystyle \sum _{n}({V}_{n}^{f}{\psi }_{n}^{\dagger }{\psi }_{n}+{V}_{n}^{b}{b}_{n}^{\dagger }{b}_{n}+{V}_{n}^{d}{d}_{n}^{\dagger }{d}_{n}).\end{eqnarray} \tag{ 40 }$$

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 energy difference ${{\rm{\Delta }}}_{f}$ according to ${V}_{2n}^{f}={V}_{0}^{f}-{{\rm{\Delta }}}_{f}/2$ and ${V}_{2n+1}^{f}={V}_{0}^{f}+{{\rm{\Delta }}}_{f}/2$, the fermionic potential contribution is given by

$$\begin{eqnarray}&&{H}_{V}^{f}={V}_{0}^{f}\displaystyle \sum _{n}{\psi }_{n}^{\dagger }{\psi }_{n}-\displaystyle \frac{{{\rm{\Delta }}}_{f}}{2}\displaystyle \sum _{n}{(-1)}^{n}{\psi }_{n}^{\dagger }{\psi }_{n}.\end{eqnarray} \tag{ 41 }$$

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 ${{\rm{\Delta }}}_{b,1}$ according to ${V}_{n}^{b}={V}_{0}^{b}+{(-1)}^{n}{{\rm{\Delta }}}_{b,1}/2$ and ${V}_{n}^{d}={V}_{0}^{b}-{(-1)}^{n}{{\rm{\Delta }}}_{b,1}/2$ and using ${L}_{z,n}=({b}_{n}^{\dagger }{b}_{n}-{d}_{n}^{\dagger }{d}_{n})/2$, we obtain

$$\begin{eqnarray}&&{H}_{V}^{b}={{\rm{\Delta }}}_{b,1}\displaystyle \sum _{n}{(-1)}^{n}{L}_{z,n},\end{eqnarray} \tag{ 42 }$$

where we disregarded an irrelevant constant proportional to V^b₀ which only depends on ℓ. Adding the second contribution in (40) and defining ${{\rm{\Delta }}}_{b}\equiv {{\rm{\Delta }}}_{b,0}+{{\rm{\Delta }}}_{b,1}$, we obtain

$$\begin{eqnarray}&&{H}_{V}^{b}+{{\rm{\Delta }}}_{b,0}\displaystyle \sum _{n}{(-1)}^{n}{L}_{z,n}={{\rm{\Delta }}}_{b}\displaystyle \sum _{n}{(-1)}^{n}{L}_{z,n}.\end{eqnarray} \tag{ 43 }$$

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

$$\begin{eqnarray}&&{H}_{0}={{\rm{\Delta }}}_{b}\displaystyle \sum _{n}{(-1)}^{n}\left({L}_{z,n}-\displaystyle \frac{1}{2}{\psi }_{n}^{\dagger }{\psi }_{n}\right),\end{eqnarray} \tag{ 44a }$$

$$\begin{eqnarray}&&{H}_{1}={\chi }_{{BB}}\displaystyle \sum _{n}\,{L}_{z,n}^{2}+{\rm{\Delta }}\displaystyle \sum _{n}{(-1)}^{n}{\psi }_{n}^{\dagger }{\psi }_{n}+\displaystyle \frac{{\chi }_{{BF}}}{2}\displaystyle \sum _{n}({\psi }_{n}^{\dagger }{L}_{+,n}{\psi }_{n+1}+{\rm{h}}.{\rm{c}}.).\end{eqnarray} \tag{ 44b }$$

Here, we introduced the detuning $2{\rm{\Delta }}\equiv {{\rm{\Delta }}}_{b}-{{\rm{\Delta }}}_{f}$, the effective boson self-interaction ${\chi }_{{BB}}$ and the boson–fermion interaction ${\chi }_{{BF}}$,

$$\begin{eqnarray}&&{\chi }_{{BB}}=\displaystyle \frac{{g}_{b,0}-{g}_{b,2}}{6}{U}_{{\bf{n}}}^{b},\end{eqnarray} \tag{ 45a }$$

$$\begin{eqnarray}&&{\chi }_{{BF}}=\displaystyle \frac{\sqrt{2}({g}_{{bf},3/2}-{g}_{{bf},1/2})}{3}{U}_{{\bf{n}}}^{{bf}}.\end{eqnarray} \tag{ 45b }$$

The index n is the same as in (23) and (29) and the overlaps ${U}_{{\bf{n}}}$ are determined in (A11) of the appendix. Upon acting with the unitary transformation $U(t)=\exp (-{{\rm{i}}{H}}_{0}t)$, it is straightforward to show that ${H}_{1}^{\prime }\ ={U}^{\dagger }(t){H}_{1}U(t)$. Performing the canonical transformation ${\psi }_{n}\to {(-{\rm{i}})}^{n}{\psi }_{n}$, we can finally identify ${H}_{1}^{\prime }$ with the quantum link Hamiltonian (8)

$$\begin{eqnarray}&&{H}_{\mathrm{QL}}=\displaystyle \frac{{g}^{2}{a}_{S}}{2}\displaystyle \sum _{n}\,{L}_{z,n}^{2}+M\displaystyle \sum _{n}{(-1)}^{n}{\psi }_{n}^{\dagger }{\psi }_{n}-\displaystyle \frac{{\rm{i}}}{2{a}_{S}\sqrt{{\ell }({\ell }+1)}}\displaystyle \sum _{n}({\psi }_{n}^{\dagger }{L}_{+,n}{\psi }_{n+1}-{\rm{h}}.{\rm{c}}.)\end{eqnarray} \tag{ 46 }$$

where we introduced the abbreviations

$$\begin{eqnarray}&&\displaystyle \frac{{\chi }_{{BB}}}{{\rm{\Delta }}}\equiv \displaystyle \frac{{g}^{2}{a}_{S}}{2M},\end{eqnarray} \tag{ 47a }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\chi }_{{BF}}}{{\rm{\Delta }}}\equiv \displaystyle \frac{1}{{a}_{S}M\sqrt{{\ell }({\ell }+1)}},\end{eqnarray} \tag{ 47b }$$

and time is measured in units of M instead of Δ. We note that we have to take the limit ${\ell }\to \infty$ in order to recover the Hamiltonian formulation of lattice QED corresponding to (1).

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}={\psi }^{\dagger }{h}_{F}\psi$, where h_F is a matrix that 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

$$\begin{eqnarray}&&{\rho }_{W}(\psi ,{\varphi }_{1};{t}_{1})=\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\varphi }_{0}}{\pi }\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\eta }_{1}}{\pi }\exp \{{\eta }_{1}^{* }({\varphi }_{0}-{\varphi }_{1})-{\eta }_{1}({\varphi }_{0}^{* }-{\varphi }_{1}^{* })\}\\ &&\quad \times \exp \left\{-{\rm{i}}{\rm{\Delta }}{{tH}}_{W}\left(\psi ,{\varphi }_{1}+\tfrac{1}{2}{\eta }_{1}\right)\right\}{\rho }_{W}(\psi ,{\varphi }_{0};{t}_{0})\exp \left\{{\rm{i}}{\rm{\Delta }}{{tH}}_{W}\left(\psi ,{\varphi }_{1}-\tfrac{1}{2}{\eta }_{1}\right)\right\},\end{eqnarray} \tag{ 66 }$$

where we emphasize again that the Wigner transform ${H}_{W}(\psi ,\varphi )$ depends on both fermionic operators ψ and c-number variables φ. Iterating this expression, we obtain again a FI representation of the time evolution of the Wigner function

$$\begin{eqnarray}&&{\rho }_{W}(\psi ,{\varphi }_{N};{t}_{N})=\displaystyle \prod _{k=0}^{N-1}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\varphi }_{k}}{\pi }\displaystyle \prod _{k=1}^{N}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\eta }_{k}}{\pi }\exp \left\{\displaystyle \sum _{k=1}^{N}{\eta }_{k}^{* }({\varphi }_{k-1}-{\varphi }_{k})-{\eta }_{k}({\varphi }_{k-1}^{* }-{\varphi }_{k}^{* })\right\}\\ &&\quad \times T\exp \left\{-{\rm{i}}{\rm{\Delta }}t\displaystyle \sum _{k=1}^{N}{H}_{W}\left(\psi ,{\varphi }_{k}+\tfrac{1}{2}{\eta }_{k}\right)\right\}{\rho }_{W}(\psi ,{\varphi }_{0};{t}_{0})\bar{T}\exp \left\{{\rm{i}}{\rm{\Delta }}t\displaystyle \sum _{k=1}^{N}{H}_{W}\left(\psi ,{\varphi }_{k}-\tfrac{1}{2}{\eta }_{k}\right)\right\},\end{eqnarray} \tag{ 67 }$$

where we introduced the time ordering operator T and the anti-time ordering operator $\bar{T}$. 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

$$\begin{eqnarray}&&{H}_{W}(\psi ,{\varphi }_{k})={H}_{{BW}}({\varphi }_{k})+{H}_{{FW}}(\psi ,{\varphi }_{k}).\end{eqnarray} \tag{ 68 }$$

Furthermore, we assume that the initial density matrix factorizes into a purely bosonic part ${\rho }_{B}$ as well as a purely fermionic contribution ${\rho }_{F}$. Accordingly, the Wigner transform at initial times only affects the bosonic part, such that

$$\begin{eqnarray}&&{\rho }_{W}(\psi ,{\varphi }_{0};{t}_{0})={\rho }_{{BW}}({\varphi }_{0};{t}_{0}){\rho }_{F}(\psi ;{t}_{0})\end{eqnarray} \tag{ 69 }$$

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 ${\mathrm{Tr}}_{F}$, we obtain

$$\begin{eqnarray}&&{\mathrm{Tr}}_{F}{\rho }_{W}(\psi ,{\varphi }_{N};{t}_{N})=\displaystyle \prod _{k=0}^{N-1}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\varphi }_{k}}{\pi }\displaystyle \prod _{k=1}^{N}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\eta }_{k}}{\pi }{\rho }_{{BW}}({\varphi }_{0};{t}_{0})\\ &&\quad \times \exp \left\{\displaystyle \sum _{k=1}^{N}{\eta }_{k}^{* }({\varphi }_{k-1}-{\varphi }_{k})-{\eta }_{k}({\varphi }_{k-1}^{* }-{\varphi }_{k}^{* })-{\rm{i}}{\rm{\Delta }}t\displaystyle \sum _{k=1}^{N}\left[{H}_{{BW}}\left({\varphi }_{k}+\tfrac{1}{2}{\eta }_{k}\right)-{H}_{{BW}}\left({\varphi }_{k}-\tfrac{1}{2}{\eta }_{k}\right)\right]\right\}\\ &&\quad \times {\mathrm{Tr}}_{F}\left[T\exp \left\{-{\rm{i}}{\rm{\Delta }}t\displaystyle \sum _{k=1}^{N}{H}_{{FW}}\left(\psi ,{\varphi }_{k}+\tfrac{1}{2}{\eta }_{k}\right)\right\}\right.\left.{\rho }_{F}(\psi ,{t}_{0})\bar{T}\exp \left\{{\rm{i}}{\rm{\Delta }}t\displaystyle \sum _{k=1}^{N}{H}_{{FW}}\left(\psi ,{\varphi }_{k}-\tfrac{1}{2}{\eta }_{k}\right)\right\}\right].\end{eqnarray} \tag{ 70 }$$

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

$$\begin{eqnarray}&&{\psi }^{\dagger }{h}_{{FW}}(\varphi )\psi \equiv {H}_{{FW}}(\psi ,\varphi ).\end{eqnarray} \tag{ 71 }$$

Here, ${h}_{{FW}}(\varphi )$ is a matrix that only depends on bosonic variables but not on fermionic operators. Further, we introduce ${\varphi }_{k}^{\pm }$ as the linear combination of center field ${\varphi }_{k}$ and difference field ${\eta }_{k}$,

$$\begin{eqnarray}&&{\varphi }_{k}^{\pm }={\varphi }_{k}\pm \tfrac{1}{2}{\eta }_{k},\end{eqnarray} \tag{ 72 }$$

and we assume that the initial fermionic density matrix can be written as

$$\begin{eqnarray}&&{\rho }_{F}={{ \mathcal Z }}^{-1}\exp [{\psi }^{\dagger }{h}_{{FW}}({\varphi }_{0})\psi ],\end{eqnarray} \tag{ 73 }$$

where ${ \mathcal Z }$ is an appropriate normalization. By utilizing the identity [60]

$$\begin{eqnarray}&&{\mathrm{Tr}}_{F}[{{\rm{e}}}^{{\psi }^{\dagger }{M}_{1}\psi }\ldots {{\rm{e}}}^{{\psi }^{\dagger }{M}_{n}\psi }]=\det (1+{{\rm{e}}}^{{M}_{1}}\ldots {{\rm{e}}}^{{M}_{n}})\end{eqnarray} \tag{ 74 }$$

with matrices M_k for $k=1,\ldots ,n$, we may explicitly perform the fermionic trace in (70). Introducing the evolution matrix

$$\begin{eqnarray}&&{S}_{k,m}(\varphi )\equiv {{\rm{e}}}^{-{\rm{i}}{\rm{\Delta }}t{h}_{{FW}}({\varphi }_{k})}\ldots {{\rm{e}}}^{-{\rm{i}}{\rm{\Delta }}t{h}_{{FW}}({\varphi }_{m})},\end{eqnarray} \tag{ 75 }$$

where ${S}_{k,m}(\varphi )$ depends on the string of fields ${\varphi }_{k},\ldots ,{\varphi }_{m}$ for $k\geqslant m$, we obtain

$$\begin{eqnarray}&&{\mathrm{Tr}}_{F}[T\exp \{\ldots \}{\rho }_{F}(\psi ;{t}_{0})\bar{T}\exp \{\ldots \}]={{ \mathcal Z }}^{-1}\,\exp \,\mathrm{Tr}\,\mathrm{log}[1+{S}_{N,1}({\varphi }^{+}){{\rm{e}}}^{{h}_{{FW}}({\varphi }_{0})}{S}_{N,1}^{\dagger }({\varphi }^{-})]\end{eqnarray} \tag{ 76 }$$

for the last line of (70). For later convenience, we also define the fermionic propagator

$$\begin{eqnarray}&&{D}_{k}\equiv D(t+k{\rm{\Delta }}t)={S}_{k,1}(\varphi ){(1+{{\rm{e}}}^{-{h}_{{FW}}({\varphi }_{0})})}^{-1}{S}_{k,1}^{\dagger }(\varphi ),\end{eqnarray} \tag{ 77 }$$

whose time evolution at order ${ \mathcal O }({\rm{\Delta }}t)$ is governed by

$$\begin{eqnarray}&&{D}_{k+1}-{D}_{k}={\rm{i}}{\rm{\Delta }}t[{D}_{k},{h}_{{FW}}({\varphi }_{k+1})].\end{eqnarray} \tag{ 78 }$$

The components of D_k can be identified with the equal-time correlation function $\langle {\psi }_{m}^{\dagger }{\psi }_{n}\rangle$ as will be shown below. This discrete equation can be solved via a mode expansion of the operator ${\psi }_{m}$ as illustrated in [34]. Knowing the mode functions then allows one to compute fermionic correlation functions.

The expansion of the $\mathrm{Tr}\,\mathrm{log}$ in (76) up to linear order in the difference field ${\eta }_{k}$ then yields

$$\begin{eqnarray}&&{\mathrm{Tr}}_{F}{\rho }_{W}(\psi ,{\varphi }_{N};{t}_{N})=\displaystyle \prod _{k=0}^{N-1}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\varphi }_{k}}{\pi }\displaystyle \prod _{k=1}^{N}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\eta }_{k}}{\pi }{\rho }_{{BW}}({\varphi }_{0};{t}_{0})\exp \left\{\displaystyle \sum _{k=1}^{N}{\eta }_{k}^{* }({\varphi }_{k-1}-{\varphi }_{k})-{\eta }_{k}({\varphi }_{k-1}^{* }-{\varphi }_{k}^{* })\right\}\\ &&\quad \times \exp \left\{-{\rm{i}}{\rm{\Delta }}t\displaystyle \sum _{k=1}^{N}\left[\displaystyle \frac{\partial {H}_{{BW}}({\varphi }_{k})}{\partial {\varphi }_{k}}+\mathrm{Tr}\left(\displaystyle \frac{\partial {h}_{{FW}}({\varphi }_{k})}{\partial {\varphi }_{k}}{D}_{k}\right)\right]{\eta }_{k}\right\}\\ &&\quad \times \exp \left\{-{\rm{i}}{\rm{\Delta }}t\displaystyle \sum _{k=1}^{N}\left[\displaystyle \frac{\partial {H}_{{BW}}({\varphi }_{k})}{\partial {\varphi }^{* }}+\mathrm{Tr}\left(\displaystyle \frac{\partial {h}_{{FW}}({\varphi }_{k})}{\partial {\varphi }_{k}^{* }}{D}_{k}\right)\right]{\eta }_{k}^{* }\right\}.\end{eqnarray} \tag{ 79 }$$

Since the difference field ${\eta }_{k}$ appears only linearly in the exponent, we may integrate it out again. As compared to the purely bosonic case, see (60), the resulting complex Dirac-delta function now impose the discrete time evolution equation [61]

$$\begin{eqnarray}&&{\rm{i}}\displaystyle \frac{{\varphi }_{k}-{\varphi }_{k-1}}{{\rm{\Delta }}t}=\displaystyle \frac{\partial {H}_{W}({\varphi }_{k})}{\partial {\varphi }_{k}^{* }}+\mathrm{Tr}\left(\displaystyle \frac{\partial {h}_{{FW}}({\varphi }_{k})}{\partial {\varphi }_{k}^{* }}{D}_{k}\right),\end{eqnarray} \tag{ 80 }$$

which is implicitly solved via ${\varphi }_{k}=\tilde{g}({\varphi }_{k-1})$. As in the purely bosonic case, the Jacobi determinant does not contribute at leading order. Accordingly, the Wigner function ${\rho }_{W}({\varphi }_{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}={\psi }^{\dagger }A\psi$ with A being a matrix, we consider

$$\begin{eqnarray}&&\langle {O}_{F}(t)\rangle =\mathrm{Tr}\{{\psi }^{\dagger }A\psi \,\rho (t)\}=\int \displaystyle \frac{{{\rm{d}}}^{2}{\varphi }_{N}}{\pi }{\mathrm{Tr}}_{F}\{{\psi }^{\dagger }A\psi \,{\rho }_{W}(\psi ,{\varphi }_{N};t)\}.\end{eqnarray} \tag{ 81 }$$

While the bosonic contributions are identical to (70), the fermionic trace now takes the form

$$\begin{eqnarray}&&{\mathrm{Tr}}_{F}\left[{\psi }^{\dagger }A\psi \,T\exp \left\{-{\rm{i}}{\rm{\Delta }}t\displaystyle \sum _{k=1}^{N}{H}_{{FW}}\left(\psi ,{\varphi }_{k}+\tfrac{1}{2}{\eta }_{k}\right)\right\}\right.\left.{\rho }_{F}(\psi ;{t}_{0})\bar{T}\exp \left\{{\rm{i}}{\rm{\Delta }}t\displaystyle \sum _{k=1}^{N}{H}_{{FW}}\left(\psi ,{\varphi }_{k}-\tfrac{1}{2}{\eta }_{k}\right)\right\}\right].\end{eqnarray} \tag{ 82 }$$

Utilizing the identity

$$\begin{eqnarray}&&{\mathrm{Tr}}_{F}[{\psi }^{\dagger }A\psi \ {{\rm{e}}}^{{\psi }^{\dagger }{M}_{1}\psi }\ldots {{\rm{e}}}^{{\psi }^{\dagger }{M}_{n}\psi }]=\det (1+{{\rm{e}}}^{{M}_{1}}\ldots {{\rm{e}}}^{{M}_{n}})\mathrm{Tr}\{{(1+{{\rm{e}}}^{-{M}_{n}}\ldots {{\rm{e}}}^{-{M}_{1}})}^{-1}A\},\end{eqnarray} \tag{ 83 }$$

we may again perform the fermionic trace explicitly. Expanding the exponent up to linear order in the response field ${\eta }_{k}$ while setting it to zero otherwise, we finally obtain

$$\begin{eqnarray}&&\langle {O}_{F}({t}_{N})\rangle =\displaystyle \prod _{k=0}^{N}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\varphi }_{k}}{\pi }{\rho }_{{BW}}({\varphi }_{0};{t}_{0})\displaystyle \prod _{k=1}^{N}\delta ({\varphi }_{k}-\tilde{g}({\varphi }_{k-1}))\mathrm{Tr}\{{D}_{N}A\}.\end{eqnarray} \tag{ 84 }$$

Accordingly, the discrete time evolution of $\langle {O}_{F}\rangle$ is determined by

$$\begin{eqnarray}&&\langle {O}_{F}({t}_{k+1})\rangle -\langle {O}_{F}({t}_{k})\rangle ={\rm{i}}{\rm{\Delta }}t\mathrm{Tr}\{[{D}_{k},{h}_{{FW}}({\varphi }_{k+1})]A\}.\end{eqnarray} \tag{ 85 }$$

By explicitly introducing spatial indices and choosing ${A}_{{mn}}={\delta }_{{{mm}}_{1}}{\delta }_{{{nn}}_{1}}$, 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 $\langle {\psi }_{m}^{\dagger }{\psi }_{n}\rangle$.

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]

$$\begin{eqnarray}&&{H}_{{\rm{QL}},W}=\displaystyle \frac{{g}^{2}{a}_{S}}{4}\displaystyle \sum _{n}(| {b}_{n}{| }^{4}+| {d}_{n}{| }^{4})+M\displaystyle \sum _{n}{(-1)}^{n}{\psi }_{n}^{\dagger }{\psi }_{n}-\displaystyle \frac{{\rm{i}}}{2{a}_{S}\sqrt{{\ell }({\ell }+1)}}\displaystyle \sum _{n}({\psi }_{n}^{\dagger }{b}_{n}^{* }{d}_{n}{\psi }_{n+1}-{\rm{h}}.{\rm{c}}.),\end{eqnarray} \tag{ 86 }$$

where b_n and d_n are now c-numbers, but the ${\psi }_{n}$ and ${\psi }_{n}^{\dagger }$ 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 ${ \mathcal O }({\rm{\Delta }}{t}^{2})$. The equations of motion for the fermionic correlator D(t) are given by

$$\begin{eqnarray}&&{\rm{i}}{\partial }_{t}{D}_{m,n}(t)={h}_{{FW},{mk}}(t){D}_{k,n}(t)-{D}_{m,k}(t){h}_{{FW},{kn}}(t),\end{eqnarray} \tag{ 87 }$$

where the fermionic matrix h_FW(t) reads

$$\begin{eqnarray}&&{h}_{{FW},{mn}}=-\displaystyle \frac{{\rm{i}}}{2{a}_{S}\sqrt{{\ell }({\ell }+1)}}{b}_{m}^{* }{d}_{m}{\delta }_{m+1,n}+M{(-1)}^{m}{\delta }_{m,n}+\displaystyle \frac{{\rm{i}}}{2{a}_{S}\sqrt{{\ell }({\ell }+1)}}{d}_{m-1}^{* }{b}_{m-1}{\delta }_{m-1,n}.\end{eqnarray} \tag{ 88 }$$

The fermionic two-point function ${D}_{{mn}}=\langle {\psi }_{m}^{\dagger }{\psi }_{n}\rangle =\tfrac{1}{2}({\delta }_{{mn}}-{F}_{{nm}})$ can be expressed in terms of the statistical two-point function ${F}_{{mn}}=\langle [{\psi }_{m},{\psi }_{n}^{\dagger }]\rangle$. Accordingly, the equations of motion for the bosonic degrees of freedom are given by

$$\begin{eqnarray}&&{\rm{i}}{\partial }_{t}{b}_{n}=\displaystyle \frac{{g}^{2}{a}_{S}}{2}{b}_{n}^{* }{b}_{n}{b}_{n}+\displaystyle \frac{{{\rm{i}}{d}}_{n}}{4{a}_{S}\sqrt{{\ell }({\ell }+1)}}{F}_{n+1,n},\end{eqnarray} \tag{ 89a }$$

$$\begin{eqnarray}&&{\rm{i}}{\partial }_{t}{d}_{n}=\displaystyle \frac{{g}^{2}{a}_{S}}{2}{d}_{n}^{* }{d}_{n}{d}_{n}-\displaystyle \frac{{{\rm{i}}{b}}_{n}}{4{a}_{S}\sqrt{{\ell }({\ell }+1)}}{F}_{n,n+1}.\hspace{7pt}\end{eqnarray} \tag{ 89b }$$

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

$$\begin{eqnarray}&&{H}_{f,0}=-\displaystyle \frac{{\rm{i}}}{2{a}_{S}}\displaystyle \sum _{n}({\psi }_{n}^{\dagger }{\psi }_{n+1}-{\rm{h}}.{\rm{c}}.)+M\displaystyle \sum _{n}{(-1)}^{n}{\psi }_{n}^{\dagger }{\psi }_{n}.\end{eqnarray} \tag{ 90 }$$

Since ${H}_{f,0}$ is quadratic, the ground state and the dispersion relation can be determined analytically, see appendix C. Given an optical lattice with N elementary cells, i.e. $2N$ lattice sites, the dispersion relation is given by two bands $\pm {\omega }_{q}$ with

$$\begin{eqnarray}&&{\omega }_{q}=\sqrt{{M}^{2}+{\pi }_{q}^{2}},\quad {\pi }_{q}=\displaystyle \frac{1}{{a}_{S}}\sin \left(\displaystyle \frac{\pi q}{N}\right),\end{eqnarray} \tag{ 91 }$$

with $q\in \{0,\ldots ,N-1\}$. The corresponding mode function expansion of the fermionic field operator reads

$$\begin{eqnarray}&&{\psi }_{n}=\displaystyle \frac{1}{\sqrt{2N}}\displaystyle \sum _{q}{{\rm{e}}}^{\tfrac{{\rm{i}}\pi {qn}}{N}}\left(\displaystyle \frac{M+{(-1)}^{n}({\omega }_{q}-{\pi }_{q})}{\sqrt{2{\omega }_{q}({\omega }_{q}-{\pi }_{q})}}{a}_{q}\right.\left.-\displaystyle \frac{M-{(-1)}^{n}({\omega }_{q}+{\pi }_{q})}{\sqrt{2{\omega }_{q}({\omega }_{q}+{\pi }_{q})}}{c}_{q}^{\dagger }\right),\end{eqnarray} \tag{ 92 }$$

and the momentum space creation/annihilation operators are defined with respect to the fully filled lower band $| {\rm{vac}}\rangle$ according to ${a}_{q}| {\rm{vac}}\rangle ={c}_{q}| {\rm{vac}}\rangle =0$. The bosonic samples are prepared in an excited eigenstate of

$$\begin{eqnarray}&&{H}_{b,0}=\displaystyle \frac{{g}^{2}{a}_{S}}{2}\displaystyle \sum _{n}\,{L}_{z,n}^{2},\end{eqnarray} \tag{ 93 }$$

determined by the number of bosonic atoms on each site $2{\ell }={b}_{n}^{\dagger }{b}_{n}+{d}_{n}^{\dagger }{d}_{n}$ and the eigenvalue of the operator ${L}_{z,n}=({b}_{n}^{\dagger }{b}_{n}-{d}_{n}^{\dagger }{d}_{n})/2$. 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

$$\begin{eqnarray}&&{\rho }_{W}({\varphi }_{0})=\displaystyle \prod _{n}\delta ({b}_{n}-{{ \mathcal B }}_{0,n})\delta ({d}_{n}-{{ \mathcal D }}_{0,n}).\end{eqnarray} \tag{ 94 }$$

In fact, the corrections to the exact Wigner function are ${ \mathcal O }(1/{\ell })$ [57]. The explicit values of ${{ \mathcal B }}_{0,n}$ and ${{ \mathcal 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).

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–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 4.

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 ${\ell }=({N}_{b}+{N}_{d})/2$. 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}| =2{M}^{2}/{g}^{2}$. We first solve the equations of motion in the limit ${\ell }\to \infty$ for $g/M=0.1$, a_SM = 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 quasi-particle distribution n(k) from F_nm [32, 34, 67], see also appendix C. We display the time-evolution of the total particle number ${ \mathcal N }={\sum }_{k}n(k)$ in figure 5(a). 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]

$$\begin{eqnarray}&&\dot{{ \mathcal N }}=\displaystyle \frac{{M}^{2}E}{2\pi {E}_{c}}\exp \left(-\pi \displaystyle \frac{{E}_{c}}{E}\right).\end{eqnarray} \tag{ 95 }$$

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].

Zoom In Zoom Out Reset image size

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 ${ \mathcal O }(\delta \rho /{\ell })$ compared to QED it is interesting to investigate this error as a function of ℓ. In figure 5(a), we demonstrate the convergence of the cold atom behavior towards the QED result upon increasing ℓ. For the chosen parameters and ${\ell }=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 figure 5(b), we display the time-evolution of the electric field. Starting from QED (${\ell }\to \infty$), 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 figure 5(a).

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)=2{M}^{2}/{g}^{2}\gt 0$ until it changes sign and reaches a local minimum $I({t}_{\min })\lt 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 figure 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 ${\ell }={10}^{5}$ in figure 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 cold atom setup is found for the first oscillation period for N = 512, we still observe sizable deviation at later times.

Zoom In Zoom Out Reset image size

String breaking. The physics of confinement in the theory of 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–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 $\pm Q$ located at $\pm 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| \lt d/2$ whereas it vanishes outside of it.

We first make contact to the corresponding QED literature [33, 44] by considering the limit ${\ell }\to \infty$ and choosing $g/M=1$, a_SM = 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 figure 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 contain enough energy to produce the required charges $\pm Q$ to screen the external charges. In the latter case the string does not break completely.

Zoom In Zoom Out Reset image size

In figure 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}_{{\rm{CA}}}\lt {d}_{{\rm{QED}}}$ for the same parameters $g/M=1$, a_SM = 0.1 and N = 1024. As expected, we observe convergence towards the QED results upon increasing the value of ℓ.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{\bar{q}}_{n}={q}_{2n}+{q}_{2n+1}=-\displaystyle \frac{1}{2}({F}_{2n,2n}+{F}_{2n+\mathrm{1,2}n+1}).\end{eqnarray} \tag{ 96 }$$

In figure 9 we display the time evolution of the charge density ${\bar{q}}_{n}$. As described previously, the dynamical charges are produced on top of each other such that the total charge density vanishes initially.

Zoom In Zoom Out Reset image size

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].

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 3.4, the fermions can be prepared in the lowest band via a subsequential adiabatic ramp of ${V}_{2,\alpha }^{f}$ and then ${V}_{1,\alpha }^{f}$. 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 ⁶Li and blue detuned for ²³Na. 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 rf 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 4, 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 FI approach. This allows us to investigate the role of the experimentally relevant parameters on the physics of the Schwinger effect, especially the spin magnitude ${\ell }={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, ${\ell }\in \{10,20,\infty \}$. Since the experimental parameters allow for the exploration of the strong coupling regime $g/M\gt 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 figure 10.

Zoom In Zoom Out Reset image size

In figure 10(a) we display the time evolution of particle 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 ${\chi }_{{BF}}{N}_{B}$, which is short compared to the limiting particle losses, and thus accessible in the experiment. As to be expected from ${\chi }_{{BF}}{N}_{B}\gt {\rm{\Delta }}$ we also observe initial oscillations, which are smeared on times $t\gt h/{\rm{\Delta }}$. The particle number per lattice site then reaches a plateau with ${ \mathcal N }(t)/N\simeq {10}^{-3}$. For the given parameters, we find convergence towards the QED results already for ${\ell }=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, see also the simulations of dynamical string breaking in figure 8. Owing to the fact that realistic values of the spin magnitude are of the order of ${\ell }={ \mathcal 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 $I(t)={N}_{b}(t)-{N}_{d}(t)$ in figure 10(b). 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 figure 10. Nevertheless, the underlying physics of particle production can be already accessed from the integrated number of produced particles. According to figure 10, around 0.2 particles have to be detected on average for a lattice with N = 100 sites. Current experimental setups allow realizing $\approx 30$ copies of the cold atom QED system by employing an array of one-dimensional traps. Thus one expects integrated ${ \mathcal O }(10)$ particles which can be detected with fluorescence in a subsequent magneto-optical 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\in \{0.4,1.0,2.6\}$ for fixed parameters ${a}_{S}M=0.05$, N = 100 and ${\ell }\to \infty$ in figure 11. We note that we have to adapt the initial species imbalance $I(0)=2{M}^{2}/{g}^{2}$ in order to provide an initial critical field ${E}_{{\rm{c}}}$ 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.

Zoom In Zoom Out Reset image size