Simulating Lattice Gauge Theories within Quantum Technologies

Lattice gauge theories, which originated from particle physics in the context of Quantum Chromodynamics (QCD), provide an important intellectual stimulus to further develop quantum information technologies. While one long-term goal is the reliable quantum simulation of currently intractable aspects of QCD itself, lattice gauge theories also play an important role in condensed matter physics and in quantum information science. In this way, lattice gauge theories provide both motivation and a framework for interdisciplinary research towards the development of special purpose digital and analog quantum simulators, and ultimately of scalable universal quantum computers. In this manuscript, recent results and new tools from a quantum science approach to study lattice gauge theories are reviewed. Two new complementary approaches are discussed: first, tensor network methods are presented - a classical simulation approach - applied to the study of lattice gauge theories together with some results on Abelian and non-Abelian lattice gauge theories. Then, recent proposals for the implementation of lattice gauge theory quantum simulators in different quantum hardware are reported, e.g., trapped ions, Rydberg atoms, and superconducting circuits. Finally, the first proof-of-principle trapped ions experimental quantum simulations of the Schwinger model are reviewed.


Introduction
In the last few decades, quantum information theory has been fast developing and consequently its application to the real world has spawned different technologies thatas for classical information theory -encompass the fields of communication, computation, sensing, and simulation [1,2,3]. To date, the technological readiness level of quantum technologies is highly diverse: some quantum communication protocols are ready for the market, while, e.g., universal quantum computers -despite experiencing an incredibly fast development -are still at the first development stage [4,5].
Some particularly interesting and potentially disruptive applications of quantum information theory and of quantum technologies lay within different scientific fields, such as high-energy, nuclear, condensed matter physics or chemistry [6]. Indeed, in the last years, it became increasingly clear that concepts and tools from quantum information can unveil new directions and will most probably provide new tools to attack long-standing open problems such as the study of information scrambling in black holes [7], the solution of complex chemical or nuclear systems [8], or the study of lattice gauge theories (LGTs) -the main subject of this review.
LGTs are characterised by an extensive number of symmetries, that is, conservation laws that hold at every lattice site. They describe an incredibly vast variety of different phenomena that range from the fundamental interactions of matter at high energies [9] -the standard model of particle physics -to the low-energy behaviour of some materials with normal and/or topological order in condensed matter physics [10,11]. Moreover, recently it has been shown that most of the hard problems in computer science can be recast as a LGT [12,13]. The connection passes through the recasting of the classical problem in Hamiltonian form, which generally assumes the form of an Ising Hamiltonian with long-range disordered interactions. This class of Hamiltonians can be mapped exactly into two-dimensional LGT [14].
For all the aforementioned scenarios, quantum science provided two novel pathways to analyse them. The first one has its root in Feynman's first intuition [15] of quantum computers: having quantum hardware able to precisely reproduce another physical quantum model, allows a powerful investigation tool for computing the observables of the model, and to verify or compare its prediction with the physical system. Today, the research frontier is at the edge of having universal quantum computers and quantum simulators able to perform such investigations beyond proof of principle analysis. Thus, detailed studies and proposals have been put forward to perform quantum simulations of LGT in the near and mid-term [16]. The second pathway exploits a class of numerical methods -tensor network methods (TNM) -which have been developed in the condensed matter and quantum information coma simone.montangero@unipd.it b enrique.rico.ortega@gmail.com munities to study strongly correlated many-body quantum systems [17]. Indeed, as it has been shown recently, TNM can be exploited to study LGT going in regimes where standard approaches are severely limited [18,19].
Lattice gauge theory was originally constructed by Wilson in order to define Quantum Chromodynamics (QCD) -the relativistic SU (3) gauge field theory that describes the strong interaction between quarks and gluons -beyond perturbation theory. For this purpose, he introduced a hyper-cubic space-time lattice as a gauge invariant regulator of ultraviolet divergences, with quark fields residing on lattice sites and gluons fields residing on links connecting nearest-neighbour sites. This framework makes numerous important physical quantities accessible to first principles Monte Carlo simulations using classical computers. These include static properties, like masses and matrix elements, of baryons (such as protons and neutrons) and mesons (such as pions). Properties of the hightemperature quark-gluon plasma in thermal equilibrium are accessible as well. This includes, e.g., the critical temperature of the phase transition in which the chiral symmetry of the quarks, which is spontaneously broken at low temperatures, gets restored. After more than four decades of intensive research, lattice QCD has matured to a very solid quantitative tool that is indispensable for correctly interpreting a large variety of experiments, including the ones at the high-energy frontier of the Large Hadron Collider (LHC) at CERN.
However, there are other important aspects of the QCD dynamics, both at high baryon density (such as in the core of neutron stars) and for out-of-equilibrium real-time evolution (such as the various stages of heavy-ion collisions), where importance-sampling-based Monte Carlo simulations fail due to very severe sign or complex action problems. In these cases, reliable special purpose quantum simulators or universal quantum computers may be the only tools to successfully address these grand-challenge problems. While immediate results with quantitative impact on particle physics are unrealistic to hope for, a longterm investment in the exploration of quantum technologies seems both timely and most interesting. Lattice gauge theory has a very important role to play in this endeavour, because, besides fully-fledged lattice QCD, it provides a large class of simpler models, in lower dimensions, with simpler Abelian or non-Abelian gauge groups, or with a modified matter content, which often are interesting also from a condensed matter perspective. The real-time evolution of all these models is as inaccessible to classical simulation as the real-time evolution of QCD itself. Hence, learning how to tackle with these challenges in simpler models is a necessary and very promising step towards the ultimate long-term goal of quantum simulating QCD. Along the way, via a large variety of lattice field theory models, particle physics provides an important intellectual stimulus for the development of quantum information technology.
Validation of quantum simulation experiments is vital for obtaining reliable results. In certain cases, which are limited to equilibrium situations, importance sampling Monte Carlo simulations using classical computers can provide such validation. However, Matrix Product States (MPS) and Tensor Network (TN) calculations are often the more promising method of choice, in particular, because they can even work in some out-of-equilibrium realtime situations. This provides an important stimulus to further develop these techniques. While they work best in one (and sometimes in two) spatial dimensions, an extension to higher dimensions is not at all straightforward, but very well worth to pursue vigorously. Even if these methods should remain limited to lower dimensions, they offer a unique opportunity to gain a deep understanding of the real-time evolution of simple lattice gauge models. By quantitatively validating quantum simulators in outof-equilibrium situations, even if only in lower dimensions, MPS and TN methods play a very important role towards establishing quantum simulators as reliable tools in quantum physics.
This paper reviews the recent activities along these lines, in particular of the groups that form the QTFLAG consortium, a European project funded under QuantERA with the goal of developing novel quantum science approaches to simulate LGT and study physical processes beyond what could be done via standard tools. First, the main concepts of interest are introduced, the LGT formulation and the tools used to study them: quantum simulators on different hardware and tensor network methods. Then, the recent numerical studies of one-dimensional Abelian and non-Abelian LGTs in and out of equilibrium, at zero and finite temperature are presented. Different theoretical proposals for the implementation of LGTs on digital and analog quantum simulators in trapped ions, Rydberg atoms and in superconducting circuits are reviewed. Finally, the first experimental realisations of these ideas are also briefly mentioned.

Lattice Field Theory background
Gauge field theories are at the heart of the current theoretical understanding of fundamental processes in nature, both in condensed matter and in high-energy physics. Although their formulation appears to be simple, they potentially give rise to very intriguing phenomena, such as asymptotic scaling, confinement, spontaneous chiral symmetry breaking or (non-trivial) topological properties, which shape the observed physical world around us. Solving gauge theories from first principles has been a major goal for several decades. Their formulation on a discrete Euclidean space-time lattice, originally proposed by Wilson in the seventies [20], has provided very powerful methods to study the non-perturbative regimes of quantum field theories 1 . A most prominent example is the success of ab-initio Lattice Quantum Chromodynamics (LQCD) calculations. Here, starting from the QCD Lagrangian, the low-lying baryon spectrum could be computed on very large lattices 1 See, however, [21] for a recent alternative approach using a ζ-regularisation and extrapolated to the continuum limit [22]. Lattice QCD calculations have also provided most important insights into the structure of hadrons [23,24]; they provide information on non-perturbative contributions to electroweak processes [25] and flavour physics [26]; and they are very successful to determine thermodynamic properties [27]. Today, lattice calculations are performed on large lattices -presently of sizes around 100 3 × 200 lattice points -and directly in physical conditions. By controlling systematic errors, such as discretisation and finite-volume effects, lattice field theory, and in particular lattice QCD, is providing most important input to interpret and guide ongoing and planned experiments world-wide, such as those at the Large Hadron Collider at CERN. These most impressive results became possible by a combined progress on algorithmic and computational improvements as well as the development of new supercomputer architectures. Thus, lattice field theory computations have demonstrated the potential to characterise the most fundamental phenomena observed in nature.
The standard approach of lattice field theory relies on Monte Carlo-based evaluations of path integrals in Euclidean space-time with positive integrands. Thus it suffers from an essential limitation in scenarios that give rise to a sign problem. These include the presence of a finite baryon density, which is relevant for the early universe and for neutron stars; real-time evolution, e.g., to understand the dynamics of heavy-ion collisions; or topological terms, which could shed light on the matter-anti-matter asymmetry of the universe. There is therefore an urgent quest to find alternative methods and strategies that enable tackling these fundamental open problems in the understanding of nature.
One such alternative is the application of tensor networks (TN). Originally introduced in the context of condensed matter physics, TN can solve quasi-exactly one dimensional strongly correlated quantum many-body problems for system sizes much larger than exact diagonalisation allows. They are naturally free from the sign problem. In fact, for 1-dimensional systems a number of successful studies have demonstrated the power of TN for lattice gauge theory calculations [28]. In particular, it has been shown that TN provide accurate determinations of mass spectra and that they can map out a broad temperature region. They can also treat chemical potentials and topological terms and they can be used to study realtime dynamics. TN also allow the study of entanglement properties and the entropy (leading in turn to the determination of central charges) in gauge theories, which brings new aspects of gauge theories into focus. However, applications to higher-dimensional problems remain a challenge presently. There are well-founded theoretical formulations such as projected entangled pair states (PEPS) but their practical application is still rather limited (for a recent review see [29]). New ideas such as the ones developed in [30] could have the potential to overcome these limitations but clearly further studies and developments are necessary in order to turn them into practical tools for addressing gauge theories in higher dimensions.
Ultimately, the intrinsic quantum nature of lattice gauge theories will be a limiting factor for classical calculations, even for TN, e.g., when out-of-equilibrium phenomena of a system are to be studied. In this context, quantum simulation, i.e., the use of another well-controlled quantum system to simulate the physics of the model under study, appears as a more adequate strategy. The idea of quantum simulation, first proposed by Feynman [15], is now becoming a reality [31,32,33,34], and very different condensed matter models have already been successfully quantum simulated in cold-atom laboratories around the world [35,36,37]. Regarding the simulation of LGT, a number of proposals have been put forward in the last years [38,39,40,41,42,43,44,45,46], and have even been realised by a few pioneering experiments [47]. TN calculations have been instrumental in the definition of many of these proposals. It is in particular the approach of hybrid quantum-classical simulation schemes which can take advantage of these new concepts and there is a great potential to realise them on near-term quantum architectures.
Gauge fields on the lattice manifest themselves as parallel transporters residing on the links that connect neighbouring lattice sites. In Wilson's formulation of lattice gauge theory, the link parallel transporters take values in the gauge group [20]. As a consequence, for continuous gauge groups such as the Abelian U (1) gauge group of Quantum Electrodynamics (QED) or the non-Abelian SU (3) gauge group of QCD, the link Hilbert space is infinite dimensional. When gauge fields are treated by TN techniques or they are embodied by ultra-cold matter or quantum circuits, representing an infinite dimensional Hilbert space is challenging, because usually only a few quantum states can be sufficiently well controlled in quantum simulation experiments. There are different approaches to addressing this challenge. First, the link Hilbert space of the Wilson theory can be truncated to a finite dimension in a gauge-covariant manner. Gradually removing the truncation by a modest amount allows one to take the continuum limit.
An alternative approach is provided by quantum link models (also known as gauge magnets) [48,49,50] which work with quantum degrees of freedom with a finite-dimensional Hilbert space from the outset. For example, the parallel transporters of a U (1) quantum link model are constructed with quantum spins [51], which can naturally be embodied by ultra-cold matter. Again, when one moderately increases the spin value, one can reach the continuum limit. Both approaches are actively followed presently and it will be interesting to see in the future, which strategy will be most appropriate to treat gauge theories with TN or on quantum devices.
Even when one restricts oneself to the smallest spin value 1 2 , interesting gauge theories emerge. For example, when its Gauss' law is appropriately modified, the U (1) quantum link model turns into a quantum dimer model [52,53], which is used in condensed matter physics to model systems related to high-temperature superconductors. Kitaev's toric code [54] -a topologically protected storage device for quantum information -provides an example of a Z(2) lattice gauge theory formulated with parallel transporters consisting of quantum spins 1 2 . Quantum spin chains were among the first systems to be quantum simulated successfully. SU (N ) quantum spin ladders, i.e., systems consisting of n transversely coupled spin chains, can be quantum simulated with ultra-cold alkaline-earth atoms in optical lattices [55]. For moderate values of n, these (2 + 1)-dimensional systems dimensionally reduce to (1 + 1)-dimensional CP (N − 1) models, which are asymptotically free and thus serve as toy models for QCD. Furthermore, for odd n they have non-trivial topology, very much like non-Abelian gauge theories in four space-time dimensions. Also QCD itself can be formulated as a quantum link model [56,57]. In that case, the parallel transporters are matrices with non-commuting matrix elements, just as quantum spins are vectors with non-commuting components. Alkaline-earth atoms can again be used to encode the QCD color degree of freedom in the nuclear spin of these atoms [43]. Lattice gauge theory, either in its gauge covariantly truncated Wilson formulation or in the description of quantum link models, which nicely complement each other, provides a broad framework for upcoming quantum simulation experiments.
Whatever the most effective simulations may be in the future, classical Monte Carlo, tensor network, or quantum simulations for addressing gauge theories, there will remain a big challenge: in the end, all calculations aim at providing input for world-wide experiments, whether the ones in condensed matter physics or the large-scale collider experiments in high-energy physics. As a consequence, all results emerging from theoretical computations based on the underlying Hamiltonian or Lagrangian need to have controlled statistical and systematic errors. This will lead to a substantial, demanding effort for such calculations since many simulations at various values of the lattice spacing and lattice volumes as well as possibly other technical parameters (e.g., the bond dimension in the TN approach) have to be executed. Only by performing a controlled continuum and infinite volume (or infinite bond dimension) limit, it will become possible to rigorously attribute the obtained results to the underlying model. In this way, the underlying model can be thoroughly tested and, in turn, any significant deviation seen in experiment will thus open the door to completely new and unexplored physics.

Quantum Science and Technologies tools
In a seminal paper published in 1982, Feynman discussed in great detail the problems connected with the numerical simulation of quantum systems. He envisaged a possible solution, the so-called universal quantum simulator, a quantum-mechanical version of the usual simulators and computers currently exploited in many applications of the 'classical' world. If realised, such a device would be able to tackle many-body problems with local interactions by using the quantum properties of nature itself. Interest-ingly, even without the advent of a fully universal quantum computer, the construction of dedicated devices, also known as purpose-based quantum simulators, would already be of significant importance for the understanding of quantum physics. The basic idea is to engineer the Hamiltonian of the quantum model of interest in a highly controllable quantum system and to retrieve all of the desired information with repeated measurements of its properties. Many research fields would eventually benefit from such devices: for example, two-dimensional and three-dimensional many-body physics, non-equilibrium dynamics or lattice gauge theories.
In recent years, the scientific community has been considering several quantum technologies such as cold atoms, trapped ions or superconducting circuits as examples of the most promising candidates for the realisation of a wide variety of dedicated quantum simulations. Indeed, these platforms are genuine quantum systems where the available experimental techniques offer an impressive degree of control together with high-fidelity measurements, thus combining two fundamental requirements for a quantum simulator. Among the most recent experimental achievements, just to mention a few, such as, the observation of Anderson localisation in disordered Bose-Einstein condensates (BECs), the research on itinerant ferromagnetism with cold fermions or the reconstruction of the equation of state of fermionic matter in extreme conditions, such as in neutron stars.
The advantages of quantum simulation are numerous: first, one can use it to study physical systems which are not experimentally accessible (systems of large or small scales, for example), or to observe the physical properties of unreal physical systems, which are not known to be found in nature, but can be mapped to the simulating systems. So far, a lot of quantum simulations were suggested, and some were even experimentally implemented. The simulated systems come from almost every area of physics: condensed matter and relativistic quantum physics, gravity and general relativity, and even particle physics and quantum field theory. The last of these is the topic of this review, specifically gauge theories. While quantum simulations have been proposed (and even realised) for condensed matter models, gauge theories are a newer branch where quantum technologies might be employed. See also [58,59,60,61,62,63,64,65,36] and for a general review on quantum simulation [6].

