## Abstract

Gauge theories are of paramount importance in our understanding of fundamental constituents of matter and their interactions. However, the complete characterization of their phase diagrams and the full understanding of non-perturbative effects are still debated, especially at finite charge density, mostly due to the sign-problem affecting Monte Carlo numerical simulations. Here, we report the Tensor Network simulation of a three dimensional lattice gauge theory in the Hamiltonian formulation including dynamical matter: Using this sign-problem-free method, we simulate the ground states of a compact Quantum Electrodynamics at zero and finite charge densities, and address fundamental questions such as the characterization of collective phases of the model, the presence of a confining phase at large gauge coupling, and the study of charge-screening effects.

## Introduction

Ranging from high-energy particle physics (Standard Model)^{1,2,3} to low-temperature condensed matter physics (spin liquids, quantum Hall, and high-*T*_{c} superconductivity)^{4,5}, gauge theories constitute the baseline in our microscopical description of the universe and are a cornerstone of contemporary scientific research. Yet, capturing their many-body body behavior beyond perturbative regimes, a mandatory step before experimentally validating these theories often eludes us^{6}. One for all, the quark confinement mechanism in quantum chromodynamics, a founding pillar of the Standard Model which has been studied for almost half a century, is still at the center of current research efforts^{7,8,9,10,11,12}. Indeed, a powerful numerical workhorse such as Monte-Carlo simulations^{13,14,15}, capable of addressing discretized lattice formulations of gauge theories^{11,16,17,18}, struggles in highly interesting regimes, where matter fermions and excess of charge are concerned, due to the infamous sign problem^{19}. In recent years, a complementary numerical approach, Tensor Networks (TN) methods, have found increasing applications for studying low-dimensional Lattice Gauge Theories (LGT) in the Hamiltonian formulation^{20,21}. As tailored many-body quantum state ansätze, TNs are an efficient approximate entanglement-based representation of physical states, capable of efficiently describe equilibrium properties and real-time dynamics of systems described by complex actions, where Monte Carlo simulations fail to efficiently converge^{22}. TN methods have proven remarkable success in simulating LGTs in (1+1) dimensions^{23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41}, and very recently they have shown potential in (2+1) dimensions^{42,43,44,45,46,47,48,49,50}. To date, due to the lack of efficient numerical algorithms to describe high-dimensional systems via TNs, no results are available regarding the realistic scenario of LGTs in three spatial dimensions.

In this work, we bridge this gap by numerically simulating, via TN ansatz states, an Abelian lattice gauge theory akin to (3 + 1) Quantum Electrodynamics (QED), at zero temperature. We show that, by using the quantum link formalism (QLM) of LGTs^{51,52} and an unconstrained Tree TN (TTN), we can access multiple equilibrium regimes of the model, including finite charge densities. Precisely, we analyze the ground state properties of quantum-link QED in (3 + 1)D for intermediate system sizes, up to 512 lattice sites. The matter is discretized as a staggered spinless fermion field on a cubic lattice^{16}, while the electromagnetic gauge fields are represented on lattice links, and truncated to a compact representation of spin-*s*. Here, we present results from a nontrivial representation for lattice gauge fields (the spin-1 case), with possible generalizations to higher spin requiring only a polynomial overhead in *s*. Our picture can be similarly adapted to embed non-Abelian gauge symmetries, such as they appear in QCD^{32}. Finally, we stress that the truncation of the gauge field is a common step in quantum simulations and computations^{53,54,55,56,57,58,59,60,61,62}, making the presented numerical approach a landmark benchmarking and cross-verification tool for current and future experiments. By variationally approximating the lattice QED ground state with a TTN, we address a variety of regimes and questions inaccessible before. In the scenario with zero excess charge density, we observe that the transition between the vacuum phase and the charge-crystal phase is compatible with a second-order quantum phase transition^{47}. In the limit of zero magnetic couplings, this transition occurs at negative bare masses *m*_{0}, but as the coupling is activated, the critical point is shifted to larger, and even positive, *m*_{0} values. To investigate field-screening properties, we also consider the case where two parallel charged plates are placed at a distance (a capacitor). By studying the polarization of the vacuum in the inner volume, we observe an equilibrium string-breaking effect akin to the Schwinger mechanism. Furthermore, we address the confinement problem by evaluating the binding energies of charged particle pairs pinned at specified distances. Finally, we consider the scenario with a charge imbalance into the system, i.e., at finite charge density, and we characterize a regime where charges accumulate at the surface of our finite sample, analogously to a classic perfect conductor.

## Results

### The model

Hereafter, we numerically simulate, at zero temperature, the Hamiltonian of *U*(1) quantum electrodynamics on a finite *L* × *L* × *L* three-dimensional simple cubic lattice^{16,17,18}

with \({\bf{x}}\equiv \left(i,j,k\right)\) for 0 ≤ *i*, *j*, *k* ≤ *L* − 1 labeling the sites of the lattice and \({\hat{\square }}_{{\mu }_{\alpha },{\mu }_{\beta }}={\hat{U}}_{{\bf{x}},{\mu }_{\alpha }}{\hat{U}}_{{\bf{x}}+{\mu }_{\alpha },{\mu }_{\beta }}{\hat{U}}_{{\bf{x}}+{\mu }_{\beta },{\mu }_{\alpha }}^{\dagger }{\hat{U}}_{{\bf{x}},{\mu }_{\beta }}^{\dagger }\). Here, we adopted the Kogut–Susskind formulation^{16}, representing fermionic degrees of freedom with a staggered spinless fermion field \(\{{\hat{\psi }}_{{\bf{x}}},{\hat{\psi }}_{{\bf{x}}^{\prime} }^{\dagger }\}={\delta }_{{\bf{x}},{\bf{x}}^{\prime} }\) on lattice sites. Their bare mass *m*_{x} = (−1)^{x}*m* is staggered, as tracked by the site parity (−1)^{x} = (−1)^{i+j+k}, so that fermions on even sites represent particles with positive electric charge +*q*, while holes on odd sites represent anti-particles with a negative charge −*q*, as shown in Fig. 1. Charge \(\hat{Q}\) conservation is thus expressed as global fermion number \(\hat{N}\) conservation, since \(\hat{Q}={\sum }_{{\bf{x}}}\left({\hat{\psi }}_{{\bf{x}}}^{\dagger }{\hat{\psi }}_{{\bf{x}}}-\frac{1-{(-1)}^{{\bf{x}}}}{2}\right)=\hat{N}-{L}^{3}/2\).

The links of the 3D lattice are uniquely identified by a couple of parameters (**x**, *μ*) where **x** is any site, *μ* is one of the three positive lattice unit vectors \({\mu }_{x}\equiv \left(1,0,0\right)\), \({\mu }_{y}\equiv \left(0,1,0\right)\), and \({\mu }_{z}\equiv \left(0,0,1\right)\). The gauge fields are defined on lattice links through the pair of operators \({\hat{E}}_{{\bf{x}},\mu }\) (electric field) and \({\hat{U}}_{{\bf{x}},\mu }\) (unitary comparator) that satisfy the commutation relation

For comfort of notation, we can extend the definition to negative lattice unit vectors via \({\hat{E}}_{{\bf{x}}+\mu ,-\mu }=-{\hat{E}}_{{\bf{x}},\mu }\) and \({\hat{U}}_{{\bf{x}}+\mu ,-\mu }={\hat{U}}_{{\bf{x}},\mu }^{\dagger }\).

