Dynamics of a Lattice Gauge Theory with Fermionic Matter -- Minimal Quantum Simulator with Time-Dependent Impurities in Ultracold Gases

We propose a minimal model to study the real-time dynamics of a $\mathbb{Z}_2$ lattice gauge theory coupled to fermionic matter in a cold atom quantum simulator setup. We show that dynamical correlators of the gauge fields can be measured in experiments studying the time-evolution of two pairs of impurities, and suggest the protocol for implementing the model in cold atom experiments. Further, we discuss a number of unexpected features found in the integrable limit of the model, as well as its extensions to a non-integrable case. A potential experimental implementation of our model in the latter regime would allow one to simulate strongly-interacting lattice gauge theories beyond current capabilities of classical computers.


INTRODUCTION
Due to remarkable experimental advances offering unprecedented control of isolated quantum systems, quantum simulators are becoming a reality, allowing one to test theoretical models of strongly-interacting quantum systems in the regimes beyond classical simulations. These recent experimental breakthroughs have been witnessed in a wide range of settings, including superconducting chips [1], photonic quantum circuits [2], and most notably in experiments with cold trapped ions [3], following early theoretical proposals, see e.g. Ref. [4]. Here we will focus on cold atom quantum simulator setups, such as recently used in the studies of many-body localization phenomena in two dimensions [5].
An intriguing and natural application of quantum simulators is in the studies of lattice gauge theories (LGT) [6][7][8][9][10], and in particular in testing their dynamical properties. Gauge theories play a central role in theoretical physics, from the standard model of fundamental particles to the low-energy descriptions of condensed matter systems. Important and wellknown examples are those which are used in theoretical models of quantum chromodynamics (QCD), quantum electrodynamics (QED), as well as quantum spin liquids [11], and Kitaev's toric code [12]. While models appearing in condensed matter theory context are often lattice models, QCD and QED are continuum theories. However, they have also been modelled using LGT approaches [13,14]. For example, a lattice version of the Schwinger model, which is a famous toy model of 1+1 dimensional QED [8,15], has been recently simulated in cold ion trap experiments [3].
Digital (discrete time) quantum simulations of LGT have been so far restricted to one-dimensional systems because of requirements on high fidelity and on the large number of qubits in implementation of these simulations. Here, we propose a minimal setting for simulating the dynamics of a LGT with fermions in 2+1 dimensions. We suggest that our model can be an ideal candidate for experimental implementations. First, we present a mapping of the model to free fermions (which can be studied exactly), and therefore the theory can be benchmarked against experiments. Via duality transformations we show that dynamical correlation functions of the gauge fields can be directly mapped to local impurity quenches in the free-fermion system. Second, even in the simple version of the model, where the free-fermion mapping holds, the model shows novel phenomenology of disorder-free localization [16][17][18]. Further, it can be tuned away from this 'integrable' limit, where classical computation is no longer applicable. Third, measurements of correlation functions can be implemented using current technology in cold atomic gases [5], and we propose a simple protocol based on Ramsey interferometry [19][20][21][22] in a cold atom setting.
We note that there has been a large number of papers related to quantum simulations of gauge theories, and we refer the reader to Refs. [9,10,23] for more detailed information.