Quantum information techniques 4.1 Tensor networks for Lattice Gauge Theories
As mentioned before, classical numerical simulations are playing a leading role in the understanding of lattice gauge theories. In particular, in recent years, there has been a boost in the development of tensor network methods to simulate lattice gauge theories. There are different approaches, that range from the exploitation of mappings of some theories to spin models [66,67], to the development of gauge invariant tensor networks in the quantum link formulation [68,18,69,70,71,19]. This section reviews some of the studies that appeared in the last years, covering most of the available approaches for Abelian and non-Abelian lattice gauge theories [66,67,68,72,73,74,75].
In the following subsections, a selection of works performed along these lines is described in some detail.

Matrix Product States for Lattice Field
Theories [66,67] The Schwinger model [76,77], i.e., QED in one spatial dimension, is arguably the simplest theory of gauge-matter interaction, and yet it exhibits features in common with more complex models (like QCD) such as confinement or a non-trivial vacuum. Therefore, it constitutes a fundamental benchmark to explore the performance of lattice gauge theory techniques. In particular it has been extensively used in the last years to probe the power of TN as alternative methods to conventional Monte Carlo-based lattice techniques for solving quantum field theories in the continuum.
The first such study was carried out by Byrnes and coworkers [78] using the original DMRG formulation, and it already improved by orders of magnitude the precision of the ground state energy and vector particle mass gap, with respect to results obtained by other numerical techniques, although the precision decreased fast for higher excitations. The application of TN formulated algorithms, including extensions to excited states, time evolution and finite temperature has allowed a more systematic exploration of the model in recent years.
The discretised Hamiltonian of the model, in the Kogut-Susskind formulation with staggered fermions [79] reads where φ † n represents the creation operator of a spin-less fermion on lattice site n, and U n = e iθn is the link operator between sites n and n + 1. L n , canonical conjugate to θ n , corresponds to the electric field on the link, and α corresponds to a background field. Physical states need to satisfy Gauss' law as an additional constraint, In the continuum, the only dimensionless parameter of the model is the fermion mass m/g (expressed in terms of the coupling). The discretisation introduces one more parameter, namely the lattice spacing ag. For convenience, the Hamiltonian is often rescaled and expressed in terms of the dimensionless parameters x = 1/(ag) 2 , µ = 2 √ xm/g, with the continuum limit corresponding to x → ∞. The local Hilbert space basis for the fermionic sites can be labeled by the occupation of the mode, φ † n φ n ∈ {0, 1} (for site n), while the basis elements for the links can be labeled by the integer eigenvalues of L n , n . Using this basis, an MPS ansatz can be optimised to approximate the ground state or the excitations.  Instead of working with explicit gauge degrees of freedom, it is possible to integrate them out using Gauss' law, and to work directly in the physical subspace. This results in a Hamiltonian expressed only in terms of fermionic operators, but with non-local interactions among them. Additionally, a Jordan-Wigner transformation can be applied to map the model onto a more convenient spin Hamiltonian [80]. In [67] a systematic study of the mass spectrum in the continuum was performed using MPS with open boundary conditions, in the absence of a background field, for different values of the fermion mass. The ground state and excitations of the discrete model were approximated by MPS using a variational algorithm, and the results were successively extrapolated in bond dimension, system size (individual calculations were done on finite systems) and lattice spacing, in order to extract the continuum values of the ground state energy density and the mass gaps ( Fig. 1 illustrates the continuum limit extrapolations). These steps resemble those of more usual lattice calculations, so that also standard error analysis techniques could be used to perform the limits and estimate errors, and thus gauge the accuracy of the method. Values of the lattice spacing much smaller than the usual ones in similar Monte Carlo calculations could be explored, and very precise results were obtained for the first and second particles in the spectrum (respectively vector and scalar), beyond the accuracy of earlier numerical studies (see table  1).
Since the algorithms provide a complete ansatz for each excited state, other observables can be calculated.  An interesting quantity is the chiral condensate, order parameter of the chiral symmetry breaking, and written in the continuum as Σ = Ψ (x)Ψ (x) /g. When computed on the lattice, the condensate has a UV divergence, which is already present in the free theory. Using the MPS approximations for the ground state, the continuum limit of the condensate was extracted in [66] (some of these results were refined later in [81]). After subtracting the UV divergence, lattice effects were found to be dominated by corrections of the form a log a. Systematic fitting and error analysis techniques were applied to obtain very precise estimations of the condensate for massless and massive fermions (table 2, see also results with uniform MPS [82] and infinite DMRG [83]). In the former case the exact value can be computed analytically, but for the latter, very few numerical estimations existed in the literature. These results demonstrate the feasibility of the MPS ansatz to efficiently find and describe the low-energy part of the spectrum of a LGT in a non-perturbative manner. Moreover they show explicitly how the errors can be systematically controlled and estimated, something fundamental for the predictive power of the method, if it is to be used on theories for which no comparison to an exact limit is possible.
In this section, the general systematics of this approach is reviewed vis-à-vis the particularities that come with the simulation of gauge field theories in the continuum limit. An overview of the most important results that result from these simulations are also shown.
Continuum limit. As in the approach of both Byrnes et al. [78,90] and Bañuls et al. [67], the simulations start from a discretisation of the QFT Hamiltonian with the Kogut-Susskind prescription [79] followed by a Jordan-Wigner transformation. But different from [78,90,67], the simulations [73,82,84,85,86,87,88] are performed directly in the thermodynamic limit, avoiding the issue of finite-size scaling. From the lattice point of view, the QFT limit is then reached by simulating the model near (but not at) the continuum critical point [91]. Upon approaching this critical point the correlation length in lattice units diverges ξ/a → ∞. Large scale correlations require more real-space entanglement, specifically for the Schwinger model the continuum critical point is the free Dirac-fermion c = 1 CFT, implying that the bipartite entanglement entropy should have a UV-divergent scaling of 1/6 ln(ξ/a) [92]. This was confirmed explicitly by the numerical MPS simulations of the ground state of the Schwinger model [73,84], as shown in figure 3a. Notice the same UV scaling for ground states in the presence of an electric background field Q, leading to a UV finite subtracted entropy (see the inset), that can be used as a probe of the QFT IR physics [84].
For the MPS simulations this UV divergence of the entanglement requires bond dimensions D that grow polynomially with the inverse lattice spacing, D ∼ a −n . But it turns out that, despite this polynomial growth one can simulate the Schwinger model sufficiently close to its con-  [84].
tinuum critical point a → 0, with a relatively low computational cost. In the different papers simulations were performed up to a ≈ 1/(30g), corresponding to a correlation length ξ/a ≈ 15 − 35, depending on the particular ratio of the fermion mass and gauge coupling m/g in the Hamiltonian. The simulations at different decreasing values of a then allow for very precise continuum extrapolations, as illustrated in figure 3b for the dispersion relation of the (lowest lying) excitation [73]. Notice that in contrast to e.g. d = 3 + 1 QCD, d = 1 + 1 QED is a superrenormalizable theory, with a finite continuum extrapolation of the particle excitation masses m phys in terms of the bare parameters (m, g) of the theory: m phys (a) = m phys (0) + O(a). A further study demonstrated that even simulations with lattice spacings a ≥ 1/(10g) (implying a smaller computational cost) are already sufficient for continuum extrapolations with four digit precision [88].
Truncating the gauge field. The numerical Hamiltonian MPS simulations require finite local Hilbert spaces, which is in apparent conflict with the bosonic gauge degrees of freedom that come with the continuous U (1) group of the Schwinger model. As became evident in the work of Buyens et al., these bosonic fields can be efficiently truncated in the electric field basis, leading to an effective finite local Hilbert space appropriate for the simulations. In figure 4a the distribution of the Schmidt values 2 is shown over the different electric field eigenvalues q for a particular ground state simulation. Notice that the electric field values are discrete in the compact QED formulation. As one can see from the figure, the contribution from the 2 By the singular value decomposition, any matrix M can be decomposed in a positive semidefinite diagonal matrix D and two unitaries matrices U and V such that M = U DV † . The diagonal elements of the matrix D are called the Schmidt values or Schmidt spectrum higher electric field values decays rapidly, in fact exponentially, and it was shown that this decay remains stable towards the continuum limit [88]. For a given Schmidt precision one can therefore indeed safely truncate in q. Most simulations used q ∈ [−3, 3].
Gauge invariance. As was discussed already in previous sections, the Kogut-Susskind set-up starts from the Hamiltonian QFT formulation in the time-like axial gauge A 0 = 0, with the physical states obeying the Gauss constraint ∇E = ρ. This is indeed equivalent to requiring the physical states to be invariant under local gauge transformations. The resulting lattice Hamiltonian then operates on a Hilbert space of which only a subspace of gauge invariant states, obeying a discretised version of Gauss' law, is actually physical. The simulations of Buyens et al. exploited this gauge invariance by constructing general gauge invariant MPS states [73] and simulating directly on the corresponding gauge invariant manifold. As shown in [82], for ground state simulations, working with explicit gauge invariant states leads to a considerable reduction in the computation time. The reason lies in the sparseness of the matrices appearing in gauge invariant MPS states; but also in the fact that the full gauge variant Hilbert space contains pairwise excitations of non-dynamical point charges, separated by short electric field strings of length L ∼ a. In the continuum limit this leads to a gapless spectrum for the full Hilbert space, whereas the spectrum on the gauge invariant subspace remains gapped. Such a nearly gapless spectrum requires many more time steps before convergence of the imaginary time evolution towards the proper ground state. As such, these test simulations on the full Hilbert space [82] are consistent with Elitzur's theorem [93], which states that a local gauge symmetry cannot be spontaneously broken, ensuring the same gauge invariant ground state on the full gauge variant Hilbert space.
Results. Using the Schwinger model [77,94] as a very nice benchmark model for numerical QFT simulations, the results of the numerical simulations [73,82,84,85,86,87,88] were verified successfully against these analytic QFT results in the appropriate regimes. In addition, where possible, the results were compared with the numerical work of [78,90,67], and found to be in perfect agreement within the numerical precision. Taken together, the tensor network simulations of Byrnes, Bañuls, Buyens et al., form the current state of art of numerical results on the Schwinger model. Now, a selection of the results of Buyens et al. are discussed: Ground state and particle excitations. By simulating the ground state and constructing ansatz states on top of the ground state, MPS techniques allow for an explicit determination of the approximate states corresponding to the particle excitations of the theory [95]. For the Schwinger model three particles were found [73]: two vector particles (with a quantum number C = −1 under charge conjugation) and one scalar particle (C = +1). For each of these particles, the obtained dispersion relation is perfectly consistent with an effective Lorentz symmetry at small momenta, as illustrated in Fig. 3b Table 3: Energy density and masses of the one-particle excitations (in units g = 1) for different m/g. The last column displays the result for the heavy vector boson [73].
tor excitation was uncovered for the first time, confirming prior expectations from strong coupling perturbation theory [77,94]. See the extrapolated mass values obtained for the scalar and first vector particle in absence of a background field in Table 3. Furthermore, in [85] the excitations were studied in presence of a background electric field. By extrapolating towards a vanishing mass gap for a half-integer background field, this allowed for a precise determination of the critical point (m/g) c = 0.3308 in the phase diagram [88]. String breaking. By probing the vacuum of a confining theory with a heavy charge/anti-charge pair, one can investigate the detailed physics of string formation and breaking, going from small inter-charge distances to larger distances. In the latter case the heavy charges get screened by the light charged particles that are created out of the vacuum. This string breaking picture was studied in detail for the Schwinger model in [84]. Figure 4b shows one of the results on the light particle charge density for different distances between the heavy charges. At small distances there is only a partial screening, whereas at large distances the screening is complete: for both fully integrated clouds, the total charge is exactly ±1. For large values of Lg, the string is completely broken and the ground state is described by two free particles, i.e. mesons. Notice the red line in the plot which depicts the corresponding analytic result of the ground state charge distribution for the nonrelativistic hydrogen atom in d = 1 + 1. Finally, also fractional charges were studied in [84], explicitly showing for the first time the phenomenon of partial string breaking in the Schwinger model.