The Hamiltonian of Eqs. (1a)–(1c) consists of four terms: the parallel transporter (1a) describes the creation and annihilation of a particle–antiparticle pair, shifting the gauge field in-between to preserve local gauge symmetries. The staggered mass and the electric energy density (1b) are completely local. Finally, the plaquette terms (1c) capture the magnetic energy density and are related to the smallest Wilson loops along the closed plaquettes along the three planes *x–y*, *x*–*z*, *y*–*z* of the lattice. In dimensionless units (*ℏ* = *c* = 1), the couplings in Eqs. (1a)–(1c) are not independent: They can be expressed as *t* = 1/*a*, *m* = *m*_{0}, \({g}_{{\rm{e}}}^{2}={g}^{2}/a\), \({g}_{{\rm{m}}}^{2}=8/({g}^{2}a)\), where *a* is the lattice spacing, *g* is the coupling constant of QED, and *m*_{0} is the bare mass of particles/antiparticles. The numerical setup allows us to consider the couplings (*t*, *m*, *g*_{e}, *g*_{m}) as mutually independent. We then recover the physical regime of QED by enforcing \({g}_{{\rm{e}}}{g}_{{\rm{m}}}=2\sqrt{2}t\)^{16}. We also fix the energy scale by setting *t* = 1.

The local *U*(1) gauge symmetry of the theory is encoded in Gauss’s law, whose generators

are defined around each lattice site **x**. The sum in Eq. (3) involves the six electric field operators on the links identified by ±*μ*_{x}, ±*μ*_{y}, ±*μ*_{z}. Each \({\hat{G}}_{{\bf{x}}}\) commutes with the Hamiltonian \(\hat{H}\). In the absence of static (background) charges, the gauge-invariant Hilbert space consists of physical many-body quantum states \(\left|{{\Phi }}\right\rangle \) satisfying \({\hat{G}}_{{\bf{x}}}\left|{{\Phi }}\right\rangle =0\) at every site **x**.

As stressed in the standard Wilson’s formulation of lattice QED^{11}, faithful representations of the \((\hat{{\bf{E}}},\hat{{\bf{U}}})\) algebra are infinite-dimensional. A truncation to a finite dimension becomes therefore necessary for numerical simulations with TN methods, which require a finite effective Hilbert dimension at each lattice site. We use the quantum link model (QLM) approach in which the gauge field algebra is replaced by *S**U*(2) spin algebra, i.e., \({\hat{E}}_{{\bf{x}},\mu }\equiv {\hat{S}}_{{\bf{x}},\mu }^{z}\) and \({\hat{U}}_{{\bf{x}},\mu }\equiv {\hat{S}}_{{\bf{x}},\mu }^{+}/s\) for a spin-*s* representation. This substitution keeps the electric field operator hermitian and preserves Eq. (2), but \(\hat{{\bf{U}}}\) is no longer unitary. Throughout this work, we will select *s* = 1, the smallest representation ensuring a nontrivial contribution of all the terms in the Hamiltonian (see also Fig. 1). This truncation introduces a local energy cutoff based on \({g}_{{\rm{e}}}^{2}\), which in turn requires larger spin *s* to accurately represent weaker coupling regimes, still potentially accessible via TNs^{31}.

### Transition at zero charge

We focus on the zero charge sector, i.e., \({\sum }_{{\bf{x}}}{\hat{\psi }}_{{\bf{x}}}^{\dagger }{\hat{\psi }}_{{\bf{x}}}=\frac{{L}^{3}}{2}\), and periodic boundary conditions. As shown in Fig. 2 (upper panel), for \({g}_{{\rm{m}}}^{2}=0\) the system undergoes a transition between two regimes, analogously to the (1 + 1)D and (2 + 1)D cases^{25,37,47}: for large positive masses, the system approaches the bare vacuum, while for large negative masses, the system is arranged into a crystal of charges, a highly degenerate state in the semiclassical limit (*t* → 0) due to the exponential number of electric field configurations allowed. We track this transition by monitoring the average matter density \(\rho =\frac{1}{{L}^{3}}{\sum }_{{\bf{x}}}<{\rm{GS}}| {\hat{n}}_{{\bf{x}}}| {\rm{GS}}> \), where \({\hat{n}}_{{\bf{x}}}=\frac{1+{(-1)}^{{\bf{x}}}}{2}-{(-1)}^{{\bf{x}}}{\hat{\psi }}_{{\bf{x}}}^{\dagger }{\hat{\psi }}_{{\bf{x}}}\) is the matter occupation operator and the many-body ground state \(\left|{\rm{GS}}\right.> \) has been computed by TTN algorithm (see the “Methods” section for details). Figure 2b displays the result for different sizes *L* (and \({g}_{{\rm{e}}}^{2}/2=t=1\)), portraying the transition. Panels (a) and (c) display local configurations of matter \(<{\hat{n}}_{{\bf{x}}}> \) and gauge sites \(\langle {\hat{E}}_{{\bf{x}},\mu }\rangle \) for *m* = −3.0 and *m* = +3.0 respectively. In the former regime, the algorithm seems to favor a single allowed configuration of gauge fields rather than a superposition of many configuations: This is due to the fact that, when \({g}_{{\rm{m}}}^{2}=0\), the matrix element that rearranges the configurations occurs at very high perturbative order in ∣*t*/*m*∣, and is numerically neglected. A finite-size scaling analysis of the transition (as detailed in the Methods’ subsection “Critical points: scaling analysis”) yields results compatible with a II-order phase transition, with the critical point occurring at negative bare masses *m*.

The same transition appears to be more interesting when we activate the magnetic coupling, by setting \({g}_{{\rm{m}}}^{2}=8{t}^{2}/{g}_{{\rm{e}}}^{2}=4\) (physical line). The phase at large negative *m* now appears to be a genuine superposition of many configurations of the electric field, as they are coupled by matrix elements of the order \(\sim\! {g}_{{\rm{m}}}^{2}\), kept as numerically relevant by the algorithm. Moreover, the transition is still compatible with an II-order phase transition, and the critical point is shifted to larger *m* values. This can lead to a critical bare mass *m*_{c} that is positive (as we observed *m*_{c} ≈ + 0.22 for the case \({g}_{{\rm{e}}}^{2}/2=t=1\)), ultimately making the transition physically relevant.

### Quantum capacitor

To investigate field-screening and equilibrium string-breaking properties, we analyze the scenario where two charged plates (an electric capacitor) are placed at the opposite faces of a volume, with open boundary conditions (OBC). In our simulations, we achieve this regime by setting large local chemical potentials on the two boundaries. We expect that for small positive masses *m*, the vacuum inside the plates will spontaneously polarize to an effective dielectric, by creating particle and antiparticle pairs to screen the electric field from the plates, into an energetically favorable configuration.

We observe this phenomenon by monitoring the charge density function along the direction *μ*_{x} orthogonal to the plates \({q}_{{\rm{c}}}(d)=\frac{2}{{L}^{2}}\mathop{\sum }\nolimits_{j,k = 1}^{L}<{\rm{GS}}| {(-1)}^{{\bf{x}}}{\hat{\psi }}_{(d,j,k)}^{\dagger }{\hat{\psi }}_{(d,j,k)}| {\rm{GS}}> \) as well as the electric field amplitude along *μ*_{x}, \({E}_{{\bf{x}},{\bf{x}}+{\mu }_{x}}^{{\rm{c}}}(d)=\frac{2}{{L}^{2}}\mathop{\sum }\nolimits_{j,k = 1}^{L}<{\rm{GS}}| {\hat{E}}_{(d,j,k),(d+1,j,k)}| {\rm{GS}}> \), as presented in Fig. 3.

A transition from a vacuum regime to a string-breaking dielectric regime is observed, when driving *m* from negative to positive. However, here the critical point occurs at positive masses (*m*_{c} > 0) even at zero magnetic coupling \({g}_{{\rm{m}}}^{2}=0\), analogously to the (1 + 1)D case^{25}. In conclusion, the charged capacitor can make the phase transition physical even when *g* can not be tuned.

The observed behavior can be interpreted as an equilibrium counterpart to the *Schwinger mechanism*, a real-time dynamical phenomenon in which the spontaneous creation of electron-positron pairs out of the vacuum is stimulated by a strong external electric field^{63}. This could either be potentially verified in experiments or quantum simulations, by means of adiabatic quenches, ramping up the capacitor voltage.

### Confinement properties

The (3 + 1)-dimensional pure compact lattice QED predicts a confining phase at large coupling *g*^{11,64,65,66,67}. This phase, where the magnetic coupling is negligible, is characterized by the presence of a linear potential between static test charges, and is expected to survive at the continuum limit. By decreasing *g*, the system undergoes a phase transition to the Coulomb phase where the magnetic terms are not negligible and the static charges interact through the 1/*r* Coulomb potential at distance *r*^{68}. When the gauge field is coupled to dynamical matter (*t* ≠ 0 and finite *m*), new possible scenarios emerge, such as the string-breaking mechanism. Nevertheless, the transition between confined and deconfined phases is still expected to occur^{69}.

We can investigate this specific scenario with our TN method: we consider a 16 × 4 × 4 lattice and pin two opposite charges via large local chemical potentials at distance *r* along direction *μ*_{x}. The energy *E*(*r*) = *V*(*r*) − *V*(*∞*) + 2*ϵ*_{1} + *E*_{0} of this ground state comprises: the work *V*(*r*) − *V*(*∞*) needed to bring two charges from infinity to distance *r*, plus twice the excitation energy *ϵ*_{1} of an isolated pinned charge, on top of the dressed-vacuum energy *E*_{0}. Therefore we can estimate the interaction potential as *V*(*r*) = *E*(*r*) − *E*_{0} + *ξ* where the additive constant *ξ* does not scale with the volume (while *E*(*r*) and *E*_{0} separately do).

The presence of dynamical matter heavily impacts the strong-coupling picture (\({g}_{{\rm{m}}}^{2} \sim 0\)), as it can be extrapolated in the semiclassical limit (*t* ~ 0). Here, a particle-antiparticle pair at distance *r* with, a field-string between them, has an energy

that scales linearly with *r* (here \({g}^{2}=a{g}_{{\rm{e}}}^{2}\)). By contrast, two mesons (neighboring particle–antiparticle pairs) have a flat energy profile

Thus, for any mass *m*, there is critical distance *r*_{0} above which the string is broken, and the formation of two mesons is energetically favorable.

We observe this transition at finite *t*, as shown in Fig. 4 (bottom panel, *g*^{2} = 4). The crossover from the short-range to long-range behavior is still relatively sharp, and the distance *r*_{c} at which it occurs strongly depends on the bare mass *m*. This is in contrast to the weak-coupling regime (top panel, *g*^{2} = 1/4), where the potential profile *V*(*r*) is smoothly increasing with *r*, and its slope at short distances disagrees with the string tension ansatz *r**g*^{2}/2 + const.. Thus our simulations highlight visibly different features between confined and deconfined regimes, even with the dynamical matter.

### Finite density

One of the most important features of our numerical approach is the possibility to tackle finite charge-density regimes. In fact, by exploiting the global *U*(1) Fermion-number symmetry, implemented in our TTN algorithms, we can inject any desired charge imbalance into the system, while working under OBC. Figure 5 shows the results for charge density *ρ* = *Q*/*L*^{3} = 1/4. In the vacuum phase (\(m\gg {g}_{{\rm{e}}}^{2}/2\approx t\)), we obtain configurations as displayed in panel (a), where the charges are expelled from the bulk and stick to the boundaries to minimize the electric field energy of the outcoming fields. To quantify this effect, which can also be interpreted as a field-screening phenomenon, we introduce the surface charge density