Z2 LATTICE GAUGE THEORY WITH FERMIONIC MATTER
One of the goals in the field of quantum simulation is to be able to test models of interacting quantum field theories, including the ones used in QCD, and QED. In this paper we suggest a lattice description of a version of a QED Hamiltonian with gauge fields coupled to fermions, and we focus on a 2+1D case, although the model can be studied in any dimensions. The continuum version of the model reads (hereh,K are coupling constants, m the fermion mass, A the vector potential, and E, B are the electric and magnetic field strengths) (1) Note a non-standard (divE) 2 , in other words the energy depends only on the divergence of the electric field, but not on the field strength, see also a discussion in [24]. arXiv:1803.06575v3 [cond-mat.quant-gas] 28 Aug 2018 We discretise the model (1) by placing fermions on the sites of a 2D square lattice and introduce a discrete Z 2 vector potential and electric field on the links in the standard way, see e.g. [8][9][10]24], and arrive at the discrete Z 2 lattice gauge theory version of Eq. (1), where h, K, J are coupling constants,f † i are spinless fermion creation operators on sites i, σ z ij are the Pauli matrices defined on the links between neighbouring i and j sites. The starÂ i and plaquetteB p operators (which are well-known from the context of the toric code), which live on the sites i and plaquettes p are correspondingly defined aŝ The model Eq.
(2) is a natural two-dimensional generalization of the model studied by the authors in Ref. [18] in the context of disorder-free localization. The latter localization mechanism was later studied in the context of the 1D U (1) lattice Schwinger model in Ref. [25]. In this work we focus on dynamics of the gauge fields in the 2D case. One of the central results of this work is a protocol for quantum simulation of time-dependent gauge-field correlators after a quantum quench. We are interested in measuring the connected spin correlators

DUALITY MAPPING TO FREE FERMIONS AND DYNAMICAL CORRELATION FUNCTIONS
The model in Eq. (2) is the Kitaev Toric code with the additional term describing dynamics of free spinless fermions coupled to gauge fields via minimal coupling. Below we consider a quantum quench from an initial state |Ψ which is invariant under the application of all plaquette operators, i.e. B p |Ψ = |Ψ . Since these operators are conserved under dynamics we consider only the sector with eigenvalues B p = 1 for every p. More general quenches can be studied in a similar way, where one will have to take into account the disconnected sectors labelled by B p = ±1, see Refs. [18,24] for more details. Up to a constant, the Hamiltonian of Eq. (2) for a fixed sector takes the following form In order to bring the model into a solvable form we introduce a transformation of degrees of freedom. First, we perform a standard duality transformation for the spin degrees of freedom, which is well-known in the context of the Ising model,τ whereτ x i ,τ z i are the Pauli matrices defined on the sites i of the lattice. We can now identify conserved chargesq j = τ z j e iπf † jfj with the eigenvalues ±1, and introduce the following transformations for the fermion operatorsĉ j =τ x jf j . Thê τ and theĉ operators obey standard commutation relations. The Hamiltonian commutes with the operators of conserved charges [Ĥ,q j ] = 0, and the charges commute between themselves [q j ,q k ] = 0. In terms of theseĉ andτ operators, the Hamiltonian assumes the form (see [16][17][18]24]) The model in Eq. (7) is a free-fermion model, and can be studied using standard techniques, which we outline below on the example of a quantum quench problem. Note that the conserved charges in this case will be defined by the initial state of the system.

Quantum Quench Protocol
Below we focus on a quantum quench problem where the spins and the fermions are prepared in some initial state, and we will calculate the dynamics of spin-correlation functions at time t after the quench. An interesting, and simple initial state is given by a tensor product of the spins polarized along the z-axis and fermions in a Slater determinant state at half-filling |Ψ = | ↑↑ · · · ⊗ |ψ . The fermion Slater determinant describes a Fermi-sea for the HamiltonianĤ FS In terms of the chargesq and the fermionsĉ these initial states take the form where the primed sum is over all charge configurations ifi is the fermion filling, see Refs. [16][17][18] for more details.
Let us now discuss the calculation of the dynamical correlators. The simplest component of the connected spin-spin correlator is the average the local magnetisation In order to simplify notation we introduce a Hamiltonian for a fixed configuration of charges {q i } = ±1, together with the short hand notationĤ jk (q) for the Hamiltonian having the sign of q j , q k flipped relative toĤ(q), then the magnetization at time t after a quench is given by the following correlator  7). The fermions experience a binary potential set by charges {qi} = ±1, except at particular sites where we put spin-1/2 impurities (yellow). Impurity spins are localized in separately controlled much deeper wells, and generate an external potential on these sites, with the sign of the potential for spin up/down impurities being positive/negative correspondingly. To calculate spin correlators, corresponding to gauge-field correlators, one has to control four impurity spins, with pairs of spins located on the sites sharing the bond associated with the gauge spins. The correlators are then calculated using the Loschmidt echo protocol, defined in the main text, which involves a π/2-rotation and measurement of the impurity spins.
This has the form of binary disorder-averaged Loschmidt echo, where the chargesq at sites j, k having opposite signs in the forward and backward time-evolution. Repeating the same arguments we obtain the expression for the connected spin-correlator at time t after the quench where one has to exchange signs of four charges at sites denoted by jklm between forward and backward evolution. Eq. (12) is one of the central new results of this work. The dynamical correlation function of the gauge field directly correspond to a local quantum quench of a (free) fermionic Hamiltonian. In the following section we discuss how the latter can be efficiently simulated in a cold atom setup.
The Loschmidt echo appearing in these expressions can be efficiently computed numerically using fermion determinants. For example, for the expression Eq. (11) these take the form where U (q) = e −iH(q)t is the exponential of the matrix H(q), and similarly for U jk (q) and H jk (q), and V is a rectangular matrix which has as its columns the N/2 filled eigenvectors of the HamiltonianĤ FS = − ij f † if j . The determinants for the two-point correlators (12) take a similar form but with H jklm (q) instead of H jk (q).