Tensor networks for Lattice Gauge Theories and Atomic Quantum Simulation[74]
In [74], an exact representation of gauge invariance of quantum link models, Abelian and non-Abelian, was given in terms of a tensor network description. The starting point for the discussion are LGTs in the Hamiltonian formulation, where gauge degrees of freedom U x,y are defined on links of a lattice, and are coupled to the matter ones ψ x , defined on the vertices. In the quantum link formulation, the gauge degrees of freedom are described by bilinear operators (Schwinger representation). This feature allows one to solve exactly, within the tensor network represen- More concretely, in a fermionic Schwinger representation of a non-Abelian U (N ) quantum link model, the gauge operators U ij x,y that live on the links x, y of a ddimensional lattice, with color indices i, j are expressed as a bilinear of fermionic operators, U ij x,y = c i x c j † y . In this link representation, the number of fermions per link is a con- In models with matter, at every vertex x of the lattice, there is a set of fermionic modes ψ i x with color index i. The left and right generators of the SU (N ) symmetry are defined as L a x,y = i,j c i † x λ a i,j c j x and R a x,y = i,j c i † y λ a i,j c j y , with λ a i,j the group structure constants. Hence, the non-Abelian generators of the gauge symmetry are given by x , withk the different directions in the lattice. There are also similar expressions for the Abelian part of the group G x .
The "physical" Hilbert subspace is defined as the one that is annihilated by every generator, i.e., G x |phys = G a x |phys = 0 ∀x, a. A particular feature of quantum link models is that, these operators being of bosonic nature (they are bilinear combinations of fermionic operators), the spatial overlap between operators at different vertices x or y is zero, i.e., G a x G b y = 0, ∀a, b and x = y, even between nearest-neighbours. In this way, (i) the gauge invariant Hilbert space (or Gauss' law) is fixed by a projection, which is defined locally A [s x ] on the "physical" subspace {|s x } with A [s] nc,n ψ = s|n i c , n j ψ , where n i c , n j ψ is some configuration of occupations of fermionic modes c i and ψ j .
Finally, (ii) the second gauge symmetry is controlled by the fermionic number on the link, which is ensured by the product of the nearest-neighbour projectors A [s x ] being non-zero only when N = i n i c,y + n i c,x . The U (1) gauge invariant model in (1 + 1) dimensions is defined by the Hamiltonian, where ψ x are spin-less fermionic operators with staggered mass term µ living on the vertices of the one-dimensional lattice. The bosonic operators E x,x+1 and U x,x+1 , the electric and gauge fields, live on the links of the one-dimensional lattice.
The Hamiltonian is invariant under local U (1) symmetry transformations, and also it is invariant under the discrete parity transformation P and charge conjugation C. The total electric flux, E = x E x,x+1 /L is the order parameter and locates the transition. It is zero in the disordered phase, non-zero in the ordered phase, and changes sign under the C or P symmetry, i.e., P E = C E = −E.
In this framework, in [74] the phase diagram of (1+1)D quantum link version of the Schwinger model is characterised in an external classical background electric field: the quantum phase transition from a charge and parity ordered phase with non-zero electric flux to a disordered one with a net zero electric flux configuration is described by the Ising universality class (see Fig. 5). The thermodynamical properties and phase diagram of a onedimensional U (1) quantum link model are characterised, concluding that the model with half-integer link representation has the same physical properties as the model with integer link representation in a classical background electric field E 0 = 1 2 .

Tensor Networks for Lattice Gauge Theories with continuous groups[72]
The main difference between lattice gauge theories and generic many-body theories is that they require to work on an artificially enlarged Hilbert space, where the action of the group that generates the local invariance can be defined. The physical Hilbert space [96] is then embedded into the tensor product Hilbert space of the constituents by restricting it to those states that fulfill the Gauss law, Figure 6: The Hilbert space H of a quantum many body system (represented here by a 3D box) is exponentially large, since it is the tensor product of the Hilbert spaces of the constituents. Gauge symmetry allows to identify a smaller space, called the physical Hilbert space H P . This is the subspace spanned by those states that fulfill all the local constraints imposed by the gauge symmetry and is represented by a membrane inside H.
that is to those states that are gauge invariant (see Fig. 6 for a graphical description). A generic gauge transformation is built out of local operators A s (g) that represent the local rotation at site s corresponding to a certain element of the group g. The physical Hilbert space (or gauge invariant Hilbert space) H p is defined as the space spanned by all those states that are invariant under all A s (g), where s are the sites of the lattice Λ, L is the number of links, and g is an arbitrary group element.
In [97] the group algebra C(G) is considered as the local Hilbert space, as suggested in the original Hamiltonian description of lattice gauge theories [79,98]. In [97] by exploiting the locality of the operators A s (g) and the fact that they mutually commute, it is shown that the projection onto H p is compatible with a tensor network structure. In particular, the projector is built as hierarchical tensor networks such as the MERA [99] and the Tree Tensor Network [100]. While the MERA is computationally very demanding, a hybrid version of it has been built, that allows to construct the physical Hilbert space by using a MERA and then use a Tree Tensor network on the physical Hilbert space as a variational ansatz. In the same paper, it is also highlighted how the construction of a physical Hilbert space can be understood as a specific case of a duality such as the well known duality between the Z(2) gauge theory and the Ising model [101].
The idea [97] is very flexible and general but strongly relies on using C(G) as the local Hilbert space for every constituent. Since the group algebra contains an orthogonal state for every distinct group element, g, the local Hilbert space becomes infinite dimensional in the case of continuous groups such as e.g. U (1) and SU (N ).
Furthermore, the numerical results with iPEPS in the context of strongly correlated fermions in two dimensions were very promising [102], and thus it was decided to generalise the construction to PEPS tensor networks in [72]. There, it was understood that there is a unifying frame- Figure 7: The projector on the gauge invariant states defined through the contraction of the two tensors C that copy the physical Hilbert space onto the auxiliary Hilbert space and G that selects only configurations fulfilling the gauge invariance condition. The case of a 4 × 4 square lattice with PBC is presented.
work for all the Hamiltonian formulations of lattice gauge theories that can be based on a celebrated theorem in group theory, stating that the group algebra can be decomposed as the sum of all possible irreducible representations C(G) = ⊕ r (r ⊗r), where r is an irreducible representation andr is its conjugate (see Figs. 7 and 8). If the group is compact, the irreducible representations are finite dimensional.
By decomposing C(G) into the direct sum of all the irreducible representations and truncating the sum to only a finite number of them, a formulation of LGT is obtained on finite dimensional Hilbert spaces. For Abelian gauge theories furthermore this procedure [39] leads to the already known gauge magnets or link models [48,49,50].
With this group theoretical picture in mind, it is very easy to directly construct both the projector onto the physical Hilbert space as a tensor network, and tensor network ansatz for states defined on it. The general recipe is given in [72]. Here for concreteness, the construction is shown for a two-dimensional square lattice. The tensor network is composed of two elementary tensors. The first one, C α,j i,β , a four-index tensor that has all elements zero except for those corresponding to α = i = β = j. C is applied to each of the lattice sites and acts as a copy tensor that transfers the physical state of the links (encoded in the leg i) to the auxiliary legs α, β.
The two auxiliary legs are introduced to bring the information to the two sites of the lattices that the link connects. Thus, the copy tensor C allows the decoupling of the gauge constraint at the two sites and to impose the Gauss law individually.
This operation is performed at each site by the second type of tensor, G α1α2 α3α4 , onto the trivial irreducible representation contained in the tensor product Hilbert space The contraction of one C for every link with one G for every site gives rise to the desired projector onto H p with the structure of a PEPS.
Alternatively, the projector onto H p can be incorporated into a variational iPEPS ansatz for gauge invariant states, by promoting each of its tensor elements to a degeneracy tensor along the lines used to build symmetric tensor network states first introduced in [103]. The gauge invariant tensor network can thus be interpreted as an iPEPS with a fixed tensor structure dictated by the gauge symmetry, where each element is again a tensor. These last tensors collect the variational parameters of the ansatz.

Phase diagram and dynamical evolution of Lattice Gauge Theories with tensor networks
Despite their impressive success, the standard LGT numerical calculations based on Monte Carlo sampling are of limited use for scenarios that involve a sign problem, as is the case when including a chemical potential. This constitutes a fundamental limitation for LQCD regarding the exploration of the QCD phase diagram at non-zero baryon density. In contrast, TNS methods do not suffer from the sign problem, which makes them a suitable alternative tool for exploring such problems, although, it is challenging to simulate high-dimensional systems.
In this section, it is shown how tensor network techniques could go beyond Monte Carlo calculations, in the sense, of being able to perform real-time calculations and phase diagrams with finite density of fermions. Examples of these achievements appear in [81,104,105,106,107,108,109,110].

Real-time Dynamics in U(1) Lattice Gauge Theories with Tensor Networks[105]
One of the main applications of tensor network methods is real-time dynamics. Motivated by experimental proposals to realise quantum link model dynamics in optical lattice experiments, Ref. [105] studied the quench dynamics taking place in quantum link models (QLMs) when starting from an initial product state (which is typically one of the simplest experimental protocols). In particular, the model under investigation was the U(1) QLM with S = 1 variables as quantum links, whose dynamics is defined by the Hamiltonian where ψ x defines staggered fermionic fields, and E x,x+1 = S z x,x+1 are quantum link spin variables; while the three Hamiltonian terms describe minimal coupling, mass, and electric field potential energy, respectively.
Several types of time evolutions were investigated. Fig. 9 presents the time evolution corresponding to string breaking dynamics: the initial state, schematically depicted on the top of the main panel, consists of a charge and anticharge separated by a string of electric field (red region), and surrounded by the bare vacuum (light yellow). After quenching the Hamiltonian dynamics (in this specific instance, with m = g = 0), the string between the two dynamical charges breaks (as indicated by a mean value of the electric field around 0 after a time τ t 2), and the charges spread in the vacuum region. For this specific parameter range, an anti-string is created at intermediate time-scales τ t 4. Such string dynamics has also a rather clear signature in the entanglement pattern of the evolving state: in particular, it was shown how the speed of propagation of the particle wave-front extracted from the local value of the electric field was in very good agreement with the one extracted from the bipartite entanglement entropy.
With the same algorithm, it is possible to simulate the time evolution of a rather rich class of initial states up to intermediate times. As another example, Ref. [105] also investigated the scattering taking place between cartoon meson states at strong coupling g 2 1 (at smaller coupling, investigating scattering requires a careful initial state preparation, where a finite-momentum eigenstate is inserted ad hoc onto the MPS describing the dressed vacuum). Rather surprisingly, even these simplified scattering processes were found to generate a single bit of entanglement in a very precise manner.
Closer to an atomic physics implementation with Rydberg atoms is the work [111] where different string dynamics are explored to infer information about the Schwinger model.

4.2.2
Finite-density phase diagram of a (1+1)-d non-Abelian lattice gauge theory with tensor networks [106] By means of the tools introduced in Sec. 4.1.3, in [106,112] the authors have studied the finite-density phase diagram In particular, they introduced a quantum link formulation of an SU (2) gauge invariant model by means of the Hamiltonian where the first term introduces the coupling between gauge fields and matter, as j+1,s + H.c., (6) where c and U are the matter field and parallel transport operators, j ∈ {1, · · · , L − 1} numbers the lattice sites, and s ∈ {↑, ↓}. The second gauge term accounts for the energy of the free field written in terms of the fermion occupation n j,s , where g 1 = g 0 3/8. The last term H break has to be introduced to resolve the undesired accidental local conservation of the number of fermions s=↑,↓ n j,s around every site j, that results in a U (2) theory. This last term breaks this invariance and thus, the final theory is an SU (2) gauge invariant one [106]. Finally, by means of a finite-size scaling analysis of correlation functions, the study of the entanglement entropy and fitting of the central charge of the corresponding conformal field theory, the authors present the rich finite-density phase diagram of the Hamiltonian (5), as reported in Fig. 11. In particular, they identify different phases, some of them appearing only at finite densities and supported also by some perturbative analysis for small couplings. At unit filling the system undergoes a phase transition from a meson superfluid, or meson BCS, state to a charge density wave via spontaneous chiral symmetry breaking. At filling two-thirds, a charge density wave of mesons spreading over neighbouring sites appears, while for all other fillings explored, the chiral symmetry is restored almost everywhere, and the meson superfluid becomes a simple liquid at strong couplings.
Very recently, also a one-dimensional SU (3) gauge theory has been studied with the same approach. In [112], working on and extending the results reviewed in the previous paragraph, the authors present an SU (3) gauge invariant model in the quantum link formulation, and perform an extended numerical analysis on the different phases of the model. For space reasons, the model Hamiltonian is not displayed here and the interested reader is referred to the original publication. However, the main results are: at filling ν = 3/2, the Kogut-Susskind vacuum, a competition between a chiral and a dimer phase separated by a small gapless window has been reported. Elsewhere, only a baryonic liquid is found. The authors also studied the binding energies between excess quarks on top of the vacuum, finding that a single baryon state (three excess quarks) is a strongly bound one, while two baryons (six quarks) weakly repel each other. The authors concluded that the studied theory -differently from three dimensional QCD -disfavours baryon aggregates, such as atomic nuclei.

Density Induced Phase Transitions in the Schwinger Model: A Study with Matrix Product States[107]
The potential of TNS methods to deal with scenarios where standard Monte Carlo techniques are plagued by the sign problem was explicitly demonstrated in [107] by studying the multi-flavour Schwinger model, in a regime where conventional Monte Carlo suffers from the sign problem. The Hamiltonian of the Schwinger model (1) can be modified to include several fermionic flavours, with independent masses and chemical potential, that do not interact directly with each other but only through the gauge field. In the case of two-flavours with equal masses studied in [107] the model has an SU (2) isospin symmetry between the flavours. For vanishing fermion mass and systems of fixed volume, the analytical results [113,114] demonstrate the existence of an infinite number of particle number sectors, characterised by the imbalance between the number of fermions of both flavours. The different phases are separated by first order phase transitions that occur at fixed and equally separated values of the (rescaled) isospin chemical potential, independent of the volume, so that the isospin number of the ground state varies in steps as a function of the chemical potential (see left panel of Fig.  12).
The numerical calculations in [107] used MPS with open boundary conditions, and followed the procedure described in [67], with the exception that the lattices considered had constant volume, Lg = N/ √ x. Thus there was no need for a finite-size extrapolation of the lattice results (N → ∞), but a sufficiently large physical volume Lg had to be considered. The results reproduced with great accuracy the analytical predictions at zero mass, as shown in the left panel of Fig. 12. While the analytical results can only cope with massless fermions, the MPS calculation can be immediately extended to the massive case, for which no exact results exist. For varying fermion masses, the phase structure was observed to vary significantly. The location of the transitions for massive fermions depends on the volume, and the size of the steps is no longer constant. The right panel of Fig. 12 shows the mass vs. isospin phase diagram for volume Lg = 8. The MPS obtained for the ground state allows further investigation of its properties in the different phases. In particular, the spatial structure of the chiral condensate, which was studied in [107] and [115]. While for ∆N = 0, the condensate is homogeneous, for non-vanishing isospin number it presents oscillations, with an amplitude close to the zero density condensate value and a wave-length that, for a given volume, decreases with the isospin number or imbalance, but decreases with Lg.

Efficient Basis Formulation for (1+1)-Dimensional
SU (2) Lattice Gauge Theory: Spectral calculations with matrix product states [108] A non-Abelian gauge symmetry introduces one further step in complexity with respect to the Schwinger model, even in 1+1 dimensions. The simplest case, a continuum SU (2) gauge theory involving two fermion colours, was studied numerically with MPS in [108], using a lattice formulation and numerical analysis in the spirit of [67].
The discrete Hamiltonian in the staggered fermion formulation reads [79] The link operators U n are SU (2) matrices in the fundamental representation, and can be interpreted as rotation matrices. The Gauss law constraint is now non-Abelian, of the non-Abelian charge at site m (if there are external charges, they should be added to Q m ). L τ and R τ , τ ∈ {x, y, z}, generate left and right gauge transformations on the link, and the colour-electric flux energy is The Hilbert space of each link is analogous to that of a quantum rotor, and its basis elements, for the m-th link, can be labeled by the eigenvalues of J m , L z m and R z m , as |j m m m . As in the case of the Schwinger model, it is possible to truncate the gauge degrees of freedom. This can be achieved in a gauge invariant manner [46] and was applied to study the string breaking phenomenon in the discrete theory in [116]. However, in order to attain precise results that permit the extraction of a continuum limit, it is convenient to work in a more efficient basis, in which gauge degrees of freedom are integrated out. A first step to reduce the number of spurious variables is the color neutral basis introduced in [117,118]. In [108], building on that construction, a new formulation of the model on the physical subspace is introduced in which the gauge degrees of freedom are completely integrated out. Nevertheless, it is still possible to truncate the maximum colour-electric flux at a finite value j max in a gauge invariant manner, and analyse the effect of this truncation on the physics of the model. This is relevant, for instance, to understand how to extract continuum quantities from a potential quantum simulation of the truncated theory. To this end, different quantities were computed and extrapolated to the continuum, including the ground state energy density, the entanglement entropy in the ground state, the vector mass gap and its critical exponent for values of the maximum colour-electric flux j max = 1/2, 1, 3/2, 2.
The results demonstrated that, while a small truncation is enough to obtain the correct continuum extrapolated ground state energy density, the situation varies for the mass gap. In particular (see left panel of Fig.13), if the truncation is too drastic, it fails to produce a reliable extrapolation, and only j max > 1 allowed for precise estimations of the vector mass, and the extraction of a critical exponent as the gap closes for massless fermions.
Particularly interesting is the study of the entanglement entropy of the vacuum, which can be easily computed from the MPS ansatz, as already demonstrated in [84] for the Schwinger model. The gauge constraints are not local with respect to a straightforward bipartition of the Hilbert space [119,120], and different contributions to the entropy can be identified, of which only one is distillable, while the others respond only to the gauge invariant structure of the state. In [108], these different contributions were computed and their scaling analysed (see Fig.  14 a). For a massive relativistic QFT, as is the case here for non-vanishing fermion mass, the total entanglement entropy is predicted to diverge as S = (c/6) log 2 (ξ/a) [92], where c is the central charge of the conformal field theory describing the system at the critical point, in the SU (2) case, c = 2. This effect could also be studied from the MPS data ( Fig. 14 d), and it was also found to be sensitive to the truncation, leading to the conclusion that truncations Final value of the vector mass gap (after continuum extrapolation) as a function of the fermion mass using truncations j max = 1 (red triangles), 3/2 (green squares) and 2 (magenta diamonds), with the yellow stars showing results from a strong coupling expansion [118]. The blue circles correspond to j max = 1/2, although the continuum estimation is not reliable in that case. The dotted lines represent the best fit of the form γ(m/g) ν . Right: Central charges extracted from the scaling of the entanglement entropy (see panel d in Fig. 14) for different fermion masses and j max = 1/2 (blue circles), 1 (red triangles), 3/2 (green squares), and 2(magenta diamonds).  (2) LGT. In (a) the different contributions (distillable in blue circles) are shown for the entanglement entropy of a chain of 200 sites, for fermion mass m/g = 0.8 and using j max = 2. Panels (b) and (c) show the extrapolations in bond dimension and system size, and panel (d) shows the continuum limit for the total entropy, exhibiting the divergent log ag term.
of j max ≤ 1 would not recover the continuum theory in the limit of vanishing lattice spacing, as shown by Fig. 13.

4.2.5
Gaussian states for the variational study of (1+1)-dimensional lattice gauge models [109] Gaussian states [121,122,123], whose density matrix can be expressed as the exponential of a quadratic function in the creation and annihilation operators, are widely used to describe fermionic as well as bosonic quantum manybody systems. They fulfill Wick's theorem and thus can be completely described in terms of a covariance matrix, with a dimension that scales only linearly in the system size. This provides a very efficient representation of the quantum many-body state, which can be used as a variational ansatz. But in systems with interacting bosons and fermions, as is the case for lattice gauge theories with gauge and matter degrees of freedom, Gaussian states present a severe limitation, since they cannot describe any correlations between the two types of fields. However, as was recently shown [124], generalised ansätze that combine non-Gaussian unitary transformations with a Gaussian ansatz in the suitable basis, can be successfully used to approximate static and dynamic properties of systems containing fermions and bosons, also in higher dimensions. In [109] this approach was shown to work for (1+1)-dimensional lattice gauge theories. More explicitly, a set of unitary transformations was introduced that completely disentangle the gauge and matter degrees of freedom for any gauge symmetry given by a compact Lie group and a unitary representation. This allows for new ways of studying these lattice gauge theories.
The particular cases of U(1) and SU (2) were explicitly studied in [109]. For U(1), the resulting Hamiltonian is the one proposed by Hamer, Weihong, and Oitmaa in [125], which has been used for numerical calculations [125,67,104], and has been experimentally implemented with trapped ions in a pioneering quantum simulation [47]. The general character of the decoupling transformations thus provides alternative formulations of other lattice gauge theories which can be suitable for experimental implementation with the advantage of being directly defined in the physical space and not requiring the explicit realisation of any gauge degrees of freedom.
With a numerical perspective, [109] addressed the decoupled formulation using a Gaussian variational ansatz, and used it to investigate static and dynamical aspects of string breaking in the Abelian U(1) and non-Abelian SU (2) gauge models. In the U(1) case, the formulation directly allows the study of the real-time string breaking phenomenon in the presence of static external or dynamic charge. In the SU (2) case, only the case of static external charges was studied, using another unitary transformation that decouples them from the dynamical fermions. The Gaussian approach was capable of capturing the essential features of the phenomenon, both the static properties and the out-of-equilibrium dynamics (see Fig. 15). The results showed excellent agreement with previous TNS simulations over a broad range of the parameter space, despite the number of variational parameters in the Gaussian ansatz being much smaller. The approach could be extended and used for further non-equilibrium simulations of other LGTs.  (2) (right column) LGT (from [109]). The upper row of plots shows the static potential between two external charges as a function of the distance, at fixed lattice spacing, for different values of the fermion mass, computed with MPS (crosses on the left, circles on the right) and a Gaussian ansatz (solid lines on the left, triangles and asterisks on the right). The left plot corresponds to two unit charges in the U (1) LGT, and several fermion masses, while the right plot shows the corresponding calculation in the SU (2) case for a pair of external static charges carrying s = 1/2. The lower row shows the real-time evolution of a string created on top of the interacting vacuum. On the left, for the U(1) case, the edges are dynamical charges and can propagate, while on the right, for SU(2), they are static. In both cases, the color indicates the chromoelectric flux on each link as a function of time.
4.2.6 Thermal evolution of the Schwinger model [81,104] TNS ansätze can also describe density operators, in particular thermal equilibrium states, and can therefore be used to study the behaviour of a LGT at finite temperature. This approach was followed in [104,81], which employed a purification ansatz [126] to represent the thermal equilibrium state as a matrix product operator (MPO). At infinite temperature, gβ = 0, the thermal equilibrium state is maximally mixed, and has an exact representation as a simple MPO. By applying imaginary time evolution on this MPO [126,127], a whole range of temperatures can be studied. A relevant observable to analyse in the case of the Schwinger model is again the chiral condensate (see Fig.16). In the massless case, the chiral symmetry is broken (due to an anomaly), and the condensate has a nonzero value in the ground state. The symmetry is smoothly restored at infinite temperature, as demonstrated analytically in [128].
In [104,81] a finite size MPO ansatz was used in the physical subspace, i.e., after integrating out the gauge degrees of freedom. The imaginary time (thermal) evolution was applied in discrete steps, making use of a Suzuki-Trotter approximation, and after each step a variational  Figure 16: Chiral condensate in thermal equilibrium in the continuum as a function of temperature for massless (left, from [104]) and massive fermions (right, from [81]). The solid line on the left shows the analytical prediction for the restoration of the chiral symmetry, while the data points show results obtained with MPO using two different approximations for the evolution operators. For the massive case on the right the horizontal line shows the value at zero temperature (computed numerically with MPS) and the solid blue line the approximation by Hosotani and Rodriguez [129] (which is exact only at very small masses).
optimisation was used to truncate the bond dimension of the ansatz. In the physical subspace, the long-range interactions among charges pose a problem for standard TN approximations of the exponential evolution operators. Two alternatives were considered to apply the discrete steps as MPO. In one of them, the long-range exponential was approximated by a Taylor expansion. In the other, an exact MPO expression of the exponential of the long-range term was used, taking advantage of the fact that it is diagonal in the occupation basis. For the application to be efficient, a truncation was introduced in this MPO, which was equivalent to a cut-off in the electric flux that any link in the lattice can carry. The second approach was found to be more efficient, and a small cut-off was sufficient for convergence over the whole range of parameters explored. With this method, the chiral condensate in the continuum was evaluated from inverse temperature gβ = 0 to gβ ∼ O(10) both for massless and massive fermions. For non-vanishing fermion masses, the condensate diverges in the continuum, and a renormalisation scheme has to be adopted that subtracts the divergence. In [81] this was achieved by subtracting the value of the condensate at zero temperature in the non-interacting case, after the finitesize extrapolation. The continuum limit was performed for each value of the temperature in a manner similar to [67]. The width of the time step in the Trotter approximation introduced an additional source of error, and required an additional extrapolation. However, the form of the step width extrapolation is given by the order of the Suzuki-Trotter approximation, and this step did not affect the final precision, which again turned out to be controlled by the continuum limit.
All in all, the technique allowed for reliable extrapolations in bond dimension, step width, system size and lattice spacing, with a systematic estimation and control of all error sources involved in the calculation. Notably, although the large temperature regime of the lattice model is easier to describe by a MPO, the lattice effects are also more important, which resulted in larger errors after the continuum extrapolation. As the temperature decreases, the errors from the lattice effects become less relevant, but the truncation errors from the MPO approximation accumulate, so that they dominate the low-temperature regime. In conclusion, these results further validate the TNS approach as a tool to investigate the phase diagram of a quantum gauge theory.

4.2.7
Finite temperature and real-time simulation of the Schwinger model [73,86,87] Finite temperature. In [86] different aspects of the finite temperature physics of the Schwinger model were studied. For different temperatures the appropriate gauge invariant Gibbs states were obtained from imaginary time evolution on a purification of the identity operator. Among the different results here the computation of the temperature dependent renormalised chiral condensate is quoted, in agreement with the analytical result for m/g = 0 and the numerical results of [104] for m/g = 0. Furthermore, the study of the temperature-dependence of the energy density in an electric background field allowed for the study of an effective deconfinement transition. For half-integer background fields the expected restoration of the C symmetry at non-zero temperature was also verified.
Real-time simulations. Finally, some of the most relevant results on the real-time simulations of [73,87] are: an intriguing effect in the Schwinger model concerns the non-equilibrium dynamics after a quench that is induced by the application of a uniform electric field onto the ground state at time t = 0. Physically, this process corresponds to the so-called Schwinger pair creation mechanism [89] in which an external electric field separates virtual electron-positron dipoles to become real electrons and positrons. The original derivation [89] involved a classical background field, neglecting any back-reaction of the created particle pairs on the electric field. In [130,131] this back-reaction was taken into account at the semi-classical level. The real-time MPS simulations of [73,87] provide the first full quantum simulation of this non-equilibrium process. In Fig. 17 a sample result is shown, for the case of large initial electric fields. After a brief initial time interval of particle production, a regime can be observed with damped oscillations of the electric field and particle densities going to a constant value. This is in line with the semi-classical results, but as can be seen in figure 17(b) there is a quantitative difference, with a stronger damping, especially for the smaller fields. Finite temperature MPS simulations also allow a comparison with the purported equilibrium thermal values ( figure 17 (a)).

Phase Diagram and Conformal String Excitations of Square Ice using Gauge Invariant Matrix Product
States [110] The examples discussed above widely demonstrate the computational capabilities of tensor network methods in dealing with (1+1)-d lattice gauge theories. Ref. [110] reports instead results on a two-dimensional U (1) gauge theory, the (2+1)-d quantum link model, also known as square ice (for tensor network results on a theory with discrete gauge group, see Ref. [97]).
The main difference between square ice and a conventional U (1) LGT is that the gauge fields now span a two-dimensional Hilbert space, and parallel transporters are replaced by spin operators. The system Hamiltonian reads: where the summation goes over all plaquettes of a square lattice, and the plaquette operator f = σ + µ1 σ + µ2 σ − µ3 σ − µ4 + H.c. flips the spins on the links µ 1 , ..., µ 4 of oriented plaquettes. The first term corresponds to the magnetic field interaction energy, while the second term is a potential energy for flippable plaquettes. There is no direct electric field energy since the spin representation is S = 1/2.
The phase diagram of the model has been determined using a variety of methods, including exact diagonalisation [132] and quantum Monte Carlo [51]. There are two critical points: one is the so-called Rokshar-Kivelson point at λ = 1, where the ground state wave function is factorised into an equal weight superposition of closed loops [52]. This points separates a columnar phase at λ > 1 from a resonating valence-bond solid (RVBS). The latter is separated from a Néel phase by a weak first order transition point at around λ 0.36 [132,51]. All of these phases are confining. The richness of its phase diagram and the possibility of carrying out precise MC simulations make this an ideal model for testing tensor network techniques for (2+1)-d lattice gauge theories.
Ref. [110] presents an analysis based on several observables computed in L x ×L y cylinder geometries to mitigate entanglement growth as a function of the system size. The method of choice was an iTEBD algorithm applied on an MPS ansatz. Beyond simplicity and numerical stability of the algorithm, the main technical advantage of this approach is that re-arranging the MPS in columns allows the integration of the Gauss law in a relatively simple manner.
In the first part, conventional LGT diagnostics, such as the scaling of order parameters for the ordered phase, and the decay of Wilson loops, were analysed. The main conclusion is that TN methods can reach system sizes well beyond ED with the necessary accuracy for determining order parameters and correlation functions. However, the system sizes achieved (up to 600 spins) were smaller when compared to the ones accessible with QMC: this prevented, for instance, a systematic study of Wilson loops, that were found to be particularly sensitive to finite volumes and open boundary conditions. The second part of the work instead focused on entanglement properties of string states, which are generated by introducing two static charges in the system at given distance . Some results on the entropy difference between the string and ground states are depicted in Fig. 18a-b: in the region of space between the two charges (n ∈ [11,20] in the panels), this entropy difference is compatible with a conformal field theory scaling, with a central charge that is compatible with 1 in a large parameter regime within the RVBS phase (Fig. 18c), where finite volume effects are moderate. These results represent an entanglement proof of the conformal behaviour of string excitations in a confining theory, which is a well established fact in non-Abelian (3+1)-d cases as determined from the string energetics [133].

Quantum computation and digital quantum simulation
There are two avenues towards quantum simulations -analog and digital. In analog simulations, the degrees of freedom of the original system and the dynamical evolutions, are mapped to the simulating system. In digital simulations, the simulating system is evolving forward in time stroboscopically, by applying a sequence of short quantum operations. In this section, the digital quantum simulation approach to high-energy models is reviewed, while the analogue quantum simulation is described in the following one.

Quantum and Hybrid Algorithms for Quantum Field Theories
Quantum information, in general, and quantum computation, in particular, have brought new tools and per- spectives for the calculation and computation of strongly correlated quantum systems. Understanding a dynamical process as a quantum circuit or the action of a measurement as a projection in a Hilbert space are just two instances of this quantum framework. In this section, two relevant articles [134,135] are described where these new approaches are used.

Quantum Algorithms for Quantum Field Theories[134]
Quantum computers can efficiently calculate scattering probabilities in φ 4 theory to arbitrary precision at both weak and strong coupling with real-time dynamics, contrary to what is achieved in LGT where the scattering data can also be computed in Euclidean simulations [136,137,138]. In [134], Jordan et al. developed a (constructive) quantum algorithm that could compute relativistic scattering probabilities in a massive quantum field theory with quartic self-interactions (φ 4 theory) in space-time of four and fewer dimensions and solve the equations of QFT efficiently that can be compared with the data from particle accelerators. The proposed algorithm is polynomial in the number of particles, their energy, and the desired precision. In the limit of the so-called strong coupling of QFT, the algorithm actually provides an exponential acceleration with respect to the best known classical algorithms. This work is based on three important technical achievements. First, continuous fields can be accurately represented by a finite number of qubits whose coordinates form a lattice. This achievement is highly non-trivial, because QFTs are contaminated by infinite values of various quantities that must be cured by using renormalisation and regularisation methods, both of which feature naturally in the discussed algorithm. Second, one bottleneck for an efficient implementation of the simulation, the preparation of the initial state, is achieved by a preparation of particles in the form of wave packets. Third, the time evolution is split into the action of local quantum gates. This procedure works well for local theories (field theories discretised on a finite space-time mesh) whose accuracy and convergence must be controlled. Hence, quantum computation for continuous fields can be achieved in a controlled way and with an exponential quantum speedup.
More concretely, the scalar field Hamiltonian in D − 1 spatial dimensions reads H = H π + H φ with The conjugate variables φ(x) and π(x) obey the canonical commutation relations [φ (x) , φ (y)] = [π (x) , π (y)] = 0, [φ (x) , π (y)] = i a D−1 δ (x, y). If the coefficient λ 0 vanishes, then the Hamiltonian is quadratic in the variables φ and π. In that case, the theory is Gaussian, describes a massive non-interacting particle and it can be solved exactly. The complete Hamiltonian is the sum of two terms, one is diagonal in the π basis, while the other is diagonal in the φ basis. Choosing a small time step , then exp It is easy to simulate time evolution governed by the diagonal Hamiltonian H π or H φ , evolving the system using the field Fourier transform to alternate back and forth between the π and φ bases, and applying a diagonal evolution operator in each small time step.

Simulations of Subatomic Many-Body Physics on a Quantum Frequency Processor[135]
The emerging paradigm for solving optimisation problems using near-term quantum technology is a kind of hybrid quantum-classical algorithm. In this scheme, a quantum processor prepares an n-qubit state, then all the qubits are measured and the measurement outcomes are processed using a classical optimiser; this classical optimiser instructs the quantum processor to alter slightly how the n-qubit state is prepared. This cycle is repeated many times until it converges to a quantum state from which the approximate solution can be extracted. When applied to classical combinatorial optimisation problems, this procedure goes by the name Quantum Approximate Optimisation Algorithm (QAOA) [139]. But it can also be applied to quantum problems, like finding low-energy states of many-particle quantum systems (large molecules, for example). When applied to quantum problems this hybrid quantum-classical procedure goes by the name Variational Quantum Eigensolver (VQE).
Hence, VQE algorithms provide a scalable path to solve grand challenge problems in subatomic physics on quantum devices in the near future. One way to implement VQE optically is using a quantum frequency processor (QFP). A variety of basic quantum functionalities have recently been demonstrated experimentally in this approach. A QFP is a photonic device that processes quantum information encoded in a comb of equi-spaced narrow band frequency bins. Mathematically, the QFP is described by a unitary mode transformation matrix V that connects input and output modes.
In [135] it was demonstrated how augmenting classical calculations with their quantum counterparts offers a roadmap for quantum-enabled subatomic physics simulations. More concretely, a subatomic system can be simulated as a collection of nucleons with effective field theory (EFT) parameters input from experimental data or ab-initio calculations. Using a photonic QFP, the ground state energies of several light nuclei using experimentally determined EFT parameters were computed in [135].
The VQE algorithm calculates the binding energies of the atomic nuclei 3 H, 3 He, and 4 He. Further, for the first time, a VQE was employed to determine the effective interaction potential between composite particles directly from an underlying lattice quantum gauge field theory, the Schwinger model. This serves as an important demonstration of how EFTs themselves can be both implemented and determined from first principles by means of quantum simulations.

Digital quantum simulation with trapped ions
Due to the high degree of quantum control of trapped ions platforms, they can be seen as prototypes of universal quantum simulators. In the following sections, two applications are described that realise experimentally the idea of a quantum simulator for high-energy processes.

5.2.1
Real-time dynamics of lattice gauge theories with a few-qubit quantum computer [47,140] The article [47] reports on the first digital quantum simulation of a gauge theory from high-energy physics and entails a theoretical proposal along with its realisation on a trapped ion quantum computer [141,142]. This simulation addresses quantum electrodynamics in one spatial and one temporal dimension, the so-called Schwinger model [143,144].
In [47,140], the coherent real-time dynamics of spontaneous particle-anti-particle pair creation has been studied. Such dynamical processes cannot be addressed with conventional Markov Chain Monte Carlo methods due to the sign problem [145] preventing the simulation of time evolutions. The experiment performed in [47,140] realises a Trotter time evolution [146] based on the Kogut-Susskind Hamiltonian formulation of the Schwinger model [147]. The simulation protocol used in this demonstration is customtailored to the experimental platform and based on eliminating the gauge degrees of freedom. The Schwinger model entails dynamical matter and gauge fields (electromagnetic fields). The gauge degrees of freedom are analytically integrated out [148], which leads to a pure matter model that can be cast in the form of an exotic spin model that features two-body terms and long-range interactions. The gauge bosons do not appear explicitly in the description but are included implicitly in the form of long-range interactions. The resulting encoded model is gauge invariant at all energy scales and allows one to simulate the full infinite-dimensional Hamiltonian. This approach is well matched to simulators based on trapped ions [141,142], which naturally feature a long-range interaction and hence allow for a very efficient implementation of the encoded Schwinger model. The Trotter protocol that has been devised in [47,140] can be realised in a scalable and resourceoptimised fashion. The number of Trotter steps is ideal for the required ion-ion coupling matrix and scales only linearly in the number of lattice sites N . Moreover, the protocol is designed such that Trotter errors do not lead to gauge variant contributions.
The experiment [47,140] has been carried out for N = 4 lattice sites (i.e., using four qubits) and a gate sequence comprising more than 200 gate operations. The used resources are high-fidelity local gate operations and the socalled Mølmer-Sorensen gates [149] with all-to-all connectivity between the individual ions. The qubit states are encoded in electronic sub-levels of the ions. In the experiment, a quench has been performed in which the bare vacuum (i.e., the ground state for infinite fermion mass) has been prepared, followed by a Trotterised time evolution under the encoded Schwinger Hamiltonian, which leads to the generation of particle-anti-particle pairs. This type of experiment can also be performed starting from the dressed vacuum (i.e., an eigenstate of the full Hamiltonian for finite fermion mass), which is a highly entangled state that can be prepared on a trapped ion quantum simulator using the method demonstrated in [150] (see Sec. 5.2.2 below). In this case, a quench to generate pair creation events would involve the time evolution under the Schwinger Hamiltonian including background electric fields. Including electric background fields to the encoded Schwinger Hamiltonian leads to additional local terms and therefore requires only minor modifications in the quantum simulation protocol [140]. Using trapped ions, highfidelity measurements can be made using fluoresce detection [141,142]. In [47,140] local measurements in the zbasis have been used to study the site-resolved particle number density and electric field distribution in real-time After a fast transient pair creation regime, the increased particle density favours particle-anti-particle recombination inducing a decrease of ν(t). This non-equilibrium interplay of regimes with either dominating production or recombination continues over time and leads to an oscillatory behaviour of ν(t) with a slowly decaying envelope. (c) The entanglement entropy S(t) quantifies the entanglement between the left and the right half of the system, generated by the creation of particle-anti-particle pairs that are distributed across the two halves. An increasing particle mass m suppresses the generation of entanglement. From [140] as a function of the fermion mass. This type of analysis can be directly scaled up to larger system sizes. The experiment [47,140] probed also the entanglement generated during pair creation, which can be done for small system sizes and involved the measurement of the density matrix of the spin system. As shown in [47,140], the entanglement of the encoded model corresponds to the entanglement in the original model involving both gauge fields and fermions.

Self-Verifying Variational Quantum Simulation of the Lattice Schwinger Model[150]
In this article, a quantum co-processor successfully simulated particle physics phenomena on 20 qubits for the first time. The experiment uses new methods with a programmable ion trap quantum computer with 20 quantum bits as a quantum co-processor, in which quantum mechanical calculations that reach the limits of classical computers are outsourced. For this, a sophisticated optimisation algorithm has been developed that, after about 100,000 uses of the quantum co-processor by the classi-cal computer, leads to the result. In this way, the programmable variational quantum simulator has simulated the spontaneous creation and destruction of pairs of elementary particles from a vacuum state on 20 quantum bits.
An analog quantum processor prepares trial states, quantum states that are used to evaluate physical quantities. The classical computer analyses the results of these evaluations, with the aim of optimising certain adjustable (variational) parameters on which the trial states depend. This computer then suggests improved parameters to its quantum co-worker in a feedback loop. In the study, the quantum device contains a line of atomic ions that each represent a qubit. This set-up is used to carry out quantum simulations of the ground state of electrons coupled to light, a system that is described by the theory of quantum electrodynamics in one spatial dimension.

Digital quantum simulation with superconducting circuits
This section reviews the possibility to perform digital quantum simulation of lattice gauge theories with superconducting circuits [44].

Non-Abelian SU (2) Lattice Gauge Theories in
Superconducting Circuits [44] Superconducting circuits have proven to be reliable devices that can host quantum information and simulation processes. The possibility to perform quantum gates with high fidelities, together with high coherence times, makes them ideal devices for the realisation of digital quantum simulations. In [44], a digital quantum simulation of a non-Abelian dynamical SU (2) gauge theory is proposed in a superconducting device. The proposal starts from a minimal setup, based on a triangular lattice, that can encode pure-gauge dynamics. The degrees of freedom of a single triangular plaquette of this lattice are encoded into qubits. Two implementations of this quantum simulator are described, using two different superconducting circuit architectures. A setup in which six tunable-coupling transmon qubits are coupled to a single microwave resonator is considered, and a device where six capacitively coupled Xmon qubits stand on a triangular geometry, coupled to a central auxiliary one. The experimental requirements necessary to perform the simulation on one plaquette are characterised and arguments for scaling to large lattices are also given in [44].
A minimal implementation of a pure SU (2) invariant model in a triangular lattice is considered by using triangular plaquettes. In this case, the pure-gauge Hamiltonian on a single plaquette reads (11) This interaction corresponds to the magnetic term of a gauge invariant dynamics, which acts on closed loops.
Due to gauge invariance, the local Hilbert space of a link is four dimensional, and it can be faithfully spanned by two qubits, called "position" σ a pos and "spin" qubit σ a m . In this subspace, it is useful to define the operators Γ 0 = σ x pos σ 0 m , Γ a = σ y pos σ a m , such that the total Hamiltonian is written as In order to simulate the interaction of Eq. (12), one can decompose its dynamics in terms of many-body monomials, and implement them sequentially with a digitised approximation. In a digital approach, one decomposes the dynamics of a Hamiltonian H = (here and in the following = 1), for a total of m × N gates, with an approximation error that goes to zero as the number of repetitions N grows. In a practical experiment, each quantum gate e −ih k t will be affected by a given error k . By piling up sequences of such gates, for small gate errors k 1, the total protocol will be affected by a global error, which is approximately the sum ≈ k k .
To simulate the pure-gauge interaction in a single triangular plaquette, first, a setup with six tunable-coupling transmon qubits coupled to a single microwave resonator is considered. Each tunable-coupling qubit is built using three superconducting islands, connected by two SQUID loops. Acting on these loops with magnetic fluxes, one can modify the coupling of the qubits with the resonator, without changing their transition frequencies. By threading with magnetic fluxes at high frequencies, one can drive simultaneous red and blue detuned sidebands, and perform collective gates. Each many-body operator can be realised as a sequence of collective and single-qubit gates.
A second architecture is considered where six Xmon qubits in a triangular geometry are capacitively coupled with an additional central ancillary qubit. In this case, the collective interactions can be decomposed and performed with pairwise C-phase gates, using the central ancillary qubit to mediate non-nearest-neighbour interactions. In this way, the quantum simulation of one digital step of the Hamiltonian in Eq. (12) will amount to realise 168 C-phase gates and a number of single-qubit rotations which is upper bounded by 520.

Digital quantum simulation with ultra-cold atoms
Digital quantum simulators show the possibility to achieve universal quantum computation. Among the most promising platforms are the ones built with ultra-cold atoms. In this section, several instances are shown using optical lattices and Rydberg platforms where even a completely gauge invariant simulation could be achieved [151,39,42,152,153].

A Rydberg Quantum Simulator[151]
Figure 21: Setup of the system: a) Two internal states |A i and |B i give rise to an effective spin degree of freedom. These states are coupled to a Rydberg state |R i in two-photon resonance, establishing an electromagnetically induced transparency (EIT) condition. On the other hand, the control atom has two internal states |0 c and |1 c . The state |1 c can be coherently excited to a Rydberg state |r c with Rabi frequency Ω r , and can be optically pumped into the state |0 c for initialising the control qubit. b) For the toric code, the system atoms are located on the links of a two-dimensional square lattice, with the control qubits in the centre of each plaquette for the interaction A p and on the sites of the lattice for the interaction B s . Setup required for the implementation of the color code (c), and the U (1) lattice gauge theory (d). From [151].
A universal quantum simulator is a controlled quantum device that faithfully reproduces the dynamics of any other many-particle quantum system with short-range interactions. This dynamics can refer to both coherent Hamiltonian and dissipative open-system time evolution. Cold atoms in optical lattices, which are formed by counterpropagating laser beams, represent a many-particle quantum system, where the atomic interactions and dynamics of the particles can be controlled at a microscopic level by external fields. This high level of control and flexibility offers the possibility to use these systems as quantum simulators, i.e., as devices which can mimic the behaviour of other complex many body quantum systems and allow the study of their properties, dynamics and phases.
Stored cold atoms in deep lattices, in which atoms do not hop between the lattice sites, can be used to encode quantum bits in different electronic states of the atoms. Interestingly, although the atoms sit at different sites and do not collide, it is possible to induce very strong interactions between atoms separated by distances of several micrometers. This can be achieved by exciting them to electronically high-lying Rydberg states. These Rydberg interactions offer the possibility to realise fast quantum gates between remote atoms. Motivated by and building on these achievements, a digital Rydberg simulator architecture based on sequences of fast and efficient quantum gates between Rydberg atoms is developed in [151]. This "digital" simulator offers promising perspectives for the simulation of complex spin models, which are of great interest both in quantum information science, condensed matter, and high-energy physics.
The proposed simulation architecture allows one to realise a coherent Hamiltonian as well as dissipative opensystem time evolution of spin models involving n-body interactions, such as, e.g., the Kitaev toric code, colour code and lattice gauge theories with spin-liquid phases. The simulator relies on a combination of multi-atom Rydberg gates and optical pumping to implement coherent operations and dissipative processes. Highly excited Rydberg atoms interact very strongly, and it is possible to switch these interactions on and off in a controlled way by applying laser pulses. By choosing on which atoms to shine light, the properties of the quantum simulator can be precisely tuned.
As a key ingredient of the setup, extra auxiliary qubit atoms are introduced in the lattice, which play a twofold role. First, they control and mediate effective n-body spin interactions among a subset of n system spins residing in their neighbourhood of the lattice. This is achieved efficiently, making use of single-site addressability and a parallelised multi-qubit gate, which is based on a combination of strong and long-range Rydberg interactions and electromagnetically induced transparency (EIT). Second, the auxiliary atoms can be optically pumped, thereby providing a dissipative element, which in combination with Rydberg interactions results in effective collective dissipative dynamics of a set of spins located in the vicinity of the auxiliary particle, which itself eventually factors out from the system spin dynamics. The resulting coherent and dissipative dynamics on the lattice can be represented by, and thus simulates a master equation, where the Hamiltonian is the sum of n-body interaction terms, involving a quasi-local collection of spins in the lattice. The Liouvillian term in the Lindblad form governs the dissipative time evolution, where the many-particle quantum jump operators involve products of spin operators in a given neighbourhood.

