Lattice quantum electrodynamics in (3+1)-dimensions at finite density with tensor networks

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.

Ranging from high-energy particle physics (Standard Model) [1][2][3] to low-temperature condensed matter physics (spin liquids, quantum Hall, 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 (QCD), a founding pillar of the Standard Model which have been studied for almost 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], struggles in highly interesting regimes, where matter fermions and excess of charge are concerned, due to the infamous sign problem [17]. 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 [18,19]. As tailored many-body quantum state ansätze, TNs are an efficient approximate entanglementbased 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 [20]. TN methods have proven remarkable success in simulating LGTs in (1+1) dimensions [21][22][23][24][25][26][27][28], and very recently they have shown potential in (2+1) dimensions [29][30][31][32][33]. 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.
Here, we bridge this gap by numerically simulating the Hamiltonian formulation of (3+1) Quantum Elec-trodynamics (QED) at zero temperature, via TN ansatz states. We show that, by using the quantum link formalism (QLM) of LGTs [34,35] and an unconstrained Tree Tensor Network (TTN), we can access multiple equlibrium 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 non-trivial 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 [28]. Finally, we stress that the truncation of the gauge field is a common step in quantum simulations and computations [36][37][38][39][40][41][42], making the presented numerical approach a landmark benchmarking and cross-verification tool for current and future experiments.
Hereafter, 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 [32]. In the limit of zero magnetic coupling, 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 fieldscreening 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, ( 1) i+j+k = +1 Figure 1: Scheme of the three-dimensional LGT with three electric field levels (spin-1 compact representation). Fermionic degrees of freedom are represented by staggered fermions on sites with different parity: on the even (odd) sites, a full red (blue) circle corresponds to a particle (antiparticle) with positive (negative) charge. As an illustrative example, it is shown a gauge-invariant configuration of matter and gauge fields with one particle and one antiparticle in the sector of zero total charge.
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.

I. 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]:Ĥ labelling the sites of the lattice and Here we adopted the Kogut-Susskind formulation [16], representing fermionic degrees of freedom with a staggered spinless fermion field {ψ x ,ψ † x } = δ x,x 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 negative charge −q, as shown in Fig. 1. ChargeQ conservation is thus expressed as global fermion numberN conservation, The links of the 3D lattice are uniquely identified by the couple of parameters (x, µ) where x is any site, µ is one of the three positive lattice unit vectors µ x ≡ (1, 0, 0), µ y ≡ (0, 1, 0), µ z ≡ (0, 0, 1). The gauge fields are defined on lattice links through the pair of operatorsÊ x,µ (electric field) andÛ x,µ (unitary comparator) that satisfy the commutation relation (2) For comfort of notation, we can extend the definition to negative lattice unit vectors viaÊ x+µ,−µ = −Ê x,µ and U x+µ,−µ =Û † x,µ . The Hamiltonian of Eq. (1) consists of four terms: the parallel transporter (1a) describes 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 Eq. (1) are not independent: They can be expressed as t = 1/a, m = m 0 , g 2 e = g 2 /a, g 2 m = 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 e g m = 2 √ 2t. 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 generatorŝ 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Ĝ x commutes with the HamiltonianĤ and the gauge invariant Hilbert space consists of physical many-body quantum states |Φ sat-isfyingĜ x |Φ = 0 at every site x.
As stressed in the standard Wilson's formulation of lattice QED [11], faithful representations of the (Ê,Û ) 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 SU (2) spin algebra, i.e. E x,µ ≡Ŝ z x,µ andÛ x,µ ≡Ŝ + x,µ /s for a spin-s representation. This substitution keeps the electric field operator hermitian and preserves Eq. (2), butÛ 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 2 e , which in turn requires larger spin s to accurately represent weaker coupling regimes, still potentially accessible via TNs [24].

II. TRANSITION AT ZERO CHARGE
We focus on the zero charge sector, i.e. x ψ † x ψ x = L 3 2 , and Periodic Boundary Conditions (PBC). As shown in Fig. 2 (upper panel), for g 2 m = 0 the system undergoes a transition between two regimes, analogously to the (1+1)D and (2+1)D cases [22,25,32]: 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 ρ = x ψ x is the matter occupation operator and the many-body ground state |GS has been computed by TTN algorithm (see Appendices A, B, C for details). Fig. 2(b) displays the result for different sizes L (and g 2 e /2 = t = 1), portraying the transition. Panels (a) and (c) display local configurations of matter n x and gauge sites Ê x,µ 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 2 m = 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 (see Appendix D) yields results compatible with a II-order phase transition, with the critical point occurying at negative bare masses m.
The same transition appears to be more interesting when we 'activate' the magnetic coupling, by setting g 2 m = 8t 2 /g 2 e = 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 ∼ g 2 m , kept as numerically relevant by the algorithm. Moreover, the transition is still compatible with a 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 2 e /2 = t = 1), ultimately making the transition physically relevant.