QUANTUM SIMULATION AND EXPERIMENTAL SETUP
A schematic picture for the experimental setup that we propose is shown in Fig. 1. We consider a square optical lattice with fermions prepared in a Fermi-sea at half-filling. At time t = 0 we abruptly turn-on a binary disorder potential -similar to the quantum gas microscope set up used in Ref. [5] and we also add two or four impurity spins which control the potential in order to obtain the Loschmidt echo [20][21][22]. The impurities should be localized by an external trapping potential, and one has to chose these impurity spins such that the lattice fermions should interact strongly with one of the two spin states | ↑ j and weakly with the other | ↓ j , effectively turning on/off the potential for the fermions on the impurity site.
As shown above, the correlators that we are interested in can be mapped to the measurement of the Loschmidt echo. We prescribe an implementation of this measurement using Ramsey interferometry, which we describe below. Further details on experimental implementations can be found in [19][20][21][22]. The impurity spins in the up state interact with the fermions via a local interaction, while the spins in down states and the fermions are decoupled. Now we introduce a composite two state system, which we will call a control spin. In the measurement of the average local magnetisation the control spin affects two nearest-neighbour charges, i.e. | ⇓ ↔ | ↓ j ↓ k , | ⇑ ↔ | ↑ j ↑ k , and in the measurement of two point correlators we need two pairs of impurity spins (as shown in Fig. 1 Using the average local magnetisation as a specific example, the measurement protocol is the following: Initialise: The system is initialised in the state |ψ with fermions forming a half-filled Fermi-sea for the Hamil-tonianĤ FS = − ij f † if j , and the control spin is prepared in the state | ⇓ . π/2 pulse: We perform a π/2 pulse on the control spin such the state of the fermions and control spins becomes Evolve: We let the system evolve, so that its state at time t is given by π/2 pulse: At time t we perform the reverse π/2 pulse on the control spin so that we can perform a measurement in the original basis of impurity spins.

Measure:
We measure the z-component of the control spin, which gives This procedure must be repeated for different disorder realisations and the result averaged over these measurements. In the case of the two-point correlator the control spin corresponds to two pairs of impurities which amounts to replacingĤ jk (q) byĤ jklm (q). From our calculations we find that it is sufficient to average over a small subset of random binary disorder realisations. Because these spin correlators must be real, this procedure amounts precisely to the calculations we wish to perform in Eq. (4).

NUMERICAL RESULTS
Here we present some numerical results for a 15×14 square lattice for h/J = 0.7, 2 where we used averaging over 1000 disorder realisations. In Fig. 2 we show the time dependence of the connected two-point spin correlator (4) as a function of separation between two spins along the horizontal and diagonal cuts indicated by dashed lines in Fig. 3. Two main features common to all four figures are the linear light-cone spreading and the eventual decay of all spatial correlations.
The spreading of correlations is linear in all cases and has velocity v = 2J, which is the maximal group velocity of the fermions. This light-cone regime is short-lived due to the overall decay of the correlations. A notable difference between the horizontal and diagonal cuts is that the correlations between the neighbouring spins along the diagonal grow immediately, leading to a slightly offset light-cone. This can be appreciated by looking at Fig. 3 where one can see that these two spins belong to the same star operator and thus the correlations begin to grow after the quench with the rate set by h/J, as shown in Fig. 4. Difference appears when we increase h, where we see that the peaks in the correlations become much sharper along the light-cone, which is followed by decaying oscillations. The spatial pattern of correlations also changes as shown for a time slice Jt = 1.7 in Fig. 3. The extent of the spreading is bigger for smaller h/J. Further, for small values of h/J one can notice the asymmetry due to the fact that the central bond is vertical and we don't have 90 degree rotational symmetry, whereas for h = 2J this asymmetry seems to be smaller.
While for an exact simulation of the gauge field we would need to average over all possible configurations of the potential generated by the conserved charges, in order to obtain accurate results we require only a tiny fraction of configurations. For h/J = 0.7, we see that the results have the correct symmetry and there is very little noise (due to finite number of disorder realizations). On the other hand, for h/J = 2 there seem to be more-pronounced non-physical correlations, most notably appearing as a stripe in Figs. 2(c-d) at around Jt = 1, which is responsible for partially obscuring the light-cone. We can also see a faint non-uniform random background in Fig. 3(b). In order to minimise these artefacts one has to use a larger number of disorder realisations. Despite these issues, qualitatively good results can be obtained with as few as 50 disorder configurations as shown in Fig. 4 for the nearestneighbour correlators along diagonal cut. With a very small number of samples one can extract a number of qualitative   features such as sharp growth of correlations and the decaying oscillations.

DISCUSSION
We have presented a minimal 2D model of a Z 2 lattice gauge theory coupled to fermionic matter and outlined an experimental protocol for measuring time-dependent gauge field correlation functions. We believe that this protocol should be accessible with current experimental technology. In experiments with cold atoms in optical lattices, such as discussed in Ref. [5], relatively large 2D systems have already been simulated, and the protocol that we propose is designed to require minimal extra complications. Using a duality mapping to free-fermions we are able to perform efficient numerical computations which should allow one to quantitatively compare theory and experiment. Some clear features of dynamical correlations, such as the light-cone spreading and the decay of the correlations at long times should be directly accessible to experimental measurements.
We anticipate a number of possible practical issues regarding the difference between the experimental setup and our model, which should not be very difficult to account for. These may include an actual realization of the binary potential, the shape of a trapping potential for fermions, or realization of the coupling to impurity spins. All of these features can be easily implemented in our theoretical calculations. Further, one can access relatively large system sizes in numerical calculations, which would allow one to look at the scaling in the system size and the number of disorder realisations.
While the dynamics of the model in Eq. (5) can be efficiently simulated on a classical computer, the model offers a number of generalisations which turn it into a stronglyinteracting quantum model. Specifically, one can add to the Hamiltonian additional terms which will be still commuting with the charges. This means that the mapping and experimental protocol presented above are still valid. For example, one can add nearest-neighbour density interactions between fermions jk n jnk which have the same form in the original and the dual representation. In the 1D case the model maps onto an XXZ spin chain with binary disorder potential, realizing the minimal model of MBL behaviour. While there is no such mapping in 2D, the physics of MBL may be relevant in this case as well. One can also generalise our model using spin-1/2 fermions. In this case our model can be identified with the slave-spin description of the Hubbard model [26,27], but without Gauss' law imposed. With the addition of interactions, classical computations can only access very small system sizes in two dimensions (with some exceptions). However, in experiment, these generalisations are not expected to present significant extra difficulties.
To conclude, our model offers an ideal setting for simulating LGT dynamics in two-dimensions in the regimes beyond classical computation capabilities, with the important feature of a well defined free-fermion limit which provides theoretical results in a non-trivial setting.