Optical Abelian Lattice Gauge Theories[39]
In [39], it is described how to perform a digital quantum simulation of the gauge-magnet/quantum link version of a pure U(1) lattice gauge theory with ultra-cold atoms, for a recent proposal of an analogue Rydberg simulator for the same theory see [154]. Its phase diagram has been recently characterised by numerical investigations [51]. The experiment aims at mapping the phase diagram of the spin 1/2 U (1) quantum link model by measuring the string tension of the electric flux tube between two static charges and its dependence on the distance. In the confined phase, the string tension is finite, and thus the energy of the system increases linearly with the inter-charge separation. Charges are thus bound together. In the deconfined phase the string tension vanishes and thus the charges can be arbitrarily far away with only a finite energy cost.
In the proposed quantum simulation the gauge bosons are encoded in the hyper-fine levels of Rydberg atoms. The atoms are in a Mott-insulating phase with one atom per site. Extra atoms are needed in order to collectively and coherently address several atoms at the same time. The simulation requires imposing the Gauss law and engineer the dynamics. The latter is obtained digitally decomposing unitary time evolution in elementary Trotter steps that can be performed by Rydberg gate operations. The former can be imposed by dissipation or by engineering digitally an energy penalty for the forbidden configurations. This is achieved by using the Rydberg blockade as first proposed in [155]. The key ingredient is the mesoscopic Rydberg gate in which one control atom is excited and de-excited from its Rydberg state and as a result of the blockade this affects several atoms inside its blockade radius. The setup thus requires two set of atoms, atoms encoding the gauge boson degrees of freedom (one per link of the lattice) that are called ensemble atoms. These atoms are controlled by addressing another set of atoms, the control atoms. In this setup the control atoms are used in order to imprint the desired dynamics on the ensemble atoms.
In order to simulate the U (1) quantum link model, one control atoms located at the center of every plaquette and one control atom located at every site are used. The ensemble atoms are located at the center of the links of the lattice. The lattice spacing should also be engineered in such a way that only four atoms encoding the gauge boson degrees of freedom should be contained inside the blockade radius of the control atoms. Individually addressing and manipulating the control atoms via, e.g., a quantum-gas microscope is also needed.
With this setup, an arbitrary Hamiltonian can be implemented on the atoms encoding the gauge boson degrees of freedom digitally, by decomposing it into a sequence of elementary operations, involving single-site rotations combined with the use of the mesoscopic Rydberg gate. As a result of the lattice geometry, the gate involves one control atom (either at one site or in the center of one plaquette) and the four ensemble atoms surrounding it. This architecture is indeed sufficient to perform a universal quantum simulation of Abelian lattice gauge theories [151].
The simulation requires two stages. During the first stage one starts from some trivial state and prepares the state to be studied such as, e.g., the ground state of the quantum link Hamiltonian. In a second stage, the mesoscopic Rydberg gates are reversed and the state of the system is transferred to the state of the control atoms, that if appropriately read out (through, e.g., a quantumgas microscope), allow the measurement of the physical state of the system and its properties, such as, e.g., the string tension between two static charges.
The simulations are digital, in the sense that they require applying a discrete sequence of pulses to the atoms, whose nature and duration can be found by using optimal control techniques.