where *A*(*l*) contains only sites sitting at lattice distance *l* from the closest boundary. The deeper we are in the vacuum phase, the faster the surface charge decays to zero away from the boundary (*l* = 1). By contrast, close to the transition, the spontaneous creation of charge-anticharge pairs determines a finite charge density of the bulk. Finally, for large negative *m*, the charge distribution is roughly uniform.

## Discussion

We have shown that TN methods can simulate LGT in three spatial dimensions, in the presence of matter and charge imbalance, ultimately exploring those regimes where other known numerical strategies struggle. We have investigated collective phenomena of lattice QED which stand at the forefront of the current research efforts, including quantum phase diagrams, confinement issues, and the string breaking mechanism at equilibrium. We envision the possibility of including more sophisticated diagnostic tools, such as the ’t Hooft operators^{70} which nicely fit TNs designs, to provide more quantitatively precise answers to the aforementioned open problems.

From a theoretical standpoint, our work corroborates the long-term perspective to employ TN methods to efficiently tackle non-perturbative phenomena of LGTs, in high dimensions and in regimes that are out of reach for other numerical techniques. As ansatz states with a refinement parameter chosen by the user, the bond dimension, TTNs are automatically equipped with a self-validation tool: convergence of each quantity with the bond dimension can be verified in polynomial time.

However, while TTNs perform well for small and intermediate system sizes, as the ones considered in this work (*L* = 2, 4, 8), the pathway to general LGTs analysis with large *L* is still a technical challenge. In particular, TTNs suffer from poor scalability for higher *L*, since they fail to explicitly capture area law for large systems, which denotes a possible bottleneck for further investigations toward the study of the thermodynamical limit of Abelian and non-Abelian high-dimensional LGTs. As a promising perspective, ref. ^{71} presents the *augmented* TTN ansatz which compensates this drawback offering better scalability. Further development in this direction will contribute to overcoming the current limitations of TTN in high-dimensional systems opening the pathway to the possibility of investigating realistic physics by starting from the TTN approach presented here.

Furthermore, we stress that our simulations have been performed on standard clusters by taking advantage only of OpenMP parallelization on single multi-core nodes. We have not yet exploited a full-scale parallelization on multi-node architectures. At a purely technical level, it is straightforward to upgrade our algorithms in this direction, in order to fully exploit the capabilities of high-performance computing. For instance, following the ideas presented in ref. ^{72}, each TTN variational sweep could be parallelized in a way to optimize its tensors separately on different computing nodes, so as to optimally scale the computational resources with the system size. On top of this, the implementation of tensor contractions on GPUs could be used to speed up the low-level computations as well^{73}.

In this work we consider the spin *s* = 1 representation, which leads to a local basis dimension of 267, as described in the Methods’ subsection “Fermionic compact representation of local gauge-invariant site”. Following the same theoretical construction for the local gauge-invariant sites, we estimate a local basis dimension of 1102 for the next representation of QED with *s* = 3/2, whereas, for the SU(2) Yang–Mills theory, by truncating after the first nontrivial irreducible representation and considering the spin-1/2 fermionic matter, one finds a local basis of 178 states for the cubic lattice. TTNs algorithms scale only polynomially with the local basis dimension but taking into account the aforementioned numbers, specific strategies for truncating also the local dimension in an optimal way (see for instance ref. ^{74}) could be required for studying higher representations of the gauge fields.

In conclusion, the aforementioned technical steps will be fundamental to tackle the problem of the continuum limit of realistic Abelian and non-Abelian LTGs and we foresee that, although very challenging, they are only some steps away along the path of TTNs developments presented here.

Alongside, from an experimental point of view, QLM formulations are among the most studied pathway towards the simulation of LGTs on quantum hardware^{61,75}. The recent developments in low-temperature physics and control techniques, for trapped ions, ultracold atoms, and Rydberg atoms in optical lattices, have led to the first experimental quantum simulations of one-dimensional LGTs^{55,56,57,58,59}. In this framework, numerical methods capable of accessing intermediate sizes, such as TNs, play a fundamental role as a cross-verification toolbox.

## Methods

### Fermionic compact representation of the local gauge-invariant site

In describing a framework for LGT, a common requirement of TN numerical simulations^{76,77,78,79}, as well as quantum simulations^{80,81,82,83,84,85,86}, is working with finite-dimensional local degrees of freedom. This is a hard requirement when investigating both LGT descending from high-energy quantum field theories^{87,88,89,90}, and condensed matter models with emergent gauge fields^{91,92}. While other pathways have been developed^{93,94,95,96}, in this work we adopted the well-known approach of truncating the gauge field space based on an energy density cutoff. In this section, we present the construction of the QED gauge-invariant configurations for the local sites that we exploit as a computational basis in our TN algorithm.

The use of the spin-1 representation implies that the gauge degrees of freedom on each link of the lattice is represented by three orthogonal eigenstates of the electric field operator

The parallel transporter, which is proportional to the raising operator in the spin language, acts on these states as

In the following, in order to obtain a representation of the gauge degrees of freedom that will be useful for constructing our TN ansatz, we employ the local mapping presented in ref. ^{47} (see also^{97,98}), generalizing it to the case with three spatial dimensions. This technique is related to the standard rishon formulation of QLM^{99,100,101} and allows us to encode Gauss’s law taking into account the anticommutation relations of the fermionic particles on the lattice.

Let us consider a generic link of the lattice (**x**, *μ*) between the two sites **x** and **x** + *μ*: the starting point is the splitting of the gauge field of this link into a pair of *rishon* modes so that each mode belongs to either one of the two sites. For the *s* = 1 case, we can set each rishon mode (or half-link) to be a three-hardcore fermionic field \({\hat{\eta }}_{{\bf{x}},\mu }\). Such lattice quantum fields satisfy \({\hat{\eta }}_{{\bf{x}},\mu }^{2}\,\ne\,0\) and \({\hat{\eta }}_{{\bf{x}},\mu }^{3}=0\). They mutually anticommute at different spatial positions, i.e., \(\left\{{\hat{\eta }}_{{\bf{x}},\mu },{\hat{\eta }}_{{\bf{x}}^{\prime} ,\mu ^{\prime} }^{(\dagger )}\right\}=0\) for \({\bf{x}}\, \ne\, {\bf{x}}^{\prime} \) or \(\mu\, \ne\, \mu ^{\prime} \), and also anticommute with the staggered matter fermionic fields \(\left\{{\hat{\eta }}_{{\bf{x}},\mu },{\hat{\psi }}_{{\bf{x}}^{\prime} }^{(\dagger )}\right\}=0\)^{17,18}. Then, we express the comparator on the link as \({\hat{U}}_{{\bf{x}},\mu }={\hat{\eta }}_{{\bf{x}},\mu }{\hat{\eta }}_{{\bf{x}}+\mu ,-\mu }^{\dagger }\). To explicitly build these three-hardcore fermions for each half-link, we consider two species of standard Dirac fermions \({\hat{a}}_{{\bf{x}},\mu }\) and \({\hat{b}}_{{\bf{x}},\mu }\) and we use the following relation:

where \({\hat{n}}_{{\bf{x}},\mu }^{a}\) and \({\hat{n}}_{{\bf{x}},\mu }^{b}\) are the occupation number operators for each species, i.e., \({\hat{n}}_{{\bf{x}},\mu }^{a}={\hat{a}}_{{\bf{x}},\mu }^{\dagger }{\hat{a}}_{{\bf{x}},\mu }\) and the same for \({\hat{n}}_{{\bf{x}},\mu }^{b}\). For each three-hardcore mode, these operators act on a three-dimensional local Hilbert space with basis \({\left|0\right.> }_{{\bf{x}},\mu }\), \({\left|1\right.\, > \, }_{{\bf{x}},\mu }={\hat{a}}_{{\bf{x}},\mu }^{\dagger }{\left|0\right.\, > \, }_{{\bf{x}},\mu }\), \({\left|2\right.\, > \, }_{{\bf{x}},\mu }={\hat{b}}_{{\bf{x}},\mu }^{\dagger }{\hat{a}}_{{\bf{x}},\mu }^{\dagger }{\left|0\right.\, > \, }_{{\bf{x}},\mu }\). In fact, due to the definition in Eq. (9), the algebra of the operators \({\hat{\eta }}_{{\bf{x}},\mu }\) never accesses the fourth state obtained as \({\hat{b}}_{{\bf{x}},\mu }^{\dagger }{\left|0\right.\, > \, }_{{\bf{x}},\mu }\). By using the same representation on the other half-link through the Dirac operators \({\hat{a}}_{{\bf{x}}+\mu ,-\mu }^{\dagger }\) and \({\hat{b}}_{{\bf{x}}+\mu ,-\mu }^{\dagger }\), we would obtain for the complete link a local space of dimension 9. However, the operator that counts the total number of fermions on the complete link as

defines asymmetry of the Hamiltonian since it commutes with the operators \({\hat{E}}_{{\bf{x}},\mu }\) and \({\hat{U}}_{{\bf{x}},\mu }\). Thus, we can select the sector with \({\hat{L}}_{{\bf{x}},\mu }=2\) (two rishons on each full link), reducing the link space to dimension 3 with the basis

where the minus sign in the first element allows the operator \({\hat{U}}_{{\bf{x}},\mu }\) to act correctly following the properties of Eq. (8). By using this representation, the electric field finally corresponds to the imbalance of Dirac fermions between the two halves of the link, so that

This construction in terms of 3-hardcore fermions allows us to define, for each lattice site, a local basis that directly incorporates Gauss’s law, by constraining in this way the dynamics to the physical states only. This is a crucial point for both numerical and quantum simulations since non-physical states determine an exponential increase in the complexity of the problem.

From the definition of the link basis states of Eq. (11), it follows that, within the sector with the link-symmetry constraint \({\hat{L}}_{{\bf{x}},\mu }=2\), the electric field operator is uniquely identified by taking only the half-link fermionic configuration, namely

In this way, the generators of the Gauss’s law of Eq. (3) are transformed into completely local operators acting on the site **x** only

Taking into account this property, it is possible to construct the gauge-invariant basis for the local site **x**, which is composed by the lattice site and the six half-links along with the directions ±*μ*_{x}, ±*μ*_{y}, ±*μ*_{z} (see Fig. 6)

where \({\left|\phi \right\rangle }_{{\bf{x}}}={({\hat{\psi }}_{{\bf{x}}}^{\dagger })}^{\phi }\left|0\right\rangle \) with *ϕ* = 0, 1 describes the presence or the absence of the matter/antimatter particles. The indices *k*_{j} run over {0, 1, 2} selecting a configuration of the 3-hardcore modes for each respective half-link. The presence of the factor \({\left(-1\right)}^{{\delta }_{{k}_{1},2}+{\delta }_{{k}_{2},2}+{\delta }_{{k}_{3},2}}\) allows us to satisfy the anticommutation relations of the fermionic representation recovering the correct signs of Eq. (11). The occupation numbers *ϕ* and *k*_{j} are not independent due to the constraint imposed by Gauss’s law

This equation, in the new language of matter fermions and rishons, reads

where the factor 6 is indeed the coordination number of the cubic lattice. Thus, the gauge-invariant configurations of the local basis are obtained by applying this constraint, effectively reducing the “dressed-site” (matter and 6 rishon modes) dimension from 2 × 3^{6} = 1458 to merely 267. We encode these states as building blocks of our computational representation for the TN algorithms. In Fig. 6, we show some examples of gauge-invariant configurations for even and odd sites.

The construction of the gauge-invariant local sites is particularly advantageous for our numerical purposes: in fact, it is now possible to express all the terms in the Hamiltonian of Eqs. (1a)–(1c) of the main text as the product of completely local operators that commute on different sites. Let us consider the kinetic term of the Hamiltonian and apply the representation of the gauge field in terms of the 3-hardcore fermionic modes

where the indices *α* and \(\alpha ^{\prime} \) select the right operators depending on the different directions in which the hopping process takes place. The operators \({M}_{{\bf{x}},\mu }^{\alpha }\) are genuinely local (i.e., they commute with operators acting elsewhere) as they are always quadratic in the fermionic operators (*ψ* and/or *η*). The same argument applies to the magnetic (plaquette) terms in the Hamiltonian

where the indices *α*, \(\alpha ^{\prime} \), *α*″, *α‴* depend on the plane of the plaquette (in this case *x–y*) and the links involved in the loop. The operators \({C}_{{\bf{x}}}^{\alpha }\) are genuinely local and act on the four sites at the corners of the plaquette. The decomposition is the same for the other plaquettes in the planes *x*–*z* and *y–z*. The present construction ensures that they can be treated as spin (or bosonic) operators^{97,98}, so we can exploit standard TN algorithms, without the need of explicitly implementing the fermionic parity at each site^{102,103,104}.

The mass term and the electric field energy in the Hamiltonian of Eqs. (1a)–(1c) of the main text are diagonal in the gauge-invariant basis with the rishon representation and so it is trivial to express them as local operators. These operators include the local chemical potential terms, which we use to pin charges in order to study confinement properties^{105,106,107}. In conclusion, all the operators we employ in the TTN algorithms (see the “Methods” subsection “Tensor Networks”) are genuinely local. In order to get an idea of the numerical complexity, we emphasize that the dimension of these matrices acting on the local gauge-invariant basis is 267 × 267.

### Tensor networks

In this section, we present the main concepts of TNs with a particular focus on the TTN ansatz that we exploit in this work^{108}. For a detailed and exhaustive description of the subject, please see the technical reviews and textbooks^{21,109,110}.

Let us consider a generic quantum system composed by *N* lattice sites, each of which described by a local Hilbert space *H*_{k} of finite dimension *d* and equipped with a local basis \({\left\{{\left|i\right.> }_{k}\right\}}_{1\le i\le d}\). The whole Hilbert space of the system will be generated by the tensor product of the local Hilbert spaces, that is, \({\mathcal{H}}={{\mathcal{H}}}_{1}\otimes {{\mathcal{H}}}_{2}\otimes \cdots {{\mathcal{H}}}_{N}\), with a resulting dimension equal to *d*^{N}. Thus, a generic pure quantum state of the system \(|\psi \rangle\) can be expressed as a linear combination of the basis elements of \({\mathcal{H}}\), i.e.,

In principle, the coefficients \({c}_{{i}_{1},...,{i}_{N}}\) are *d*^{N} complex numbers. As a consequence, this exact representation of the quantum state is completely inefficient from a computational point of view, since it scales exponentially with the system size *N*. In other words, the amount of information that we would need to store in memory for a computational representation of the generic quantum state of the system is exponentially large in the number of degrees of freedom.

However, if we are concerned with local Hamiltonians, which means that a lattice site interacts only with a finite set of neighboring sites and not with all sites of the lattice, it is possible to exploit rigorous results on the scaling of entanglement under a bipartition (*area law*)^{111,112} in order to obtain an efficient representation of the states in the low-energy sectors of such Hamiltonians, e.g., ground-states and first excited states. TN provide a natural language for this representation^{113,114} by decomposing the complete rank-*N* tensor \({c}_{{i}_{1},...,{i}_{N}}\) in Eq. (20) into a network of smaller-rank local tensors interconnected with auxiliary indices (bond indices). If we control the dimension of the bond indices with a parameter *χ*, called the bond dimension, the number of coefficients in the TN is of the order *O*(poly(*N*)poly(*χ*)), allowing an efficient representation of the information encoded in the quantum state. Furthermore, the bond dimension *χ* is a quantitative estimate of the number of quantum correlations and entanglement present in the TN. In fact, by varying *χ*, TNs interpolate between a product state (*χ* = 1) and the exact, inefficient, representation of the considered quantum state (*χ* ≈ *d*^{N}).

Matrix product states (MPS) for 1D systems^{115,116,117}, projected entangled pair state (PEPS) for 2D and 3D systems^{114,118,119}, multiscale entanglement renormalization ansatz^{120,121} and TTN, that can be defined in any dimension,^{108,122,123} are all important examples of efficient representations based on TNs.

MPS algorithms, such as the density matrix renormalization group^{124}, represent the state-of-the-art technique for the numerical simulation of many-body systems in 1D. MPS satisfy area law and are extremely powerful since they allow to compute scalar products between two wave functions and local observables in an exact and efficient way. This property does not hold true for higher-dimensional generalizations, such as PEPS, and the development of TN algorithms, for accurate and efficiently scalable computations, is at the center of current research efforts.