III. QUANTUM CAPACITOR
To investigate field-screening and equilibrium stringbreaking 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 c (d) = 2 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 cou-pling g 2 m = 0, analogously to the (1+1)D case [22]. In conclusion, the charged capacitor can make the phase transition physical even when g can not be tuned.
The observed behaviour can be interpreted as an equilibrium counterpart to the Schwinger mechanism, a realtime dynamical phenomenon in which the spontaneous creation of electron-positron pairs out of the vacuum is stimulated by a strong external electric field [43]. This could either be potentially verified in experiments or quantum simulations, by means of adiabatic quenches, ramping up the capacitor voltage.

IV. CONFINEMENT PROPERTIES
The (3+1)-dimensional pure compact lattice QED predicts a confining phase at large coupling g [11,[44][45][46][47]. 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 [48]. 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 [49].
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 2 m ∼ 0), as it can be extrapolated in the semiclassical limit (t ∼ 0). Here, a particleantiparticle pair at distance r with, a field-string between them, has an energy that scales linearly with r. On the contrary, 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 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 shortrange 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 rg 2 /2+ const.. Thus our simulations highlight visibly different features between confined and deconfined regimes, even with dynamical matter.

V. 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) fermionnumber symmetry, implemented in our TTN algorithms, we can inject any desired charge imbalance into the system, while working under OBC. Fig. 5 shows the results for charge density ρ = Q/L 3 = 1/4. In the vacuum phase (m g 2 e /2 ≈ t), we obtain configurations as displayed in panel (a), where the charges are expelled from the bulk,°2 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.

VI. OUTLOOK
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 [50] 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 selfvalidation tool: convergence of each quantity with the bond dimension can be verified in polynomial time.
From an experimental point of view, quantum link model formulations are the most studied pathway towards the simulation of LGTs on quantum hardware [41]. 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 the first experimental quantum simulations of one-dimensional LGTs [36][37][38][39][40]. In this framework, numerical methods capable of accessing intermediate sizes, such as TNs, play a fundamental role as a cross-verification toolbox.

VII. ACKNOWLEDGMENTS
Authors kindly acknowledge support from the Italian PRIN2017 and Fondazione CARIPARO, the Horizon 2020 research and innovation programme under grant agreementNo 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 by ATOS Bull.
In this section we present the construction of the QED gauge-invariant configurations for the local sites that we exploit as 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 are represented by three orthogonal eigenstates of the electric field operator: The parallel transporter, that in the spin language corresponds to the raising operator, acts on these states aŝ 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. [32] (see also [71,72]), generalizing it to the case with three spatial dimensions. This technique is related to the standard rishon formulation of QLM [73][74][75] and allows us to encode the 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 3-hardcore fermionic field η x,µ . Such lattice quantum fields satisfyη 2 x,µ = 0 and η 3 x,µ = 0. They mutually anticommute at different spatial positions, i.e. η x,µ ,η ( †) x ,µ = 0 for x = x or µ = µ , and also anticommute with the staggered matter fermionic fields η x,µ ,ψ ( †) x = 0 [76,77]. Then, we express the comparator on the link asÛ x,µ =η x,µη † x+µ,−µ . To explicitly build these 3-hardcore fermions for each half-link, we consider two species of standard Dirac fermionsâ x,µ andb x,µ and we use the following relation: wheren a x,µ andn b x,µ are the occupation number operators for each species, i.e.,n a x,µ =â † x,µâx,µ and the same for n b x,µ . For each 3-hardcore mode, these operators act on a three-dimensional local Hilbert space with basis |0 x,µ , x,µ |0 x,µ . In fact, due to the definition in Eq. (A3), the algebra of the operatorsη x,µ never accesses the fourth state obtained as b † x,µ |0 x,µ . By using the same representation on the other half-link through the Dirac operatorsâ † x+µ,−µ and b † x+µ,−µ , 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 aŝ L x,µ =n a x,µ +n b x,µ +n a x+µ,−µ +n b x+µ,−µ , (A4) defines a symmetry of the Hamiltonian since it commutes with the operatorsÊ x,µ andÛ x,µ . Thus, we can select the sector withL x,µ = 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 op-eratorÛ x,µ to act correctly following the properties of Eq. (A2). 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 the 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. (A5), it follows that, within the sector with the link-symmetry constraintL x,µ = 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: Due to the use of staggered-fermions, the presence/absence of a fermion in an even/odd site represents the presence of a charge/anti-charge.
Taking into account this property, it is possible to construct the gauge-invariant basis for the local site x, that is composed by the lattice site and the six half-links along the directions ±µ x , ±µ y , ±µ z (see Fig. 6): where |φ x = (ψ † x ) φ |0 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 (−1) δ k 1 ,2 +δ k 2 ,2 +δ k 3 ,2 allows us to satisfy the anticommutation relations of the fermionic representation recovering the correct signs of Eq. (A5). The occupation numbers φ and k j are not independent due to the constraint imposed by the Gauss's laŵ 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 Eq. (1) of the main text as 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 α select the right operators depending on the different directions in which the hopping process takes place. The operators M α x,µ 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 α, α , α , α depend on the plane of the plaquette (in this case x − y) and the links involved into the loop. The operators C α x 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 [71,72], so we can exploit standard TN algorithms, without the need of explicitly implementing the fermionic parity at each site [78][79][80].
The mass term and the electric field energy in the Hamiltonian of Eq. (1) 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 [81,82]. In conclusion, all the operators we employ in the TTN algorithms (see Appendix B) are genuinely local. In order to get an idea on the numerical complexity, we emphasize that the dimension of these matrices acting on the local gauge-invariant basis is 267 × 267.