Simulations of non-Abelian gauge theories with optical lattices[42]
An important and necessary step towards the quantum simulation of QCD is the simulation of simpler non-Abelian gauge theories in two dimensions to study the interplay of electric and magnetic interactions with non-Abelian local symmetry. The minimal relevant example is given by SU (2) gauge magnets or quantum link models [49,50] with static charges considered in [42]. There it is shown how to characterise confinement in the model and determine its phase diagram by simulating it digitally with Rydberg atoms.
For SU (2), the quantum link is written as the direct sum of two spins 1 2 sitting at each end of the link, see Fig. 22 a). As in [79], physical states, i.e., configurations allowed by gauge invariance, are determined through the (non-Abelian version of the) Gauss' law, and the dynamics comes from competition of electric (on each link) and magnetic (plaquette) interactions. In SU (2) gauge magnets, the charges occupying the sites of the lattice are also represented as spins and the Gauss law demams that the total spin at each site, i.e., spins 1 2 at the link ends coupled to the static charge residing at the site, is zero. Thus, for spin 1 2 charges, physical states are singlet coverings. The electric term weights them depending on the position of the singlets while the plaquette interchanges singlet coverings (or annihilates them) as shown in Fig. 22 b).
The main features that SU (2) gauge magnets share with QCD (and other non-Abelian gauge theories) are: the nature of confinement phases at weak (plaquettes dominate) and at strong coupling (electric terms dominate) and long-range entanglement between charges. To satisfy the Gauss law, the charges must form singlets with the nearby link spins, thus many singlets must be rearranged, and the allowed singlet coverings are different with respect to the ones of the vacuum, at least along a string between the charges. Such rearrangement generates long-range entanglement and costs an energy that increases linearly with the charge separation, i.e., linear confinement, see in Fig.  22 c). To target such phenomena in a quantum simulator, it is enough to consider static spin 1 2 charges [42]. Both spin 1 2 or qubits on the links and on the sites are repre-sented by ground and Rydberg states of atoms. The non-Abelian Gauss' law is converted into an energy penalty and added to the Hamiltonian. The dynamics of the generalised Hamiltonian is decomposed in a sequence of simultaneous Rabi transfers controlled by ancillary qubits and realised by Rydberg gates [155] see Fig. 22 d), in a similar fashion as done for Abelian gauge theories [151,39]. In such a simulator, the ground state is prepared with a pair of opposite static charges at distance L adiabatically (or super-adiabatically). By measuring the final state of the control qubits the energy of such a ground state can be computed with respect to the vacuum as a function of L, E(L), and thus determines the string tension σ = E(L)/L. If σ is finite for large L, there is a linearly confined phase. The proposed Rydberg simulator can probe confinement at any coupling. By inspecting quantum correlations in the prepared ground state, it is also possible to experimentally access the long-range entanglement due to confinement in non-Abelian gauge theories.