In particular, one of the main problems is related to the choice of the TN geometry for simulating higher-dimensional systems. PEPS intuitively reproduces the structure of the lattice with one tensor for each physical site and the bond indices directly follow the lattice grid. The resulting TN follows the area-law of entanglement but it contains loops, making the contractions for computing expectation values exponentially hard^{125}. Furthermore, the computational cost for performing the variational optimization of PEPS, as for instance in the ground state searching, scales as *O*(*χ*^{10}) as a function of the bond dimension. This severely limits the possibility of reaching high values of *χ*, especially for large system sizes (typical values are *χ* ≈10 for spin systems). For our purpose of simulating LGT in three-spatial dimensions, this represents a crucial problem since the local dimension of our model is extremely high, i.e., *d* = 267, and so, it becomes necessary to be able to handle high values of *χ* in order to reach the numerical convergence.

Alternative ansätze for simulating quantum many-body systems are the TTNs, which decompose the wave function into a network of tensors without loops, allowing efficient contraction algorithms with a polynomial scaling as a function of the system size. In Fig. 7, we show the typical TTN ansazts for 1D and 2D systems and our generalization to the 3D lattice. TTNs offer more tractable computational costs since the complete contraction and the variational optimization algorithms scale as *O*(*χ*^{4}), making it easier to reach high values of the bond dimension (up to *χ* ≈1000). The price to pay for using the loopless structure is related to the area law that TTNs may not explicitly reproduce in dimensions higher than one^{126}. Nevertheless, we use the TTN ansatz in a variational optimization, so we can improve the precision by using increasing values of *χ*, providing in this way a careful control over the convergence of our numerical results.

Ground state computation of our LGT model employs the TTN algorithm for variational ground state search, including the exploitation of Abelian symmetries and the Krylov subspace expansion^{110}. The algorithm is implemented to conserve the total charge through the definition of global *U*(1) symmetry sectors encoded in the TTN. Thus, we can easily access finite charge–density regimes, with an arbitrary imbalance between charges and anticharges.

Our TTN for the 3D lattice is composed entirely of tensors with three links (this structure is usually called *binary tree*). The construction of the TTN starts from merging the physical indices at the bottom, which represent two neighboring lattice sites along the *x*-direction, into one tensor. Then, these tensors are connected along the *y*-direction through new tensors in an upper layer. The tensors in this layer are then connected along the *z*-direction through a new layer of tensors. Thus, this procedure is iteratively repeated by properly setting the connections along with the three spatial directions in the upper layers of the tree. At the beginning of the simulation, we randomly initialize all the tensors in the network and the distribution of the global symmetry sectors. During the variational optimization stage, in order to improve the convergence, we perform the single-tensor optimization with subspace-expansion technique, i.e., allowing a dynamical increase of the local bond dimension and adapting the symmetry sectors^{110}. This scheme has a global computational cost of the order *O*(*χ*^{4}). The single tensor optimization is implemented in three steps: (i) the effective Hamiltonian *H*_{eff} for the tensor is obtained by contracting the complete Hamiltonian of the system with all the remaining tensors of the tree; (ii) the local eigenvalue problem for *H*_{eff} is solved by using the Arnoldi method of the ARPACK library; (iii) the tensor is updated by the eigenvector of *H*_{eff} corresponding to the lowest eigenvalue. This procedure is iterated by sweeping through the TTN from the lowest to the highest layers, gradually reducing the energy expectation value. After completing the whole sweep, the procedure is iterated again and again, until the desired convergence in the energy is reached. The precision of the Arnoldi algorithm is increased in each sweep, for gaining more accuracy in solving the local eigenvalue problems as we approach the final convergence.

TTN computations presented in this work are extremely challenging due to the complexity of LGTs in the three-dimensional scenario. They were performed on different HPC-clusters (CloudVeneto, CINECA, BwUniCluster, and ATOS Bull): a single simulation for the maximum size that we reached, an 8 × 8 × 8 lattice, can last up to five weeks until final convergence, depending on the different regimes of the model and the control parameters of the algorithms.

### Numerical convergence

With our numerical simulations, we characterize the properties of the ground state of the system as a function of the parameters in the Hamiltonian of Eqs. (1a)–(1c) of the main text. We fix the energy scale by setting the hopping coefficient *t* = 1 and we access several regimes of the mass *m*, the electric *g*_{e} and the magnetic coupling *g*_{m}. We consider simple cubic lattices *L* × *L* × *L* with the linear size *L* being a binary power; in particular, we simulate the case with *L* = 2, 4, 8, that is, up to 512 lattice sites.

As explained in the “Methods” subsection “Fermionic compact representation of local gauge-invariant site”, in order to obtain the right representation of the electric field operators, we have to enforce the extra link symmetry constraint \({\hat{L}}_{{\bf{x}},\mu }=2\) at every pair of neighboring sites. For this reason, we include in the Hamiltonian additional terms that energetically penalize all the states with a number of hardcore fermions per link different from two, namely

where *ν* > 0 is the penalty coefficient and \({\hat{\delta }}_{2,{\hat{L}}_{{\bf{x}},\mu }}\) are the projectors on the states that satisfy the extra link constraint. In this way, the penalty terms vanish when the link symmetry is satisfied and raise the energy of the states violating the constraint. In principle, the link symmetry is rigorously satisfied for *ν* → *∞*. At the numerical level, this limit translates into choosing *ν* much larger than the other simulation parameters of the Hamiltonian, i.e., \(\nu \gg \max \left\{| t| ,| m| ,| {g}_{{\rm{e}}}| ,| {g}_{{\rm{m}}}| \right\}\). However, setting *ν* too large in the first optimization steps could lead to local minima or non-physical states, since the variational algorithm would focus only on the penalty terms more than the physical ones. In order to avoid this problem and reach the convergence, we adopt a driven optimization, by varying the penalty coefficient *ν* in three steps: (i) starting from a very small value of *ν* and from a random state of the TTN, that in general does not respect the extra link symmetry, we drive the penalty term with linear growth of *ν* during the first optimization sweeps. In this stage, the optimization will focus mainly on the physical quantities, until we notice a slight rise in the energy: this effect signals that the global optimization procedure of the TTN becomes significantly sensitive to the penalty terms. (ii) Thus, we impose a quadratic growth of *ν* so that, in the immediately following sweeps, the penalty is increased at a slower rate with respect to the linear regime. (iii) After reaching the maximum desired value of *ν*, which is an input parameter of the simulation, we keep it fixed, performing the last sweeps in order to ensure the convergence of the energy. This driven optimization strategy is summarized in Fig. 8a where we show the three different stages of the penalty coefficient *ν* and typical behavior of the energy difference *δ**e*, computed with respect to the lowest final energy that we reach, as a function of the iterations.

We can also quantify the global error with respect to the link symmetry during the driven optimization sweeps, by defining