Appendix B: Tensor Networks
In this section we present the main concepts of Tensor Networks (TNs) with a particular focus on the Tree Tensor Network (TTN) ansatz that we exploit in this work [83]. For a detailed and exhaustive description of the subject, please see the technical reviews and textbooks [19,84,85].
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 {|i k } 1≤i≤d . The whole Hilbert space of the system will be generated by the tensor product of the local Hilbert spaces, that is, H = H 1 ⊗ H 2 ⊗ · · · H N , with a resulting dimension equal to d N . Thus, a generic pure quantum state of the system |ψ can be expressed as a linear combination of the basis elements of H, i.e., In principle, the coefficients c i1,...,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 possibile to exploit rigorous results on the scaling of entanglement under a bipartition (area law) [86,87] 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. Tensor Networks provide a natural language for this representation [88,89] by decomposing the complete rank-N tensor c i1,...,i N in Eq. (B1) 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 amount 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 ).
MPS algorithms, such as the Density Matrix Renormalization Group (DMRG) [99], represent the state-ofthe-art technique for the numerical simulation of manybody 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 higherdimensional systems. PEPS intuitively reproduce 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 [100]. 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 ansazts for simulating quantum manybody systems are the TTNs, that 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 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 [101]. 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.
The TTN algorithm for the ground state computation of our LTG model follows the technical implementation described in [85] and it takes into account the conservation of the total charge through the definition of global U (1) symmetry sectors encoded in the TTN. In this way we can easily access finite charge-density regimes, with any 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, that 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 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 sec- tors. 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 [85]. 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 ef f 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 ef f is solved by using the Arnoldi method of the ARPACK library; (iii) the tensor is updated by the eigenvector of H ef f 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 desidered 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, BwUni-Cluster and ATOS Bull): a single simulation for the maximum size that we reached, a 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.

Appendix C: 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 Eq. (1) 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 Appendix A, in order to obtain the right representation of the electric field operators, we have to enforce the extra link symmetry constraint L x,µ = 2 at every pair of neighboring sites. For this reason, we include in the Hamiltonian additional terms that energetically penalise all the states with a number of hardcore fermions per link different from two, namely: where ν > 0 is the penalty coefficient and δ 2,Lx,µ 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 numerical level, this limit translates into choosing ν much larger than the other simulation parameters of the Hamiltonian, i.e., ν max {|t|, |m|, |g el |, |g m |}. However, setting ν too large in the first optimisation 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 a 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 of the energy: this effect signals that the global optmiza- tion procedure of the TTN become 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 desidered value of ν, that 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. 8(a) 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 L x,µ = 2, computed over all the links of the lattice on the ground state. The typical behavior of this quantity is shown in Fig. 8(b): 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 stability of our findings. Depending from the different system sizes and regimes of physical parameters, we estimate the relative error of the energy in the range 10 −2 , 10 −4 . A typical scaling of the energy density as a function of the inverse of the bond dimension 1/χ is shown in Fig. 8(c).

Appendix D: 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 2 m = 0, the Hamiltonian of Eq. (1) results diagonal in the local basis described in Appendix A and it is trivial to prove that the system undergoes a first-order phase transition between the bare vacuum, with energy E v = −m L 3 2 and the charge-crystal phase, with energy E ch = (m + that is obtained at E v = E ch . This behaviour is clearly seen in Fig. 9(a), showing a discontinuous transition between the two configurations.
In order to understand the behavior of the system for finite t = 1 and g 2 m = 0, we observe that the density, plotted in Fig. 2(a) 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 [102] 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 ρL β ν 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 these parameter and for each set of values we fit our points ρL β ν with an highdegree polynomial f L 1 ν (m − m c ) . 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 residual sum of squares 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 2 m = 8/g 2 e = 4, we check again the presence of a critical point and the values of critical exponents through the formula of Eq. (D1). 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. 9(b). Thus, while the transition and its univer-sality 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 anticharges 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. [103]).