Digital quantum simulation of Z(2)
lattice gauge theories with dynamical fermionic matter [152,153] In a recent work [152,153], a digital scheme was introduced and its implementation with cold atoms was studied, for Z(2) and Z(3) lattice gauge theories. The scheme includes, in addition to the gauge and matter degrees of freedom, auxiliary particles that mediate the interactions and give rise to the desired gauge theory dynamics, by constructing stroboscopically the evolution from small time steps. The individual time steps respect local gauge invariance, so errors due to the digitisation will not break local gauge symmetry. Moreover, it is shown that the required three-and four-body interactions, can be obtained by a sequence of two-body interactions between the physical degrees of freedom and the ancillary particles. The construction is general in form, and valid for any gauge group. Its generality and simplicity follows from the use of a mathematical quantum mechanical object called stator [156,157].
6 Analog Quantum simulations 6.1 Analog quantum simulation of classical gauge potential The complete challenge of quantum simulating a lattice gauge theory has many interesting side products such as the study of classical gauge potential and the related topological insulators. In these cases, the gauge potential appears just as a classical configuration of the vector potential that is described by the minimal Wilson line in the lattice or the Peierls substitution of the hopping term. In the following, several theoretical proposals and an experimental realisation are reviewed [158,159,160,161]. Magnetic interactions exchange parallel singlets on plaquettes and annihilate the other configurations. c) Linear confinement induced by a pair of opposite charges at strong and weak couplings, where the electric and magnetic terms dominate, respectively. d) Implementation scheme without charges: the Gauss law and the plaquette interactions are decomposed in elementary Cnot gates that involve all physical qubits/atoms (in red) within the yellow and blue blockade areas, respectively, of the auxiliary Rydberg atoms (in blue). For the full scheme see [42].

Wilson Fermions and Axion Electrodynamics in Optical Lattices[158]
Ultra-cold Fermi gases in optical superlattices can be used as quantum simulators of relativistic lattice fermions in 3 + 1 dimensions. By exploiting laser-assisted tunnelling, an analogue of the so-called naive Dirac fermions is characterised in [158]. An implementation of Wilson fermions is shown, and it is discussed how their mass can be inverted by tuning the laser intensities. In this regime of the quantum simulator, Maxwell electrodynamics is replaced by axion electrodynamics: a three dimensional topological insulator.
In particular, a concrete proposal for the realisation of laser-assisted tunnelling in a spin-independent optical lattice trapping a multi-spin atomic gas is presented. The setup consists of a spin-independent optical lattice that traps a collection of hyperfine states of the same alkaline atom, to which the different degrees of freedom of the field theory to be simulated are then mapped. Remarkably enough, it is possible to tailor a wide range of spin-flipping hopping operators, which opens an interesting route to push the experiments beyond the standard superfluid-Mott insulator transition. The presented scheme combines bi-chromatic lattices and Raman transfers, to adiabatically eliminate the states trapped in the middle of each lattice link. These states act as simple spectators that assist the tunnelling of atoms between the main minima of the optical lattice. This mechanism is clearly supported by numerical simulations of the time evolution of the atomic population between the different opticallattice sites.
Such a device could have important applications in the quantum simulation of non-interacting lattice field theories, which are characterised, in their discrete version, by on-site and nearest-neighbour hopping Hamilto- nians. Once the fields of the theory to be simulated are mapped into the atomic hyperfine states, the desired operators correspond to population transfers between such levels. The former can be realised by standard microwaves, whereas the latter might be tailored with the laser-assisted schemes. Combining this trapping scheme with Fermi gases, this platform would open a new route towards the simulation of high-energy physics and topological insulators.

Non-Abelian gauge fields and topological insulators in shaken optical lattices[159]
A preliminary step to quantum simulating full-fledged non-Abelian gauge theory is to consider classical non-Abelian gauge fields. It is thus crucial to devise efficient experimental strategies to achieve classical synthetic non-Abelian gauge fields for ultra-cold atoms in the bulk and in optical lattices [159]. The latter situation is especially interesting as it allows for an anomalous quantum Hall effect [162] and, in combination with strong interactions, for fractional quantum Hall states with non-Abelian anyonic excitations [163]. On the lattice, synthetic non-Abelian gauge fields (in LGT language, the Wilson operators on the links) correspond to tunnelling matrices that determine the superposition each atomic species is sent to when it tunnels to a neighbouring site. In [159] it is shown for the first time how to achieve such matrices for the SU (2) gauge group from lattice shaking (for general theory of lattice shaking and Floquet driving see [164]). For simplicity, consider a spin-dependent square lattice described by Fig. 25, which in combination with a uniform magnetic field produces a sublattice-dependent energy splitting between the spin up and down states, e.g., m F = ±1 hyperfine states of 87 Rb. With a constant microwave beam Ω coupling the two spin states, the eigenstates of the on-site Hamiltonian in the two sublattices thus differ by an SU (2) rotation. Thus, an atom tunnelling between the sublattices would experience an SU (2) gauge field that is trivial: the product of the tunnelling matrices around the plaquette L is the identity, i.e., its trace -the Wilson loop-is 2. Furthermore, such tunnelling is highly suppressed as out of resonance due to the energy offset between the sublattices (if it is sufficiently large). However, the energy conservation can be restored, and thus the tunnelling, by shaking the lattice at a frequency commensurable with the energy off-sets. Such analysis can be made precise by time averaging the total Hamiltonian in the rotating frame of the driving plus the on-site Hamiltonian [159]. The main result is that for feasible experimental parameters, generic non-trivial classical SU (2) gauge field configurations can be engineered, i.e., characterised by |Tr L| < 2, as shown in Fig. 25. √ ±∆E + ∆E becomes sub-lattice dependent. In combination with a coupling Ω of both spin states realised microwave (or magnetic) fields the sublattice splitting gives local spin eigenbasis that differ in the two sub-lattices by SU (2) rotation. With lattice shaking, the SU (2) gauge transformation is converted by timeaveraging into a non-trivial synthetic non-Abelian gauge field. c) Accessible Wilson loops |Tr L| for convenient experimental parameters (K y /K x is the relative shaking amplitude in the y and x directions and ω is the frequency). Deviations from 2 imply non-Abelian physics: outside the white (black) regions, |Tr L| < 1.9 ( < 1.99). Figure 26: A convenient approach for analog quantum simulation with ultra-cold neutral atoms relies in the concept of "synthetic dimensions": a coherent coupling Ω between internal atomic states |i (induced e.g. with laser fields) mimics an effective hopping t between sites with position m i of a fictitious synthetic dimensionm.

Observation of chiral edge states with neutral fermions in synthetic Hall ribbons[160]
A powerful resource for the implementation of analog quantum simulations with ultra-cold-atomic samples is based on the manipulation of the internal atomic degrees of freedom. In this context, atoms with two valence electrons (such as alkaline-earth elements or lanthanide ytterbium) represent a convenient choice, as they provide access to several stable internal states (either nuclear or electronic), that can be initialised, manipulated with long-coherence times and read-out optically with high-fidelity. This platform is particularly suitable for implementing the concept of "synthetic dimensions", in which a manifold of internal states is mapped onto effective positions along a fictitious discretised extra-dimension, and the coherent optical coupling between the different states can be described as an effective hopping between synthetic sites (see Fig. 26 for a sketch of the general idea).
This approach, initially proposed in Ref. [165], provided a very convenient method for the realisation of tunable background gauge potentials [166]. This is explained in Fig. 27a, showing a hybrid lattice structure combining one real direction (discretised by a real optical lattice in sites with position j) with a synthetic direction composed by internal atomic states, depicted orthogonally to the real one. The hopping matrix element along the synthetic lattice is a complex quantity Ωe iφj , with an argument depending on the phase of the electric field inducing the transition between the internal states generating the synthetic dimension. As a consequence, hopping around a unit cell of the real and synthetic lattice can be described in terms of an effective geometric Aharonov-Bohm phase φ associated with the effect of a background synthetic magnetic field on effectively charged particles.
This idea was experimentally realised in Ref. [160] (and in a related experiment [167]), where different nuclearspin projection states of fermionic 173 Yb atoms were coupled coherently with a two-photon Raman transition, realising the scheme of Fig. 27a. This system is particularly suited to study edge physics, as the synthetic dimension is made up by a finite number of states/sites, Figure 27: a) Sketch of the experimental configuration employed in Ref. [160] for the generation of classical gauge potentials. b) Measured lattice momentum distribution along the different synthetic legs of the ladder, showing the emergence of steady-state chiral edge currents J. c) Experimental reconstruction of the average trajectory on the synthetic strip, evidencing a non-equilibrium dynamics induced after a quench in the tunnelling. Adapted from Ref. [160].
resulting in a ladder geometry (with a tunable number of legs). The emergence of steady-state chiral currents at the edges of the ladder was detected with a statedependent imaging technique (corresponding to a singleleg detection in momentum-space), evidencing a counterpropagating atomic motion on the two outer legs. This behaviour is shown in the graphs of Fig. 27b for a threeleg ladder made by three nuclear-spin projection states m = −5/2, −1/2, +3/2. The asymmetry of the lattice momentum distribution n(k) (along the real direction) provides a direct measurement of the steady-state edge currents J induced in the system after an adiabatic loading. Non-equilibrium dynamics was also studied after imposing a quench on the system: after preparing the fermionic particles on a single external leg, tunnelling along the rungs was suddenly activated. The ensuing dynamics was studied by reconstructing the average trajectory of the particles on the synthetic strip: the result is shown in Fig.   27c, evidencing a "skipping-orbit" motion, in which the synthetic gauge field bends the centre-of-mass trajectory in a cyclotron-like fashion, with repeated bouncing at the edge of the strip producing a net motion along the edge (akin to the steady-state chiral currents discussed above).
Synthetic dimensions are a quite general concept, that can be adapted to different implementations. In a followup of that experimental work [168], with a different implementation relying on the manipulation of the electronic state of 173 Yb (rather than the nuclear-spin states) with a single-photon optical clock transition, the strength and direction of the chiral currents was measured as a function of the synthetic magnetic flux φ, all the way from zero to above half of a quantum of flux per unit cell (a quantum of flux corresponding to φ = 2π). The measurements are summarised in Fig. 28, where the chiral current J is plotted vs. φ. It is apparent that the chiral current vanishes and then changes direction crossing the φ = π point. This can be explained by recalling the two underlying symmetries of the system: the time-reversal symmetry (broken by the magnetic field) imposing J(φ) = −J(−φ) and the periodicity condition coming from the underlying lattice structure J(φ) = J(φ + 2π). A similar conclusion could be reached by considering the behaviour of an extended two-dimensional system, which realises the Harper-Hofstadter model describing charged particles in a two-dimensional lattice with a transverse magnetic field. The topological invariant, i.e., the Chern number, resulting from that model (shown in the inset of Fig. 28 as a colour-coded shading of the Hofstadter butterfly spectrum) changes from positive to negative values when the φ = π point is crossed: this sign change is connected, by the bulk-boundary correspondence principle, to the reversal of chiral currents measured in the experimentallyrealised ladder geometry.
The concept of synthetic dimensions is promising for the realisation of quantum simulators with advanced control on topological properties and site connectivities. For instance, controlling the parameters of the atom-laser interaction gives the possibility of engineering systems with periodic boundary conditions (i.e., a compactified extradimension) [169] or different kinds of topological ladders. Furthermore, combining this idea with atom-atom interactions (featuring an intrinsic global SU (N ) symmetry in 173 Yb atoms) results in intriguing possibilities for engineering synthetic quantum systems, as interactions along the synthetic dimension have an intrinsic non-local character (as different sites of the synthetic dimension correspond to the same physical position of space). While effects of atom-atom interactions have already been studied theoretically in combination with the classical gauge potential described above, with the prediction of strongly correlated states that are reminiscent of the fractional quantum Hall effect [170,171,172], the application of this approach to the realisation of other kinds of gauge fields is still to be investigated. Figure 28: The main graph shows the magnitude and direction of the chiral current J vs. synthetic magnetic flux in an experimental configuration similar to that shown in Fig. 27b (with two legs only). The reversal of the edge current above φ = π flux is reminiscent of the change in sign of the Chern number for an extended two-dimensional system (plotted in the inset of the figure as a shading of the Hofstadter butterfly spectrum, with warm colours corresponding to positive Chern number and cold colours corresponding to negative Chern number). Experimental data taken from Ref. [168].

Quantum Simulation of Abelian Gauge
Fields with ultra-cold atoms 6.2.1 Atomic quantum simulator for lattice gauge theories and ring exchange models [173] Reference [173] presents the design of a ring exchange interaction in cold-atomic gases subjected to an optical lattice using well-understood tools for manipulating and controlling such gases. The strength of this interaction can be tuned independently and describes the correlated hopping of two bosons. This design offers the possibility for the atomic quantum simulation of a certain class of strong coupling Hamiltonians and opens an alternative approach for the study of novel and exotic phases with strong correlations. A setup is discussed where this coupling term may allow for the realisation and observation of exotic quantum phases, including a deconfined insulator described by the Coulomb phase of a three-dimensional U(1) lattice gauge theory.

Cold-atom simulation of interacting relativistic quantum field theories[174]
Dirac fermions self-interacting or coupled to dynamic scalar fields can emerge in the low-energy sector of designed bosonic and fermionic cold-atom systems. In the reference [174] this is illustrated with two examples defined in two space-time dimensions: the first one is the self-interacting Thirring model, and the second one is a model of Dirac fermions coupled to a dynamic scalar field that gives rise to the Gross-Neveu model. The proposed cold-atom experiments can be used to probe spectral or correlation properties of interacting quantum field theories thereby presenting an alternative to lattice gauge theory simulations. The Hamiltonians of these systems are supported on one spatial dimension. The necessary building blocks are the Dirac Hamiltonian, which describes relativistic fermions, and the interaction of fermions with themselves or with a dynamic scalar field. The presented models are exactly solvable and serve for demonstrating the ability to simulate important properties of the standard model such as dynamical symmetry breaking and mass generation with cold atoms.