i.e., the sum of the deviations from the exact link constraint \({\hat{L}}_{{\bf{x}},\mu }=2\), computed over all the links of the lattice on the ground state. The typical behavior of this quantity is shown in Fig. 8b: at the end of the optimization procedure, the global error results of the order of 10^{−6}. We also check the convergence of our TTN algorithms as a function of the bond dimension *χ*, by using *χ* = 450 at most to ensure the stability of our findings. Depending on the different system sizes and regimes of physical parameters, we estimate the relative error of the energy in the range \(\left[1{0}^{-2},1{0}^{-4}\right]\). A typical scaling of the energy density as a function of the inverse of the bond dimension 1/*χ* is shown in Fig. 8c.

### Critical points: scaling analysis

In this section, we show the finite-size scaling analysis for detecting the phase transition separating the charge-crystal phase and the vacuum phase and the related location of the critical points.

At *t* = 0 and neglecting the magnetic interactions, i.e., for \({g}_{{\rm{m}}}^{2}=0\), the Hamiltonian of Eq. (1a)–(1c) results diagonal in the local basis described in the “Methods” subsection “Fermionic compact representation of local gauge-invariant site” and it is trivial to prove that the system undergoes a first-order phase transition between the bare vacuum, with energy \({E}_{{\rm{v}}}=-m\frac{{L}^{3}}{2}\) and the charge-crystal phase, with energy \({E}_{{\rm{ch}}}=(m+\frac{{g}_{e}^{2}}{2})\frac{{L}^{3}}{2}\). The ground-state exhibits a level-crossing at the critical value \({m}_{{\rm{c}}}^{(0)}=-\frac{{g}_{{\rm{e}}}^{2}}{4}=-\frac{1}{2}\) that is obtained at *E*_{v} = *E*_{ch}. This behavior is clearly seen in Fig. 9a, showing a discontinuous transition between the two configurations.

In order to understand the behavior of the system for finite *t* = 1 and \({g}_{{\rm{m}}}^{2}=0\), we observe that the density, plotted in Fig. 2a of the main text, changes continuously as a function of the mass parameter and we might have a second-order phase transition. Finite-size scaling theory^{127} implies that the behavior of the system close to a critical point, i.e., for *m* ≈ *m*_{c}, can be described in terms of a universal function *λ*(*x*) such that for our observable

where *β* and *ν* are critical exponents. In particular, this relation implies that for *m* ≈ *m*_{c}, the value of \(\rho {L}^{\frac{\beta }{\nu }}\) is independent of the size of the system. We use this property to get an estimate of the values of *m*_{c}, *β*, *ν*. In particular, we consider a grid of values for this parameter, and for each set of values, we fit our points \(\rho {L}^{\frac{\beta }{\nu }}\) with an high-degree polynomial \(f\left({L}^{\frac{1}{\nu }}(m-{m}_{{\rm{c}}})\right)\). We compute the residual sum of squares (RSS) and we select the set of values which minimize this quantity, producing the best data collapse. We get for the critical point *m*_{c} ≈ −0.39 and for the critical exponents *β* ≈0.16 and *ν* ≈1.22. In Fig. 9(b) we show the collapse of our numerical results onto the same universal function *λ*(*x*) and in Fig. 9(c) a contour plot of the square root of the RSS in the (*ν*, *β*) plane for the best-fitting values.

By extending the previous considerations and the finite-size scaling analysis to the case with magnetic interactions with \({g}_{{\rm{m}}}^{2}=8/{g}_{{\rm{e}}}^{2}=4\), we check again the presence of a critical point and the values of critical exponents through the formula of Eq. (23). We obtain a universal scaling function for *m*_{c} ≈ 0.22 and the same critical exponents *β* ≈0.16, *ν* ≈1.22, as reported in the inset of Fig. 9b. Thus, while the transition and its universality remain unchanged in the presence of the magnetic coupling, the critical point is shifted toward positive values of the mass parameter, signaling that the magnetic interactions determine a visible enhancement of the production of charges and anti charges out of the vacuum.

Although a more precise determination of the numerical values of the critical exponents would require additional extensive analysis that results beyond the scope of this paper, our findings strongly indicate the presence of a phase transition at finite *m* for the three-dimensional lattice model of QED (compare with other previously investigated transitions in lattice QED, e.g., ref. ^{128}).

## Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request. Source data of the figures are attached to the manuscript. Source data are provided with this paper.

## Code availability

The authors are available for discussing details of the implementation of the computer codes developed and used in this study upon reasonable request.

## References

- 1.
Peskin, M. E. & Schroeder, D. V.

*An Introduction to Quantum Field Theory*. (Westview, Boulder, CO, 1995). - 2.
Cheng, T. & Li, L.

*Gauge Theory of Elementary Particle Physics*(Oxford University Press, 2006). - 3.
Schwartz, M. D.

*Quantum Field Theory and the Standard Model*(Cambridge University Press, 2014). - 4.
Fradkin, E.

*Field Theories of Condensed Matter Physics*(Cambridge University Press, 2013). - 5.
Zhou, Y., Kanoda, K. & Ng, T.-K. Quantum spin liquid states.

*Rev. Mod. Phys.***89**, 025003 (2017). - 6.
Brambilla, N. et al. Qcd and strongly coupled gauge theories: challenges and perspectives.

*Eur. Phys. J. C***74**, 2981 (2014). - 7.
Fiebig, H. R. & Woloshyn, R. M. Monopoles and chiral-symmetry breaking in three-dimensional lattice qed.

*Phys. Rev. D***42**, 3520–3523 (1990). - 8.
Herbut, I. F. & Seradjeh, B. H. Permanent confinement in the compact qed

_{3}with fermionic matter.*Phys. Rev. Lett.***91**, 171601 (2003). - 9.
Nogueira, F. S. & Kleinert, H. Compact quantum electrodynamics in 2 + 1 dimensions and spinon deconfinement: a renormalization group analysis.

*Phys. Rev. B***77**, 045107 (2008). - 10.
Xu, X. Y. et al. Monte carlo study of lattice compact quantum electrodynamics with fermionic matter: the parent state of quantum phases.

*Phys. Rev. X***9**, 021022 (2019). - 11.
Wilson, K. G. Confinement of quarks.

*Phys. Rev. D***10**, 2445–2459 (1974). - 12.
Alkofer, R. & Greensite, J. Quark confinement: the hard problem of hadron physics.

*J. Phys. G***34**, S3–S21 (2007). - 13.
Creutz, M. Lattice gauge theories and monte carlo algorithms.

*Nucl. Phys. B***10**, 1–22 (1989). - 14.
Bazavov, A. et al. Nonperturbative qcd simulations with 2 + 1 flavors of improved staggered quarks.

*Rev. Mod. Phys.***82**, 1349–1417 (2010). - 15.
Tanabashi, M. et al. Review of particle physics.

*Phys. Rev. D***98**, 030001 (2018). - 16.
Kogut, J. & Susskind, L. Hamiltonian formulation of wilson’s lattice gauge theories.

*Phys. Rev. D***11**, 395–408 (1975). - 17.
Kogut, J. B. An introduction to lattice gauge theory and spin systems.

*Rev. Mod. Phys.***51**, 659–713 (1979). - 18.
Susskind, L. Lattice fermions.

*Phys. Rev. D***16**, 3031–3039 (1977). - 19.
Troyer, M. & Wiese, U.-J. Computational complexity and fundamental limitations to fermionic quantum monte carlo simulations.

*Phys. Rev. Lett.***94**, 170201 (2005). - 20.
Verstraete, F., Murg, V. & Cirac, J. I. Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems.

*Adv. Phys.***57**, 143–224 (2008). - 21.
Orús, R. A practical introduction to tensor networks: matrix product states and projected entangled pair states.

*Ann. Phys.***349**, 117–158 (2014). - 22.
Haegeman, J., Lubich, C., Oseledets, I., Vandereycken, B. & Verstraete, F. Unifying time evolution and optimization with matrix product states.

*Phys. Rev. B***94**, 165116 (2016). - 23.
Byrnes, T. M. R., Sriganesh, P., Bursill, R. J. & Hamer, C. J. Density matrix renormalization group approach to the massive schwinger model.

*Phys. Rev. D***66**, 013002 (2002). - 24.
Bañuls, M. C., Cichy, K., Cirac, J. I. & Jansen, K. The mass spectrum of the schwinger model with matrix product states.

*J. High Energy Phys.***2013**, 158 (2013). - 25.
Rico, E., Pichler, T., Dalmonte, M., Zoller, P. & Montangero, S. Tensor networks for lattice gauge theories and atomic quantum simulation.

*Phys. Rev. Lett.***112**, 201601 (2014). - 26.
Buyens, B., Haegeman, J., Van Acoleyen, K., Verschelde, H. & Verstraete, F. Matrix product states for gauge field theories.

*Phys. Rev. Lett.***113**, 091601 (2014). - 27.
Bañuls, M. C., Cichy, K., Cirac, J. I., Jansen, K. & Saito, H. Thermal evolution of the schwinger model with matrix product operators.

*Phys. Rev. D***92**, 034519 (2015). - 28.
Kühn, S., Zohar, E., Cirac, J. I. & Bañuls, M. C. Non-abelian string breaking phenomena with matrix product states.

*J. High Energy Phys.***2015**, 1–26 (2015). - 29.
Buyens, B., Haegeman, J., Verschelde, H., Verstraete, F. & Van Acoleyen, K. Confinement and string breaking for qed

_{2}in the hamiltonian picture.*Phys. Rev. X***6**, 041040 (2016). - 30.
Buyens, B., Verstraete, F. & Van Acoleyen, K. Hamiltonian simulation of the Schwinger model at finite temperature.

*Phys. Rev. D***94**, 085018 (2016). - 31.
Buyens, B., Montangero, S., Haegeman, J., Verstraete, F. & Van Acoleyen, K. Finite-representation approximation of lattice gauge theories at the continuum limit with tensor networks.

*Phys. Rev. D***95**, 094509 (2017). - 32.
Bañuls, M. C., Cichy, K., Cirac, J. I., Jansen, K. & Kühn, S. Efficient basis formulation for (1 + 1)-dimensional su(2) lattice gauge theory: spectral calculations with matrix product states.

*Phys. Rev. X***7**, 041046 (2017). - 33.
Bañuls, M. C., Cichy, K., Cirac, J. I., Jansen, K. & Kühn, S. Density induced phase transitions in the schwinger model: a study with matrix product states.

*Phys. Rev. Lett.***118**, 071601 (2017). - 34.
Buyens, B., Haegeman, J., Hebenstreit, F., Verstraete, F. & Van Acoleyen, K. Real-time simulation of the schwinger effect with matrix product states.

*Phys. Rev. D***96**, 114501 (2017). - 35.
Kull, I., Molnar, A., Zohar, E. & Cirac, J. I. Classification of matrix product states with a local (gauge) symmetry.

*Ann. Phys.***386**, 199–241 (2017). - 36.
Sala, P. et al. Gaussian states for the variational study of (1+1)-dimensional lattice gauge models. in

*The 36th Annual International Symposium on Lattice Field Theory, 22–28 July*, 230 (2018). - 37.
Ercolessi, E., Facchi, P., Magnifico, G., Pascazio, S. & Pepe, F. V. Phase transitions in

*Z*_{n}gauge models: towards quantum simulations of the schwinger-weyl qed.*Phys. Rev. D***98**, 074503 (2018). - 38.
Magnifico, G. et al. Symmetry-protected topological phases in lattice gauge theories: topological qed

_{2}.*Phys. Rev. D***99**, 014503 (2019). - 39.
Huang, Y.-P., Banerjee, D. & Heyl, M. Dynamical quantum phase transitions in u(1) quantum link models.

*Phys. Rev. Lett.***122**, 250401 (2019). - 40.
Magnifico, G. et al. Real time dynamics and confinement in the \({{\mathbb{Z}}}_{n}\) Schwinger-Weyl lattice model for 1+1 QED.

*Quantum***4**, 281 (2020). - 41.
Funcke, L., Jansen, K. & Kühn, S. Topological vacuum structure of the schwinger model with matrix product states.

*Phys. Rev. D***101**, 054507 (2020). - 42.
Tagliacozzo, L. & Vidal, G. Entanglement renormalization and gauge symmetry.

*Phys. Rev. B***83**, 115127 (2011). - 43.
Tagliacozzo, L., Celi, A., Zamora, A. & Lewenstein, M. Optical abelian lattice gauge theories.

*Ann. Phys.***330**, 160–191 (2013). - 44.
Tagliacozzo, L., Celi, A. & Lewenstein, M. Tensor networks for lattice gauge theories with continuous groups.

*Phys. Rev. X***4**, 041024 (2014). - 45.
Zohar, E., Burrello, M., Wahl, T. B. & Cirac, J. I. Fermionic projected entangled pair states and local u(1) gauge theories.

*Ann. Phys.***363**, 385–439 (2015). - 46.
Zohar, E. & Cirac, J. I. Combining tensor networks with monte carlo methods for lattice gauge theories.

*Phys. Rev. D***97**, 034510 (2018). - 47.
Felser, T., Silvi, P., Collura, M. & Montangero, S. Two-Dimensional Quantum-Link Lattice Quantum Electrodynamics at Finite Density.

*Phys. Rev. X***10**, 041040 (2020). - 48.
Emonts, P., Bañuls, M. C., Cirac, I. & Zohar, E. Variational monte carlo simulation with tensor networks of a pure \({{\mathbb{z}}}_{3}\) gauge theory in (2 + 1)D.

*Phys. Rev. D***102**, 074501 (2020). - 49.
Haegeman, J., Van Acoleyen, K., Schuch, N., Cirac, J. I. & Verstraete, F. Gauging quantum states: from global to local symmetries in many-body systems.

*Phys. Rev. X***5**, 011024 (2015). - 50.
Meurice, Y., Sakai, R. & Unmuth-Yockey, J. Tensor field theory with applications to quantum computing. Preprint at

*arXiv*https://arxiv.org/abs/2010.06539 (2020). - 51.
Horn, D. Finite matrix models with continuous local gauge invariance.

*Phys. Lett. B***100**, 149–151 (1981). - 52.
Chandrasekharan, S. & Wiese, U.-J. Quantum link models: a discrete approach to gauge theories.

*Nucl. Phy. B***492**, 455–471 (1997). - 53.
Wiese, U.-J. Ultracold quantum gases and lattice systems: quantum simulation of lattice gauge theories.

*Ann. Phys.***525**, 777–796 (2013). - 54.
Zohar, E., Cirac, J. I. & Reznik, B. Quantum simulations of gauge theories with ultracold atoms: Local gauge invariance from angular-momentum conservation.

*Phys. Rev. A***88**, 023617 (2013). - 55.
Martinez, E. A. et al. Real-time dynamics of lattice gauge theories with a few-qubit quantum computer.

*Nature***534**, 516–519 (2016). - 56.
Klco, N. et al. Quantum-classical computation of schwinger model dynamics using quantum computers.

*Phys. Rev. A***98**, 032331 (2018). - 57.
Mil, A. et al. A scalable realization of local u(1) gauge invariance in cold atomic mixtures.

*Science***367**, 1128–1130 (2020). - 58.
Yang, B. et al. Observation of gauge invariance in a 71-site Bose-Hubbard quantum simulator. Preprint at

*arXiv*https://arxiv.org/abs/2003.08945 (2020). - 59.
Schweizer, C. et al. A scalable realization of local u(1) gauge invariance in cold atomic mixtures.

*Science***367**, 1128–1130 (2020). - 60.
Kasper, V., Hebenstreit, F., Jendrzejewski, F., Oberthaler, M. K. & Berges, J. Implementing quantum electrodynamics with ultracold atomic systems.

*New J. Phys.***19**, 023030 (2017). - 61.
Mathis, S. V., Mazzola, G. & Tavernelli, I. Toward scalable simulations of lattice gauge theories on quantum computers.

*Phys. Rev. D***102**, 094501 (2020). - 62.
Jordan, S. P., Lee, K. S. M. & Preskill, J. Quantum algorithms for quantum field theories.

*Science***336**, 1130–1133 (2012). - 63.
Schwinger, J. On gauge invariance and vacuum polarization.

*Phys. Rev.***82**, 664–679 (1951). - 64.
Polyakov, A. Compact gauge fields and the infrared catastrophe.

*Phys. Lett. B***59**, 82–84 (1975). - 65.
Banks, T., Myerson, R. & Kogut, J. Phase transitions in abelian lattice gauge theories.

*Nucl. Phys. B***129**, 493–510 (1977). - 66.
Drell, S. D., Quinn, H. R., Svetitsky, B. & Weinstein, M. Quantum electrodynamics on a lattice: a Hamiltonian variational approach to the physics of the weak-coupling region.

*Phys. Rev. D***19**, 619–638 (1979). - 67.
Jersák, J., Neuhaus, T. & Zerwas, P. U(1) lattice gauge theory near the phase transition.

*Phys. Lett. B***133**, 103–107 (1983). - 68.
Guth, A. H. Existence proof of a nonconfining phase in four-dimensional u(1) lattice gauge theory.

*Phys. Rev. D***21**, 2291–2307 (1980). - 69.
Kondo, K.-I. Existence of a confinement phase in quantum electrodynamics.

*Phys. Rev. D***58**, 085013 (1998). - 70.
Hooft, G. On the phase transition towards permanent quark confinement.

*Nucl. Phys. B***138**, 1–25 (1978). - 71.
Felser, T., Notarnicola, S. & Montangero, S. Efficient Tensor Network Ansatz for High-Dimensional Quantum Many-Body Problems.

*Phys. Rev. Lett*.**126**, 170603 (2021). - 72.
Stoudenmire, E. M. & White, S. R. Real-space parallel density matrix renormalization group.

*Phys. Rev. B***87**, 155137 (2013). - 73.
Efthymiou, S., Hidary, J. & Leichenauer, S. Tensornetwork for machine learning. Preprint at

*arXiv*https://arxiv.org/abs/1906.06329 (2019). - 74.
Stolpp, J. et al. Comparative study of state-of-the-art matrix-product-state methods for lattice models with large local hilbert spaces. Preprint at

*arXiv*https://arxiv.org/abs/2011.07412 (2020). - 75.
Atas, Y. et al. Su(2) hadrons on a quantum computer. Preprint at

*arXiv*https://arxiv.org/abs/2102.08920 (2021). - 76.
Silvi, P., Rico, E., Calarco, T. & Montangero, S. Lattice gauge tensor networks.

*New J. Phys.***16**, 103015 (2014). - 77.
Silvi, P., Sauer, Y., Tschirsich, F. & Montangero, S. Tensor network simulation of an su(3) lattice gauge theory in 1d.

*Phys. Rev. D***100**, 074512 (2019). - 78.
Pichler, T., Dalmonte, M., Rico, E., Zoller, P. & Montangero, S. Real-time dynamics in u(1) lattice gauge theories with tensor networks.

*Phys. Rev. X***6**, 011023 (2016). - 79.
Silvi, P., Rico, E., Dalmonte, M., Tschirsich, F. & Montangero, S. Finite-density phase diagram of a (1 + 1) −

*d*non-abelian lattice gauge theory with tensor networks.*Quantum***1**, 9 (2017). - 80.
Feynman, R. P. Simulating physics with computers.

*Int. J. Theor. Phys.***21**, 467–488 (1982). - 81.
Kokail, C. et al. Self-verifying variational quantum simulation of lattice models.

*Nature***569**, 355–360 (2019). - 82.
Bañuls, M. C. et al. Simulating lattice gauge theories within quantum technologies.

*Eur. Phys. J. D***74**, 165 (2020). - 83.
Paulson, D. et al. Towards simulating 2d effects in lattice gauge theories on a quantum computer. Preprint at

*arXiv*https://arxiv.org/abs/2008.09252 (2020). - 84.
Zohar, E. & Burrello, M. Formulation of lattice gauge theories for quantum simulations.

*Phys. Rev. D***91**, 054506 (2015). - 85.
Zohar, E., Cirac, J. I. & Reznik, B. Quantum simulations of lattice gauge theories using ultracold atoms in optical lattices.

*Rep. Progr. Phys.***79**, 014401 (2015). - 86.
Dalmonte, M. & Montangero, S. Lattice gauge theory simulations in the quantum information era.

*Contemp. Phys.***57**, 388–412 (2016). - 87.
Yang, C. N. & Mills, R. L. Conservation of isotopic spin and isotopic gauge invariance.

*Phys. Rev.***96**, 191–195 (1954). - 88.
Goldstone, J., Salam, A. & Weinberg, S. Broken symmetries.

*Phys. Rev.***127**, 965–970 (1962). - 89.
Englert, F. & Brout, R. Broken symmetry and the mass of gauge vector mesons.

*Phys. Rev. Lett.***13**, 321–323 (1964). - 90.
Higgs, P. W. Broken symmetries and the masses of gauge bosons.

*Phys. Rev. Lett.***13**, 508–509 (1964). - 91.
Baskaran, G. & Anderson, P. W. Gauge theory of high-temperature superconductors and strongly correlated fermi systems.

*Phys. Rev. B***37**, 580–583 (1988). - 92.
Kleinert, H.

*Gauge Fields in Condensed Matter*(World Scientific, 1989). - 93.
Haase, J. F. et al. A resource efficient approach for quantum and classical simulations of gauge theories in particle physics.

*Quantum***5**, 393 (2021). - 94.
Kaplan, D. B. & Stryker, J. R. Gauss’s law, duality, and the hamiltonian formulation of u(1) lattice gauge theory.

*Phys. Rev. D***102**, 094515 (2020). - 95.
Bender, J. & Zohar, E. Gauge redundancy-free formulation of compact qed with dynamical matter for quantum and classical computations.

*Phys. Rev. D***102**, 114517 (2020). - 96.
Bender, J., Emonts, P., Zohar, E. & Cirac, J. I. Real-time dynamics in 2 + 1

*d*compact qed using complex periodic gaussian states.*Phys. Rev. Res.***2**, 043145 (2020). - 97.
Zohar, E. & Cirac, J. I. Eliminating fermionic matter fields in lattice gauge theories.

*Phys. Rev. B***98**, 075119 (2018). - 98.
Zohar, E. & Cirac, J. I. Removing staggered fermionic matter in

*u*(*n*) and*s**u*(*n*) lattice gauge theories.*Phys. Rev. D***99**, 114511 (2019). - 99.
Brower, R., Chandrasekharan, S. & Wiese, U.-J. Qcd as a quantum link model.

*Phys. Rev. D***60**, 094502 (1999). - 100.
Brower, R., Chandrasekharan, S., Riederer, S. & Wiese, U.-J. D-theory: field quantization by dimensional reduction of discrete variables.

*Nucl. Phys. B***693**, 149–175 (2004). - 101.
Banerjee, D. et al. Atomic quantum simulation of

**U**(*n*) and SU(*n*) non-abelian lattice gauge theories.*Phys. Rev. Lett.***110**, 125303 (2013). - 102.
Kraus, C. V., Schuch, N., Verstraete, F. & Cirac, J. I. Fermionic projected entangled pair states.

*Phys. Rev. A***81**, 052338 (2010). - 103.
Corboz, P. & Vidal, G. Fermionic multiscale entanglement renormalization ansatz.

*Phys. Rev. B***80**, 165129 (2009). - 104.
Barthel, T., Pineda, C. & Eisert, J. Contraction of fermionic operator circuits and the simulation of strongly correlated fermions.

*Phys. Rev. A***80**, 042333 (2009). - 105.
Polyakov, A. Quark confinement and topology of gauge theories.

*Nucl. Phys. B***120**, 429–458 (1977). - 106.
Greensite, J.

*An Introduction to the Confinement Problem*. Vol. 821 (2011). - 107.
Bender, J., Emonts, P., Zohar, E. & Cirac, J. I. Real-time dynamics in 2 + 1

*d*compact qed using complex periodic gaussian states.*Phys. Rev. Res.***2**, 043145 (2020). - 108.
Gerster, M. et al. Unconstrained tree tensor network: an adaptive gauge picture for enhanced performance.

*Phys. Rev. B***90**, 125154 (2014). - 109.
Montangero, S.

*Introduction to Tensor Network Methods*(Springer International Publishing, 2018). - 110.
Silvi, P. et al. The tensor networks anthology: simulation techniques for many-body quantum lattice systems.

*SciPost Phys. Lect. Notes***8**https://scipost.org/SciPostPhysLectNotes.8 (2019). - 111.
Amico, L., Fazio, R., Osterloh, A. & Vedral, V. Entanglement in many-body systems.

*Rev. Mod. Phys.***80**, 517–576 (2008). - 112.
Eisert, J., Cramer, M. & Plenio, M. B. Colloquium: area laws for the entanglement entropy.

*Rev. Mod. Phys.***82**, 277–306 (2010). - 113.
Verstraete, F. & Cirac, J. I. Matrix product states represent ground states faithfully.

*Phys. Rev. B***73**, 094423 (2006). - 114.
Verstraete, F., Wolf, M. M., Perez-Garcia, D. & Cirac, J. I. Criticality, the area law, and the computational power of projected entangled pair states.

*Phys. Rev. Lett.***96**, 220601 (2006). - 115.
Fannes, M., Nachtergaele, B. & Werner, R. F. Finitely correlated states on quantum spin chains.

*Commun. Math. Phys.***144**, 443–490 (1992). - 116.
Klumper, A., Schadschneider, A. & Zittartz, J. Equivalence and solution of anisotropic spin-1 models and generalized t-J fermion models in one dimension.

*J. Phys. A Math. Gen.***24**, L955–L959 (1991). - 117.
Klümper, A., Schadschneider, A. & Zittartz, J. Matrix product ground states for one-dimensional spin-1 quantum antiferromagnets.

*Europhys. Lett.***24**, 293 (1993). - 118.
Verstraete, F. & Cirac, J. I. Renormalization algorithms for quantum-many body systems in two and higher dimensions. Preprint at

*arXiv*https://arxiv.org/abs/cond-mat/0407066 (2004). - 119.
Tepaske, M. & Luitz, D. J. Three-dimensional isometric tensor networks. Preprint at

*arXiv*https://arxiv.org/abs/2005.13592 (2020). - 120.
Vidal, G. Entanglement renormalization.

*Phys. Rev. Lett.***99**, 220405 (2007). - 121.
Evenbly, G. & Vidal, G. Entanglement renormalization in two spatial dimensions.

*Phys. Rev. Lett.***102**, 180406 (2009). - 122.
Shi, Y.-Y., Duan, L.-M. & Vidal, G. Classical simulation of quantum many-body systems with a tree tensor network.

*Phys. Rev. A***74**, 022320 (2006). - 123.
Silvi, P. et al. Homogeneous binary trees as ground states of quantum critical hamiltonians.

*Phys. Rev. A***81**, 062335 (2010). - 124.
Schollwöck, U. The density-matrix renormalization group in the age of matrix product states.

*Ann. Phys.***326**, 96–192 (2011). - 125.
Haferkamp, J., Hangleiter, D., Eisert, J. & Gluza, M. Contracting projected entangled pair states is average-case hard.

*Phys. Rev. Res.***2**, 013010 (2020). - 126.
Ferris, A. J. Area law and real-space renormalization.

*Phys. Rev. B***87**, 125139 (2013). - 127.
Cardy, J.

*Scaling and Renormalization in Statistical Physics*. Cambridge Lecture Notes in Physics (Cambridge University Press, 1996). - 128.
Göckeler, M. et al. Qed—a lattice investigation of the chiral phase transition and the nature of the continuum limit.

*Nucl. Phys. B***334**, 527–558 (1990).

## Acknowledgements

Authors kindly acknowledge support from the Italian PRIN2017 and Fondazione CARIPARO, the Horizon 2020 research and innovation program under grant agreement No. 817482 (Quantum Flagship—PASQuanS), the QuantERA projects QTFLAG and QuantHEP, the DFG project TWITTER, the INFN project QUANTUM, and the Austrian Research Promotion Agency (FFG) via QFTE project AutomatiQ. We acknowledge computational resources by the Cloud Veneto, CINECA, the BwUniCluster, and ATOS Bull.

## Author information

### Affiliations

### Contributions

G.M. and T.F. implemented the TTN algorithms for the three-dimensional model, basing on a code previously developed by T.F. for 2D LGT; G.M. performed the numerical simulations and analyzed the results; P.S. provided the basis for the theoretical strategy of mapping the lattice gauge theory into the proper computational framework. The interpretation of the results was mainly done by G.M., P.S., and S.M.; G.M. and P.S. wrote the paper; S.M. conceived, supervised, and directed the project. All authors discussed the results, contributed to refining the paper, and approved it.

### Corresponding author

## Ethics declarations

### Competing interests

The authors declare no competing interests.

## Additional information

**Peer review information** *Nature Communications* thanks the anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supplementary information

## Source data

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Magnifico, G., Felser, T., Silvi, P. *et al.* Lattice quantum electrodynamics in (3+1)-dimensions at finite density with tensor networks.
*Nat Commun* **12, **3600 (2021). https://doi.org/10.1038/s41467-021-23646-3

Received:

Accepted:

Published:

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.