Confinement and lattice QED electric flux-tubes simulated with ultra-cold atoms[175]
The effect of confinement is known to be linked with a mechanism of electric flux tube formation that gives rise to a linear binding potential between quarks. While confinement is a property of non-Abelian gauge theories including QCD, it has been shown that the Abelian model of compact quantum electrodynamics (cQED), also gives rise to similar phenomena. Particularly, it has been shown long ago by Polyakov that in cQED models in 2+1 dimensions, the effect of confinement persists for all values of the coupling strength, in both the strong and the weak coupling regimes, due to non-perturbative effects of instantons [176]. Lattice cQED in 2+1 dimensions, hence provides one of the simplest play grounds to study confinement in a setup that contains some essential properties of full fledged QCD, in a rather simple Abelian system.
It was first suggested in [175] that the well-known gauge invariant Kogut-Susskind Hamiltonian, which describes cQED in 2+1-d, can be simulated by representing each link along the lattice by a localised BEC, properly tuning the atomic scattering interactions, and using external lasers. In this proposal, the cQEDs' degrees of freedom on each link are the condensate's phases which correspond to the periodic cQED vector potentials and the atomic number operator to the quantised electric field on the link. Charged sources are non-dynamical, and introduced to the model by coupling the fields to external static charges. To obtain the Kogut-Susskind Hamiltonian, it is shown that particular two-and four-body interactions between the condensates provide manifest local gauge invariance. Furthermore, to avoid the hopping processes of an ordinary Bose-Hubbard model, one introduces a four-species two-dimensional setup, Raman transitions and two-atom scattering processes in order to obtain particular 'diagonal' hopping and nonlinear interactions. It is then shown that a certain choice of parameters gives rise to gauge invariance in the low-energy sector, hence compact QED emerges. In the strong coupling limit, the atomic system gives rise to electric flux-tubes and confinement. To observe such effects one needs to measure local density de- viations of the BECs. It is reasoned that the effect should persist also outside the perturbatively calculable regime.

Atomic Quantum Simulation of Dynamical Gauge
Fields coupled to Fermionic Matter: From String Breaking to Evolution after a Quench [40] In [40], a quantum simulator for lattice gauge theories is proposed, where bosonic gauge fields are coupled to fermionic matter, which allows the demonstration of experiments for phenomena such as time-dependent string breaking and the dynamics after a quench. Using a Fermi-Bose mixture of ultra-cold atoms in an optical lattice, a quantum simulator is constructed for a U (1) gauge theory coupled to fermionic matter. The construction is based on quantum links which realise a continuous gauge symmetry with discrete quantum variables. At low energies, quantum link models with staggered fermions emerge from a Hubbard-type model which can be quantum simulated. This allows one to investigate string breaking as well as the real-time evolution after a quench in gauge theories, which are inaccessible to classical simulation methods. While the basic elements behind the model have been demonstrated individually in the laboratory, the combination of these tools and the extension to higher dimensions remain a challenge to be tackled in future generations of optical lattice experiments.

Simulating Compact
Quantum Electrodynamics with ultra-cold atoms: Probing confinement and non-perturbative effects [38] An alternative for simulating cQED, that relies on single atoms arranged on a lattice, (rather than small BECs), has been proposed in ref. [38]. In this work the idea is to use single atoms, described in terms of spin-gauge Hamiltonians, in an optical lattice, carrying 2l+1 internal levels. As l increases, it is found that the spin-gauge Hamiltonian manifests a rather fast convergence to the Kogut-Susskind cQED Hamiltonian model. This then enables the simulation in both the strong and weak regime of cQED in 2+1 dimensions. Hence the model can be used to simulate confinement effects in the non-perturbative regime, with a rather compact system and with a modest value of l. The case l = 1, corresponds to the lowest value, which is sufficient for demonstrating confinement and flux tubes between external charges, and is here explicitly constructed. This implementation with single atoms might be experimentally easier compared to the BEC approach [175], as it does not rely on the overlaps of local BECs wave functions and thus involves larger Hamiltonian energy scales.
6.2.6 Quantum simulations of gauge theories with ultra-cold atoms: local gauge invariance from angular momentum conservation [177] Further development in realising compact QED in (1+1) and (2+1)-d is described in [177], which also provides a rather systematic discussion of the structure and needs of quantum simulations of lattice gauge theory in the Hamiltonian form of Kogut and Susskind. In particular, in section IV of the article, a systematic method is described for constructing the required fermion-gauge boson interaction terms needed on the links for the particular cases involving a U (1) interaction, and a Z(N ) elementary interaction. The key point is that one can in fact reduce the problem of gauge invariance to that of conservation of angular momentum in elementary fermion-boson scatterings, by making a clever choice of the internal fermion and boson levels on the vertices and links. Then, while the fermion hops from one vertex to another, it interacts and scatters from bosonic species situated on the link. By using an adequate selection of the internal fermionic and bosonic states one can then guarantee that the resulting interaction is gauge invariant. The general plaquette interaction term tr(U U U † U † ) has been obtained through an auxiliary particle that goes around loops, in section VI of the article. This method can be used for all gauge field interactions of the minimal form ψ † n U ψ n+1 , with U being a precise unitary, by producing such an interaction between an auxiliary particle and the physical matter. Given that that is indeed possible, it is shown how the "loop method" gives rise to plaquette interactions for the cases of U(1) and SU (N ) gauge fields.

Quantum Spin Ice and dimer models with Rydberg atoms[178]
Abelian gauge theories also play a rather prominent role in the theory of frustrated quantum magnets [179]. Refs. [178,180] report two approaches aimed at realising the dynamics of frustrated quantum magnets with direct gauge theoretic interpretation, utilising Rydberg atoms trapped in lattices (either optical or tweezers).
Ref. [178] introduces a toolbox to realise frustrated quantum magnets in two dimensions via Rydberg dressing. The main feature of this type of potentials, which is generated by off-resonantly coupling ground states to Rydberg states (see Fig. 30a), is that, different from the bare Rydberg potential, it exhibits a plateau up to a critical distance of order of the Condon radius, as depicted in Fig. 30b. This feature is already sufficient to generate gauge theory dynamics on a variety of lattices that can be decomposed in triangular unit cells, including Kagome and squagome ones.
In addition, depending on the type of Rydberg states one is coupling to, it is possible to exploit the angular dependence of the Rydberg-Rydberg interactions (encoded in the angular dependent factor A(ϑ)) to engineer Gauss' law constraints in other geometries. For the square lattice case, the resulting Hamiltonian is known as square Ice. The interactions needed to impose the corresponding constraint (the Gauss law, also known as ice rule) are both anisotropic and position dependent. This can be achieved by coupling ground state atoms to p-states: an example of the resulting interaction pattern is depicted in Fig. 30c.
The resulting system dynamics, despite some additional features due to the long-range character of the Rydbergdressing potentials, is able to reproduce the quantum ice rule, including a ground state with resonating valencebond solid order. In addition, even in the absence of quantum fluctuations, an interesting thermal transition to an (imperfect) Coulomb phase takes place. This work has also served as a stimulus for tensor network simulations of twodimensional systems, as reported in Ref. [110].
In Ref. [180] (see also the related work in Ref. [181] in the context of spin-1 models), the methods discussed above were considerably expanded. In particular, it was shown that the full dynamics within the Rydberg manifold can be exploited for the realisation of almost arbitrary spin-spin interactions, including U (1) and Z(2) invariant spin exchanges. Some of these terms could be additionally tuned exploiting quantum interference effects, or local AC Stark shifts. The possibility of utilising Rydberg atoms to engineer constrained models can also be directly applied to dynamics where the matter or gauge fields are integrated out explicitly. For instance, in [182], it was shown how recent experiments in Rydberg atoms arrays have already realised [36] quantum simulations of gauge theories at large scales. In [111], it is shown how to implement a Rydberg-atom quantum simulator to study the non-equilibrium dynamics of an Abelian (1 + 1)-D lattice gauge theory, the implementation locally codifies the degrees of freedom of a Z(3) gauge field, once the matter field is integrated out by means of the Gauss local symmetries.

Toolbox for Abelian lattice gauge theories with synthetic matter[183]
In [183], it is described what can be achieved by using the simplest possible implementation, taking a mixture of two Bose-gases on an optical lattice and working with lattice potentials and their time modulation (shaking). The Hamiltonian governing the system is thus the standard Bose-Hubbard Hamiltonian in which the atoms can hop between neighbouring sites around x gaining an energy J(x) and their interactions lead to an energy cost U . The two species of bosons behave quite differently: a-bosons are strongly interacting, say in the hardcore limit, while b-bosons need to be non-interacting. Their mutual interaction is described by U ab . Upon fast modulation of this inter-species interaction and lattice potentials, an effective Hamiltonian is achieved (for the details of the implementation, refer to [183]), in a Floquet approximation, for the slow degrees of freedom in which the hopping of the a-bosons is modulated in amplitude and phase by the difference of occupations of the b-bosons. This means that the a-bosons behave as charged particles under a vector potential generated by the difference in occupations of b-bosons.
In this setting, an interesting scenario is obtained by considering initially the a-bosons hopping on a dimerised lattice as presented in the left-hand panel of Fig. 31b where in a first approximation they can only hop horizontally from even to odd sites with amplitude J a (h). Thus, at half-filling, there is an insulating phase where each dimerised link contains exactly one delocalised a-boson as in Fig. 31b. By switching on perturbatively both the hopping of b-bosons J b and the vertical hopping of a-bosons J a (v), applying second order perturbation theory with respect to the hopping ratios J 2 b /J a (h) and J 2 a (v)/J a (h) an effective Hamiltonian is achieved that can be described in terms of an exotic gauge theory.
The gauge invariant variables are plaquettes construct from a coarse grained lattice, along the horizontal direction on whose links the following Hamiltonian is obtained, The B p represent the standard plaquette magnetic term and E is the electric field. N represents the rank of the Abelian group Z(N ) and the U (1) theory is obtained in the limit N → ∞. Even in that limit, the above Hamiltonian differs from the standard two-dimensional QED Hamiltonian, since the electric field part of the Hamiltonian acts on neighbouring links, rather than on a single link as in the standard case. This specific pattern of the electric field implies that the low-energy sector of the model can be described by a gas of oriented magnetic dipoles rather than a gas of magnetic monopoles. As a consequence, even when the dipoles condense, they cannot screen the electric field and thus the model contains new exotic deconfined phases. The sketch of the conjectured phase diagram is contained in Fig. 32. 6.2.9 Many-body localisation dynamics from gauge invariance [184] The experimental realisation of a synthetic lattice gauge theory reported in Ref. [47] has immediately stimulated a novel interest in the real-time evolution of low-dimensional lattice gauge theories, following related works in the context of statistical mechanics models [185]. A particularly active area of research in the latter field has been the study of disordered, interacting systems, which, under certain conditions, may fail to thermalise -a phenomenon referred to as many-body localisation (MBL) [186].
In Ref. [184] it was shown that, in contrast to statistical mechanics models, LGTs may display MBL even in the absence of any disorder. The feature at the basis of this phenomenon is the presence of Hilbert space sectors, which are a characteristic feature in LGTs. For Abelian theories, these subspaces are simply characterised by a given background charge distribution.
An arbitrary translationally invariant state is typically spanning several of these super-selection sectors: this implies that the resulting time evolution is actually sensitive to several, distinct background charge distributions, mimicking the time evolution of a system under an inhomogeneous potential. This does not map into uncorrelated disordered patterns, at least for the gauge group and the class of initial states considered in Ref. [184].
This type of dynamics was then analysed numerically using large-scale exact diagonalisation methods for the lattice Schwinger model, starting from initial states with gauge fields in an equal weight superposition of E = (−1, 0, 1), and staggered fermions. A sample of these results is depicted in Fig. 33: absence of relaxation for two local observables is indicated in the upper two panels. The bottom two report a finite-size scaling analysis for a thermalising and non-thermalising case. The entanglement evolution of this MBL-type dynamics is however rather peculiar, as signalled by the evolution of the half-chain entanglement entropy: while this typically increases logarithmically with time, in the present context, a double-logarithmic increase was observed up to numerically accessible time-scales.
A related phenomenology was also observed in theories without confinement in Ref. [187], and was also proposed as an implementation route for non-interacting models in Ref. [188] 6.2.10 Discretizing Quantum Field Theories [189,190] The majority of platforms for quantum simulations: ultracold atoms in optical lattices, Rydberg atoms in traps or trapped ions deal with discretised versions of the considered Quantum Field Theory. As always in such cases, this can be viewed as a nuisance or an opportunity. In fact, to simulate LGT, discrete lattice models and their implementations are needed. But, studying and quantum simulating discrete versions of continuous QFT models might also lead to new fascinating physics. In [190] a Gross-Neveu-Wilson model was studied and its correlated symmetryprotected topological phases revealed. A Wilson-type discretisation of the Gross-Neveu model, a fermionic N -flavour quantum field theory displaying asymptotic freedom and chiral symmetry breaking, can serve as a playground to explore correlated symmetry-protected phases of matter using techniques borrowed from high-energy physics. A large-N study, both in the Hamiltonian and Euclidean formalisms, was used and analysed.
In [189] renormalisation group flows for Wilson-Hubbard matter and the topological Hamiltonian were studied. The aim was to understand the robustness of topological phases of matter in the presence of interactions, a problem that poses a difficult challenge in modern condensed matter physics, showing interesting connections to high-energy physics (see also other works facilitating quantum simulations with ultra-cold atoms and ions: [191] for exploration of interacting topological insulators with ultra-cold atoms in the synthetic Creutz-Hubbard model, [192] for symmetry-protected topological phases in lattice gauge theories, [193] for the study of the topological Schwinger model, [194] for Z(N ) gauge theories coupled to topological fermions. These connections present an analysis of the continuum long-wavelength description of a generic class of correlated topological insulators of Wilson-Hubbard type, feasible for experiments with quantum simulators.

6.2.11
Interacting bosons on dynamical lattices [195,196,194,197] For ultra-cold atoms in optical lattices quantum simulations with bosons are more experimentally friendly than those with fermions. For these reasons there is a clear tendency in recent years to design and study lattice models in which bosons replace fermions. Particularly fruitful and fascinating are such models, dubbed as dynamical lattices, where some objects live on the lattice links, like in LGTs. Inspired by the so-called fluctuating bond superconductivity for fermionic Hubbard-like models that include quantised phonons on the links (see for instance [198]), a bosonic version of such models was formulated where phonons on the links are replaced by two states of a spin-1/2, which provides a minimal description of a dynamical lattice. These Z(2) Bose-Hubbard models are extraordinarily rich and lead to bosonic Peierls transitions [195], intertwined topological phases [196,199], symmetry-protected topological defects [194], etc. Thus even in the absence of the gauge invariance, these models allow one to explore interesting strongly-correlated topological phenomena in atomic systems. Similar models, now with exact Z(2) gauge symmetry, were studied theoretically [200] and realised recently in a cold-atom experiment [201]. Peierls states and phases were also studied in models of Floquet-engineered vibrational dynamics in a two-dimensional array of trapped ions [202]. Recently, density-dependent Peierls phases were realised coupling dynamical gauge fields to matter [203].
The culmination of this effort is the paper [197], in which the bosonic Schwinger model was studied and confinement and lack of thermalisation after quenches was observed. The vacuum of a relativistic theory of bosons coupled to a U (1) gauge field in 1+1 dimensions (bosonic Schwinger model) was excited out of equilibrium by creating a spatially separated particle-anti-particle pair connected by an electric flux string. During the evolution, a strong confinement of the bosons is observed witnessed by the bending of their light cone, reminiscent of what was observed for the Ising model. As a consequence, for the time scales amenable to simulations, the system evades thermalisation and generates exotic asymptotic states. These states extend over two disjoint regions, an external deconfined region that seems to thermalise, and an inner core that reveals an area-law saturation of the entanglement entropy. Qualitative phase diagram of the Hamiltonian H gauge in the g-N plane. The upper shaded region denotes a gapped phase in the strong coupling limit. In the weak coupling limit, for low N , the system is gapped but deconfined (region shaded in light blue). In the U (1) limit, the system becomes gapless and an exotic dipole liquid phase emerges (region shaded in dark blue). In the intermediate region, the system is effectively in a one-dimensional gapless Bose liquid phase, extracted from [183].

Abelian Quantum Simulation with trapped ions and superconducting circuits
6.3.1 Quantum simulation of a lattice Schwinger model in a chain of trapped ions [204] Ref. [204] represents the first attempt to identify gauge theory dynamics in trapped ion systems. From the point of view of analog quantum simulation, the main difference between atom and ion experiments resides in the type of degrees of freedom one can harness. While in former settings one is typically dealing with itinerant fermions and bosons, in the latter, the dynamics is dictated by the interactions between the ions' internal degrees of freedom (either nuclear or electronic spin) and the phonon modes generated by the ion crystal.
The internal levels of the ions and the phonon modes are then coupled via external light fields: upon integration of the phonon degrees of freedom the resulting dynamics is then given by a spin chain, typically with S = 1/2. Ref. [204] exploits the fact that such interactions can be made spatially anisotropic. The main idea is that one can identify ions on one sub-lattice (A) as 'matter' fields, and ions on the other sub-lattice (B) as gauge fields in the quantum link formulation discussed above.
Within this setting, gauge invariance is then enforced by energy punishment. The latter is generated using a combination of two laser-assisted (e.g., Raman) interactions between qubits, that can be made spatially inhomogeneous either utilising local shifts, or by properly shaping the Rabi frequency of the light beams. The effective quan-3 sector, we analytically integrate out the gauge fields [41]. We assumeL 0 = 0 to minimize boundary effects, and apply a Jordan-Wigner transformations to the fermionic fields in order to re-cast the dynamics as a spin model, † n n = ( z n + 1)/2. The gauge fields can be sequentially integrated out by noting that which simply describes the fact that, given the value of the ingoing electric field on the left side of a site, and given the values of the dynamical and static charge, the value of the outgoing electric field is unambiguously fixed. After integration, the resulting Hamiltonian dynamics crucially depends on {q ↵ } -that is, states in different superselection sectors will evolve according to different Hamiltonians H {q↵} . This is a direct consequence of the fact that, because of Gauss law, different superselection sectors describe dynamics subject to different static charge configurations. In each sector, the corresponding Hamiltonian is made of two contributions, In . The first one describes electrongauge coupling, which is not sensitive to background charges, and is given by: The second term originates from the electric field potential term, which is now a function of the fermionic populations only and reads: This term can be conveniently decoupled into two parts, one containing two-body terms and one containing only one-body terms. The former results in: Crucially, this last part of the Hamiltonian depends explicitly on the superselection sector via {q`}. As such, starting from initial states of the form in Eq. (5), the system dynamics is dictated by a charge distribution average, that is: which is effectively describing a disorder average, since the terms in H {q↵} Z effectively act as a (correlated) disorder for the spin dynamics. The observation in Eq. (11) is applicable to arbitrary Abelian Wilson theories, and even non-Abelian QLM, in one-dimensional systems. The same reasoning can be applied to two-dimensional systems, in particular, to the evolution of quenched gauge theories. It explicitly shows that, contrary to conventional spin models, the dynamics of systems endowed by a gauge symmetry can naturally lead to phenomena related to disordered systems, even in the absence of any disorder (both on the initial state, and in the dynamics).
Many-body localization dynamics in U(1) lattice gauge theories. -We now turn to a numerical investigation of the system dynamics described by Eq. (11). We address this by employing a computationally optimized approach to the method of Krylov subspaces -for details, see Ref. [42].
As a first figure of merit, we focus on the staggered occupation of the fermions: which, by a Jordan-Wigner transformation, can be expressed as ⌫(t) = 1 2N P N n=1 h( 1) n z n + 1i by transforming fermionic fields to spin operators. Since we employ tum link model dynamics is then generated in second order perturbation theory, similar to several atomic schemes.
Here the main source of experimental imperfections is the presence of long-ranged terms in the spin-spin interactions. While the latter are gauge invariant, their presence might be detrimental to observe physical phenomena: as such, one has to find a balance between enforcing the constraint, while still keeping the effective QLM dynamics sufficiently fast. In [204], this question was addressed in the context of ground state preparation: there, it is possible to find an optimal parameter regime that allows the observation of the Ising transition present in the QLM.
Following this first work, other proposals have been presented to realise LGT dynamics in trapped ion systems. Ref. [205] dealt instead with two-dimensional models, mostly focusing on condensed matter realisations of Z(2) quenched gauge theories which are known to exhibit topological quantum spin liquid behaviour. The main idea there was to utilise localised phonon modes generated in two-dimensional lattices either via Rydberg excitations, or by laser-pinning a subset of the trapped ions. At the few plaquette level, the model corresponds to the frus-

Superconducting Circuits for Quantum Simulation of Dynamical Gauge Fields[207]
The essential building blocks of a superconducting quantum simulator are described in [207] for dynamical lattice gauge field theories, where the basic physical effects can already be analysed with an experimentally available number of coupled superconducting circuits. This proposal analyses a one-dimensional U (1) quantum link model, where superconducting qubits play the role of matter fields on the lattice sites and the gauge fields are represented by two coupled microwave resonators on each link between neighbouring sites. A detailed analysis of a minimal experimental protocol for probing the physics related to string breaking effects shows that, despite the presence of decoherence in these systems, distinctive phenomena from condensed-matter and high-energy physics can be visualised with state-of-the-art technology in small superconducting-circuit arrays.

Loops and Strings in a Superconducting Lattice
Gauge Simulator [208] An architecture for an analog quantum simulator of electromagnetism in 2 + 1 dimensions is proposed in [208], based on an array of superconducting fluxonium devices. The simulator can be tuned between intermediate and strong coupling regimes, and allows non-destructive measurements of non-local, space-like order and disorder parameters, resolving an outstanding gap in other proposals. Moreover, a physical encoding of the states is provided, where local electric field terms are non-trivial. The devices operate in a finite-dimensional manifold of lowlying eigenstates, to represent discrete electric fluxes on the lattice. The encoding is in the integer (spin 1) representation of the quantum link model formulation of compact U(1) lattice gauge theory. In the article, it is shown how to engineer Gauss' law via an ancilla mediated gadget construction, and how to tune between the strongly coupled and intermediately coupled regimes. The protocol is rather robust to inhomogeneities, allowing for implementations in superconducting arrays, and numerical evidence is presented that lattice QED in quasi-(2 + 1) dimensions exhibits confinement. Beyond ground state characterisation, the simulator can be used to probe dynamics and measure the evolution of non-local order or disorder parameters. The witnesses to the existence of the predicted confining phase of the model are provided by non-local order parameters from Wilson loops and disorder parameters from 't Hooft strings. In [208], it is shown how to construct such operators in this model and how to measure them non-destructively via dispersive coupling of the fluxonium islands to a microwave cavity mode. Numerical evidence was found for the existence of the confined phase in the ground state of the simulation Hamiltonian on a ladder geometry. Using ultra-cold alkaline-earth atoms in optical lattices, [43] constructs a quantum simulator for U (N ) and SU (N ) lattice gauge theories with fermionic matter based on quantum link models. These systems share qualitative features with QCD, including chiral symmetry breaking and restoration at non-zero temperature or baryon density. The proposal builds on the unique properties of quantum link models with rishons representing the gauge fields: this allows a formulation in terms of a Fermi-Hubbard model, which can be realised with multi-component alkaline-earth atoms in optical lattices, and where atomic physics provides both the control fields and measurement tools for studying the equilibrium and non-equilibrium dynamics and spectroscopy.
6.4.2 A cold-atom quantum simulator for SU (2) Yang-Mills lattice gauge theory [41] A realisation of a non-Abelian lattice gauge theory, which is an SU (2) Yang-Mills theory in (1+1)-dimensions is pro- Figure 35: U (1) quantum link model engineered in a fluxonium array. (a)"Electric"Ê α;β , and "magnetic'Û α;β , degrees of freedom are associated with links < α; β > of a square lattice. The link degrees of freedom (red circles) are encoded in eigenstates of the fluxonia. The ancillae (blue diamonds) on vertices are inductively coupled to neighbouring link islands to mediate the Gauss constraint and plaquette interactions are obtained via link nearestneighbour capacitive coupling. (b) Superconducting circuit elements used to build and couple components of the simulation. The link devices have local phaseφ link and capacitive, inductive, and flux-biased Josephson energies E C , E L , and E J , respectively, and similarly for the ancilla devices. The capacitive and inductive coupling energies are E c C and E c L . (c) A minimal quasi-1D "ladder" implementation embedded in a microwave cavity (black box), in which a 't Hooft string of link fluxonia (green circles) can be measured via anancilla coupled to the cavity (green triangle). Figure taken from Ref. [208]. posed in [41]. The system one wants to simulate involves a non-Abelian gauge field and a dynamical fermionic matter field. As in ordinary lattice gauge theory the fermions are located at the vertices and the gauge fields on the links. Using staggered fermion methods, in the (1+1)-d case the Hamiltonian is equivalent, to the non-Abelian Schwinger model, with a minimal-coupling interaction of the form ψ † n U n ψ n+1 , involving a 2 × 2 unitary matrix. In order to simulate the non-Abelian SU (2) Hamiltonian, this paper uses a particular realisation of the group elements and "left" and "right" generators, using a Jordan-Schwinger map, that connects harmonic oscillators (bosons) and angular momentum. Mapping SU (N ) to a bosonic system allows one to express the gauge fields by means of bosonic atoms in the prepotential method [209]. For the SU (2) group, one then needs four bosonic species for each link.
Remarkably, gauge invariance arises as a consequence of conservation of angular momentum conservation and thus is fundamentally robust. However, the effective Hamiltonian obtained here is not valid beyond the strong coupling limit. 6.4.3 Constrained dynamics via the Zeno effect in quantum simulation: Implementing non-Abelian lattice gauge theories with cold atoms [210] Implementing non-Abelian symmetries poses qualitatively new challenges from the theory side with respect to the Abelian case. In particular, energy punishment strategies, which are widespread for U(1) theories, are not immediately viable. The main reason is the following: consider a set of local generators G α x , and a microscopic Hamiltonian of the type: where H 0 is a gauge invariant term, and H 1 is gauge variant. Here, differently from the Abelian case, where a single generator is considered, one deals with several constraints, that have to be satisfied in a symmetric manner -that is, U α x = U x . Typically, this requires fine-tuning, making energy punishment strategies not immediately viable.
In Ref. [210], an alternative procedure was proposed based on quantum Zeno dynamics [211]. The basic idea is to consider a set of classical noise sources ξ α x (t) coupled to the generators of the gauge symmetry, described by the effective microscopic Hamiltonian: In the limit of independent, white noise sources, the system dynamics is described by a master equation, whose corresponding effective Hamiltonian reads: In the large noise limit, where κ is much larger than all microscopic scales involved in H 1 , the effective dynamics is constrained to the gauge invariant subspace H P . In this limit, any coupling term between the gauge invariant subspace to the gauge variant one H U is suppressed by the noise terms (see Fig. 37a). Depending on the nature of H 1 , one can use a limited number of noise sources, with the condition that neighbouring blocks are subject to distinct noise terms, as illustrated in Fig. 37b. From a theoretical viewpoint, this scheme can be seen as a classical analogue of the quantum Zeno effect [211], which has also been discussed in the context of quantum many-body systems as a source of local interaction [212,213,214].
The main feature of this scheme is that, different from energy penalty schemes, here the constraints are induced by just coupling to simple operator terms (there is no microscopic term (G α x ) 2 that shall be engineered), that can be easily controlled by means of external fields. This drastically simplifies the conditions required to achieve gauge invariant dynamics. This comes at the price of needing to realise the desired dynamics H 0 without the help of perturbation theory: two applications in the context of both Abelian and non-Abelian theories are discussed in Ref. [210] in the context of cold atoms in optical lattices. 6.4.4 SO(3) "Nuclear Physics" with ultra-cold Gases [215] An ab initio calculation of nuclear physics from QCD, the fundamental SU (3) gauge theory of the strong interaction, remains an outstanding challenge. In this proposal, the emergence of key elements of nuclear physics using an SO(3) lattice gauge theory as a toy model for QCD is discussed. This model is accessible to state-of-the-art quantum simulation experiments with ultra-cold atoms in an optical lattice. The phase diagram of the model is investigated, and showed that it shares some fundamental properties with QCD, most prominently confinement, the spontaneous breaking of chiral symmetry, and its restoration at finite baryon density, as well as the existence of few-body bound states, which are characteristic features of nuclear physics. The most critical step in implementing a gauge theory on a quantum simulator is to make sure that the gauge symmetry remains intact during the simulation. It is shown how the lattice gauge model, which has a non-Abelian gauge symmetry, can be realised in a quantum simulator platform by encoding the operators directly in the gauge invariant subspace, thus guaranteeing exact gauge invariance. This encoding strategy is generally applicable to the whole class of quantum link models, which are extensions of Wilson's formulation of lattice gauge theories. Quantum link models permit encoding to local spin Hamiltonians. This not only makes the implementation feasible on different platforms, such as ultracold atoms and molecules trapped in optical lattices, but also establishes a novel connection between non-Abelian lattice gauge theories including matter fields and quantum magnetism.
Concretely, it is shown in [215] how (1 + 1)-d SO(3) lattice gauge theories can be naturally realised using ultracold atoms in optical lattices. The phase diagram of these models features several paradigmatic phenomena, e.g. the presence of stable two-body bound states, phases where chiral symmetry is spontaneously broken, where the broken chiral symmetry is restored at finite baryon density, and emergent conformal invariance. The dynamics in the gauge invariant sector can be encoded as a spin S = 3/2 Heisenberg model, i.e., as quantum magnetism, which has a natural realisation with bosonic mixtures in optical lattices, and thus sheds light on the connection between non-Abelian gauge theories and quantum magnetism. This encoding technique has the dramatic advantage of realising gauge invariance exactly, and at the same time bypassing the complex interaction engineering which is required for non-Abelian theories.

Conclusions
At the time of writing this review, there is a coordinated effort to study the possible applications of quantum technologies to the study of gauge theories. The global aim is to develop novel methods and techniques, namely classical and quantum hardware and software, that eventually could be applied to study open problems in different fundamental and applied fields of science ranging from materials science and quantum chemistry to astrophysics and that will impact fundamental research and our everyday life. This challenge is a highly collaborative advanced multidisciplinary science and cutting-edge engineering project with the potential to initiate or foster new lines of quan- tum technologies. Along the way, the scientific community that is growing around this topic, is developing a deeper fundamental and practical understanding of systems and protocols for manipulating and exploiting quantum information; at enhancing the robustness and scalability of quantum information processing; identifying new opportunities and applications fostered through quantum technologies and enhancing interdisciplinarity in crossing traditional boundaries between disciplines in order to enlarge the community involved in tackling these new challenges.
The results of this endeavour serve as benchmarks for the first generation of quantum simulators and will have far reaching consequences, e.g., in the long run this research will enable the study and design of novel materials with topological error correcting capabilities, which will play a central role in the quest for building a scalable quantum computer. In particular, in the review, the most advanced quantum simulation of a lattice gauge theory achieved so far is presented, a digital ion-trap quantum variational optimisation applied to find the ground state of the Schwinger model. Moreover, it is also reviewed how novel tensor network methods are being developed and applied to study such systems in one and two dimensions, and finally, how new directions are being proposed for quantum simulations of more complex theories.
Lattice gauge theories provide both motivation and a framework for pushing forward the interdisciplinary advancement of quantum technologies. While the quantum simulation of classical intractable aspects of QCD (such as its real-time evolution or its phase diagram at high baryon density) remain long-term goals with a potentially large impact on particle physics, a wide class of lattice gauge theories, often with applications in condensed matter physics or quantum information science, suggests itself for theoretical investigation and experimental realisation. There is a lot of interesting physics to be discovered along the way towards developing powerful hard-and software for the fast growing field of quantum simulation and computation technology. All the authors acknowledge the participation in the EU-QUANTERA project QTFLAG.

Authors contributions
All the authors were involved in the preparation of the manuscript. All the authors have read and approved the final manuscript.