Tensor network simulation of QED on infinite lattices: learning from (1+1)d, and prospects for (2+1)d

TL;DR

This paper advances tensor network methods for simulating lattice gauge theories, specifically QED, in the thermodynamic limit, starting from 1+1 dimensions and proposing a framework for 2+1 dimensions with gauge-invariant PEPS.

Contribution

It demonstrates gauge-invariant tensor network simulations of the Schwinger model and proposes a novel PEPS ansatz for (2+1)d QED incorporating gauge symmetry and matter fields.

Findings

Benchmarking of Schwinger model with good agreement

Intuitive insights for (2+1)d simulation strategies

Proposed PEPS ansatz includes gauge symmetry and matter fields

Abstract

The simulation of lattice gauge theories with tensor network (TN) methods is becoming increasingly fruitful. The vision is that such methods will, eventually, be used to simulate theories in $(3+1)$ dimensions in regimes difficult for other methods. So far, however, TN methods have mostly simulated lattice gauge theories in $(1+1)$ dimensions. The aim of this paper is to explore the simulation of quantum electrodynamics (QED) on infinite lattices with TNs, i.e., fermionic matter fields coupled to a $U(1)$ gauge field, directly in the thermodynamic limit. With this idea in mind we first consider a gauge-invariant iDMRG simulation of the Schwinger model -i.e., QED in $(1+1)d$-. After giving a precise description of the numerical method, we benchmark our simulations by computing the substracted chiral condensate in the continuum, in good agreement with other approaches. Our simulations of the Schwinger model allow us to build intuition about how a simulation should proceed in $(2+1)$ dimensions. Based on this, we propose a variational ansatz using infinite Projected Entangled Pair States (PEPS) to describe the ground state of $(2+1)d$ QED. The ansatz includes $U(1)$ gauge symmetry at the level of the tensors, as well as fermionic (matter) and bosonic (gauge) degrees of freedom both at the physical and virtual levels. We argue that all the necessary ingredients for the simulation of $(2+1)d$ QED are, a priori, already in place, paving the way for future upcoming results.

Figures

Tables

Table 1. Table 1: Simulation parameters in the one-site iDMRG algorithm.

Parameter	Description
${\chi }_c$	Bond dimension of charge index
${\chi }_d$	Bond dimension of degeneracy index
$p_{max}$	Gauge boson truncation
$N$	Number of added sites
$x$	Inverse coupling
$m/g$	Dimensionless fermion mass

Table 2. Table 2: Comparison: subtracted chiral condensate in the continuum. The extrapolation is in the regime x ∈ [ 10 , 300 ] 𝑥 10 300 x\in[10,300] .

$m/g$	One-site iDMRG	Ref.SchwingerDMRG	Ref.gaugeMPS	exact
0	0.15900	0.15993	0.15992	0.15992
0.125	0.09425	0.09202	0.09201	-
0.25	0.06838	0.06666	0.06664	-
0.5	0.04293	0.04238	0.04234	-
0.75	-	-	0.03062	-
1	-	-	0.02385	-
2	-	-	0.01246	-

Equations

$$\mathcal{L}=\overline{\psi }\left(i{\partial }_{\mu }{\gamma }^{\mu }-m\right)\psi -\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }-g\overline{\psi }{A}_{\mu }{\gamma }^{\mu }\psi ,$$

$$F^{\mu \nu }={\partial }^{\mu }{A}^{\nu }-{\partial }^{\nu }{A}^{\mu }.$$

$${\psi }^{\prime }=e^{ig\chi }\psi ,A_{\mu }^{\prime }=A_{\mu }+{\partial }_{\mu }\chi ,$$

$$\left\{{\gamma }^{\mu },{\gamma }^{\nu }\right\}=2g^{\mu \nu },$$

$${\gamma }^{\mu }\left(i{D}_{\mu }-m\right)\psi =0,$$

$${\partial }_{\mu }{F}^{\mu \nu }=g{j}^{\nu },$$

$$\left\{{\psi }^{†}\left(x,t\right),\psi \left(x,t\right)\right\}$$

$$\left\{{\psi }^{†}\left(x,t\right),{\psi }^{†}\left(x,t\right)\right\}$$

$$\left[E\left(x,t\right),A_1\left(y,t\right)\right]=i\delta \left(x-y\right),$$

$$E=-F^{01}=F^{10}.$$

$$\frac{\partial E}{\partial x}=g{j}^0\equiv g\rho ,\text{(Gauss’ law)}$$

$$-$$

$$\mathcal{H}=-i\overline{\psi }{\gamma }^1\left({\partial }_1-ig{A}_1\right)\psi +m\overline{\psi }\psi +\frac{1}{2}{E}^2,$$

$$H$$

$${\varphi }_n\leftrightarrow \{\begin{array}{l}{\psi }_{upper}n\text{ even}\\ {\psi }_{lower}n\text{ odd}\end{array},\psi =\left(\begin{array}{c}{\psi }_{upper}\\ {\psi }_{lower}\end{array}\right).$$

$${\theta }_n=-ag{A}_n^1.$$

$$g{L}_n=E_n.$$

$$e^{\pm i{\theta }_n}|p_n⟩=|p_n\pm 1⟩.$$

$$L_n-L_{n-1}={\varphi }_n^{†}{\varphi }_n-\frac{1}{2}\left(1-{\left(-1\right)}^n\right),$$

$$H$$

$$L_n-L_{n-1}=1/2\left({\sigma }_n^z+{\left(-1\right)}^n\right),$$

$$H$$

$${\psi }^{\prime }=e^{i\alpha }\psi \to {\mathcal{L}}^{\prime }=\mathcal{L},$$

$$j^{\mu }=\overline{\psi }{\gamma }^{\mu }\psi$$

$${\partial }_{\mu }{j}^{\mu }=0.$$

$${\psi }^{\prime }=e^{i\alpha {\gamma }_5}\psi .$$

$${\gamma }_0=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right),{\gamma }_1=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right),{\gamma }_5\equiv {\gamma }_0{\gamma }_1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right).$$

$$j_5^{\mu }=\overline{\psi }{\gamma }^{\mu }{\gamma }_5\psi .$$

$${\partial }_{\mu }{j}_5^{\mu }=\frac{g}{2\pi }{ϵ}_{\mu \nu }{F}^{\mu \nu },$$

$$\Sigma \equiv ⟨\overline{\psi }\psi ⟩$$

$${\Sigma }_0=\frac{e^{\gamma }}{2{\pi }^{\frac{3}{2}}}\approx 0.159929,$$

Full text

I introduction

Gauge field theories QFT are currently our deepest level of understanding of how fundamental interactions emerge from local symmetry principles. The standard model is a gauge theory, where different gauge symmetries orchestrate all known interactions except for gravity, which can be seen itself also as a gauge theory. The structure of gauge theories is so complex that, sometimes, it is wise to discretize them on a lattice in order to simulate their properties on a computer. Even if bumpy at its historical origins, the numerical simulation of lattice gauge theories lgt has become one of the main tools to understand our universe. This is particularly true for quantum chromodynamics (QCD), the theory of strong interactions, where lattice simulations allowed to, e.g., understand the spectrum of hadrons observed in particle accelerators.

Still, many questions concerning gauge theories remain open, and in particular for QCD. For instance, what is its phase diagram at finite fermionic density? Or what are the dynamical properties of the theory? Usual lattice gauge theory calculations, based mostly on quantum Monte Carlo, fail to answer such questions because of fundamental algorithmic limitations. Moreover, finite-size scaling of the results relies on accurate extrapolation laws to the thermodynamic limit which need to be somehow known beforehand.

In this setting, tensor network (TN) numerical methods tn have emerged as a promising alternative. In TN methods, the wavefunction of the system is decomposed into fundamental pieces, the tensors, glued together by quantum entanglement according to some network pattern. Such methods rely on correctly reproducing the amount and structure of entanglement in the wavefunction being simulated. The methods usually target low-energy properties, but can also be adapted to compute dynamics. Moreover, one can simulate both bosons and fermions with essentially the same computational cost ftn . And on top, gauge symmetries can be implemented naturally in this framework gtn ; gaugeMPS . So all in all, TNs look like the natural option to describe the structure of quantum states present in lattice gauge theories.

Our aim with this paper is to pave the way towards higher-dimensional numerical simulations of lattice gauge theories, in particular for $(2+1)d$ quantum electrodynamics (QED), i.e., the gauge theory for electromagnetism. In order to build intuition, we first do a detailed analysis of simulations of lattice QED in $(1+1)d$, the so-called Schwinger model schwinger , using gauge-invariant Matrix Product States (MPS) mps ; gaugeMPS and a gauge-invariant version of infinite Density Matrix Renormalization Group (iDMRG) dmrg ; idmrg . This allows us to foresee how a higher-dimensional simulation should proceed. For the higher-dimensional case we discuss briefly the lattice Hamiltonian, and give a proposal for a 2d TN ansatz based on infinite Projected Entangled Pair States (iPEPS) PEPS ; iPEPS . As we shall see, such an iPEPS implements naturally fermionic matter and $U(1)$ gauge bosons. Thinking in perspective, we argue that all the necessary ingredients for a TN simulation of QED in $(2+1)d$ are a priori already there.

Previous results on the TN simulation of lattice gauge theories include a number of works. $\mathbb{Z}_{2}$ lattice gauge theories in $(1+1)d$ have been considered with DMRG sugihara . For the Schwinger model, DMRG (without MPS formulation) was considered in several works dmrgOld , whereas MPS simulations have been done to compute the chiral condensate SchwingerDMRG as well as thermal properties SchwingerThermal , the mass spectrum SchwingerSpectrum , the Schwinger effect SchwingerSchwinger , the effect of truncation in the gauge variable SchwingerTruncation , and the case of several fermionic flavours SchwingerMultiflavour . Gauge invariance in the MPS of the Schwinger model was originally considered in Ref.gaugeMPS , where the ground state was computed using the time-dependent variational principle (TDVP) tdvp . Gauge-invariant MPS were used to compute the confining potential SchwingerPot . A similar approach was used to analyze the scattering of two quasiparticles and the dynamical generation of entanglement SchwingerScat . TN simulations have also been implemented recently for non-abelian lattice gauge theories in $(1+1)d$ nonabelianTNS . For higher-dimensional systems, gauge-invariant TN ansatzs have also been proposed analytically gtn ; u1peps ; su2peps .

This paper is organized as follows: first, in Sec.II we provide a detailed introduction to QED in $(1+1)$ dimensions (the Schwinger model) in the continuum as well as its discretized version on the lattice. In this section we provide also background on the so-called chiral condensate. Then, in Sec.III we revise the infinite-DMRG algorithm. We discuss the details of the variational optimization of MPS, with one-site and two-site updates in the thermodynamic limit. In Sec.IV we explain how to do a gauge-invariant simulation of the Schwinger model using infinite DMRG. Numerical benchmarks for the chiral condensate in $(1+1)d$ are presented in Sec.V, paying attention to the continuum limit extrapolation. Based on all this, in Sec.VI we discuss the prospects for the simulation of QED in $(2+1)d$, where we consider the lattice formulation of the Hamiltonian as well as a possible TN ansatz for its ground state in terms of a 2d infinite PEPS. Finally, Sec.VII contains our conclusions.

II QED in $(1+1)d$: the Schwinger model

Let us now revise the basics of QED in $(1+1)$ dimensions, also called QED2, or the Schwinger model schwinger . We will refresh some of the properties of this theory defined in the continuum, as well as a possible formulation on the lattice, which will be the starting point of our study with TN methods. Readers who are interested in a more detailed discussion of the model and its properties are referred to, e.g., Ref. SchwingerDMRG .

II.1 Continuum formulation

The massive Schwinger model is quantum electrodynamics in two space-time dimensions. Its Lagrangian density in the continuum reads

[TABLE]

where

[TABLE]

The first term is the Dirac Lagrangian density for a free fermion and the second term corresponds to the field energy of the electric field (in $(1+1)d$ there is no “room” for a magnetic field). The third term is the interaction between the matter field and the gauge field. It has the important feature that it arises from the constraints imposed by a local gauge transformation. That means, its form is determined by demanding the invariance of the Lagrangian density under the transformation

[TABLE]

where $\chi$ is an arbitrary real function of space and time 111This is what is meant by local: one can redefine the phase of the fermionic field locally at every point in space-time., i.e. $\chi=\chi\left(x,t\right)$. The Schwinger model describes the interaction of one flavour mass-$m$ fermions $\psi$ with a $U(1)$ gauge field $A$, with coupling $g$. In $(1+1)d$ the Lorentz indices $\mu,\nu$ run from [math] to $1$ (one direction for space, and one for time), and the gamma matrices satisfy the Clifford algebra

[TABLE]

analogously to the $(3+1)d$ case. However, since there is no spin degree of freedom in one spatial dimension, these are $2\times 2$ matrices. Substituting the Lagrangian of the Schwinger model into Euler-Lagrange equations for the fields $\psi$ and $A$ results in the equations of motion

[TABLE]

and

[TABLE]

where $D_{\mu}\equiv\partial_{\mu}+igA_{\mu}$ is the gauge covariant derivative and $j^{\nu}\equiv\overline{\psi}\gamma^{\nu}\psi$ the vector current. The theory is quantized using canononical quantization by imposing anti-commutation relations on the fermion fields

[TABLE]

and by imposing commutation relations on the gauge fields

[TABLE]

where the electric field $E$ is defined by

[TABLE]

Using this definition of the electric field in Eq.(6), we get analogues to Maxwell’s equations in $(1+1)d$:

[TABLE]

Since there is “no space” for magnetic fields in one spatial dimension, we only obtain the analogue of Gauss’ law and an equation which describes the dynamics of the electric field.

II.2 Lattice formulation

Starting from the Hamiltonian density $\mathcal{H}$ in temporal gauge, $A_{0}=0$,

[TABLE]

the model can be formulated on a spatial lattice using a Kogut-Susskind staggered formulation KogutSusskind . The equivalent lattice Hamiltonian is

[TABLE]

where $a$ denotes the lattice spacing. In this formulation the correspondence between the fermionic lattice field $\phi_{n}$ on site $n$ and the continuum field $\psi$ is

[TABLE]

The gauge variables $\theta_{n}$ live on the links between the sites $n$ and $n+1$, and are related to the vector potential via

[TABLE]

Their conjugate variables $L_{n}$, with $\left[\theta_{n},L_{m}\right]=i\delta_{nm}$, are related to the electric field by

[TABLE]

Since $\theta_{n}$ is an angular variable, $L_{n}$ will have integer charge eigenvalues $p_{n}\in\mathbb{Z}$. Therefore, the local Hilbert space spanned by the corresponding eigenvectors $\left|p_{n}\right\rangle$ is infinite, and $e^{\pm i\theta_{n}}$ are the ladder operators

[TABLE]

The lattice equivalent of Gauss’ law then reads

[TABLE]

which means that excitations on odd and even sites create $\mp 1$ units of flux, corresponding to “electron” and “positron” excitations, respectively. Using a Jordan-Wigner transformation, $\phi_{n}=\Pi_{k<n}\left(i\sigma_{k}^{z}\right)\sigma_{n}^{-}$, where $\sigma^{\pm}=\frac{1}{2}\left(\sigma^{x}\pm\sigma^{y}\right)$, the fermionic degrees of freedom can be mapped to spin-1/2 degrees of freedom while keeping the Hamiltonian local, i.e.,

[TABLE]

In the above equation we introduced the parameters $x\equiv 1/\left(g^{2}a^{2}\right)$ and $\mu\equiv 2\sqrt{x}m/g$. The spins live on the sites of the lattice, with $\sigma_{n}^{z}\left|s_{n}\right\rangle=s_{n}\left|s_{n}\right\rangle$, and represent “positrons” on even sites and “electrons” on odd sites. An even site with $s_{2n}=-1$ corresponds to an empty positron state, while $s_{2n}=1$ represents an occupied positron state, and vice versa for the odd electron sites.

In $(1+1)d$, Gauss’ law can therefore be rewritten as

[TABLE]

and can in fact be used to remove the gauge degrees of freedom SchwingerDMRG . The resulting lattice Hamiltonian is then

[TABLE]

where $l$ is a possible external background charge. In this new formulation there are no gauge variables but, however, we pay the price of having a non-local, long-range interaction term in the Hamiltonian.

II.3 Chiral condensate in the continuum

Let us now revise the so-called chiral condensate. Without attempting to go into detail, we discuss two continuous symmetries of the Schwinger model of which one is broken after quantization. In this context, the chiral condensate arises as an order parameter.

The Lagrangian density of the Schwinger model is invariant under global phase transformations of the Dirac field, i.e.,

[TABLE]

where $\alpha$ is a real constant. According to Noether’s theorem (see, e.g., Ref.QFT ), there is a conserved current $j^{\mu}$ associated with every continuous symmetry. In this case the vector current

[TABLE]

is conserved, i.e.,

[TABLE]

This global $U(1)$ symmetry is known to hold in fermionic field theory models although, in principle, the vacuum state could spontaneously break it Chiral1 .

Let us now consider the case of massless ($m=0$) fermions. Then the Lagrangian of the Schwinger model has another continuous symmetry, namely the so-called chiral symmetry. This symmetry implies that the Lagrangian density is invariant if one transforms $\psi$ into $\psi^{\prime}$ as

[TABLE]

In the above equation $\gamma_{5}\equiv\gamma_{0}\gamma_{1}$ anti-commutes with $\gamma_{\mu}$ for $\mu=1,2$ and $\alpha$ is again a real constant. For example, in the Dirac representation the gamma matrices are given by 222See, e.g., Ref.Chiral1 .

[TABLE]

The associated Noether current for this symmetry is the so-called axial-vector current $j_{5}^{\mu}$ which is given by

[TABLE]

While the vector current in Eq.(22) is conserved in the quantized theory, the axial-vector current is not. This non-conservation of the axial-vector current is called chiral anomaly or axial anomaly. The divergence of the axial-vector current reads

[TABLE]

where $\epsilon_{\mu\nu}$ is the Levi-Civita symbol in two dimensions, see Refs.Chiral1 ; Chiral2 . As a consequence of this chiral symmetry breaking, the vacuum expectation value

[TABLE]

becomes non-zero. The quantity $\Sigma$ is called chiral condensate Chiral1 . In the case of the massless Schwinger model, the chiral condensate can be computed exactly (see, e.g., Ref.Wipf ), and is found to be

[TABLE]

where $\gamma$ is the Euler-Mascheroni constant. Therefore, the chiral condensate can be regarded as an order parameter signalling chiral symmetry breaking in the vacuum.

II.4 Chiral condensate on the lattice

In the lattice formulation of the massive Schwinger model, it is possible to write the chiral condensate in terms of Pauli spin operators. It is easy to see that this reads

[TABLE]

where the expectation value is computed in the ground state and $N$ is the number of lattice sites. The naively-computed chiral condensate is known to be UV-divergent. In particular, it diverges logarithmically in the continuum limit, $a\rightarrow 0$. It has been argued that this divergence comes solely from the free theory at $g=0$. In the free case the chiral condensate on the lattice $\Sigma_{free}\left(x\right)$ can be computed exactly as

[TABLE]

where $K\left(z\right)$ is the complete elliptic integral of the first kind. This result can be used to subtract the divergence from the computed chiral condensate in the interacting theory, and therefore to renormalize it. In other words, we can define a so-called subtracted chiral condensate $\Sigma_{sub}$, which allows for a continuum extrapolation, by

[TABLE]

where $\Sigma\left(x\right)$ denotes the computed chiral condensate. Details on the extrapolation procedure for numerical data will be given in the next chapters.

III Infinite DMRG

Here we review the basics of the DMRG algorithm for infinite systems. Several formulations of this algorithm have been proposed in the literature. The approach taken here is similar to that in the second paper of Ref.idmrg , where we consider both the case of one-site and two-site updates. For the sake of completeness, we also briefly review some of the basics of variational optimization algorithms over tensor networks tn .

III.1 MPS variational optimization

In general terms, we want to approximate the ground state of a Hamilitonian expressed as a Matrix Product Operator (MPO) by minimizing

[TABLE]

over the family of Matrix Product States (MPS) with bond dimension $\chi$. This can be achieved by introducing a Lagrange multiplier $\lambda$ that enforces normalization, so that the minimization reads

[TABLE]

The above minimization is performed by adjusting all tensors in the MPS for all sites in order to make the expectation value of the energy the lowest possible. In DMRG one follows a sequential approach, optimizing tensor by tensor. In terms of the chosen tensor, which we call $A$, the mimization problem defined by Eq.(34) can be written as

[TABLE]

In the above equation, all coefficients of $A$ are arranged as a vector $\vec{A}$ as shown in Fig.1(a), $\mathcal{H}_{eff}$ is an effective Hamiltonian, and $\mathcal{N}$ is a normalization matrix. The effective Hamiltonian and the normalization matrix can be considered as the environment of tensors $A$ and $A^{*}$ in the two TNs for $\left\langle\psi\right|H\left|\psi\right\rangle$ and $\left\langle\psi\right.\left|\psi\right\rangle$ respectively, but written in matrix form (see e.g. Fig.1(b)).

The minimization condition

[TABLE]

leads to the generalized eigenvalue problem

[TABLE]

Once this optimization with respect to $A$ is done, one proceeds by repeating the minimization for another tensor in the MPS. In this way, one continues sweeping through all tensors several times, until the desired convergence in expectation values is attained. Let us remark that if we start from an MPS with open boundary conditions, this algorithm is nothing else but the Density Matrix Renormalization Group (DMRG) algorithm in the language of TNs dmrg ; tn . In the case of open boundary conditions it is also always possible to choose an appropriate gauge for the tensors, e.g., a mixed canonical form with $A$ as the center site, such that $\mathcal{N}=\mathbb{I}$. Then Eq.(37) reduces to an ordinary eigenvalue problem. This is very useful for practical implementations since it avoids stability problems due to $\mathcal{N}$ being ill-conditioned, see Ref.tn . In what follows, we always consider MPS with open boundary conditions in mixed canonical form.

III.2 One-site infinite DMRG

If we start from the very beginning with an infinite system to study systems in the thermodynamic limit, we need to modify the above procedure idmrg . Let us assume that we were given an infinitely large and translationally invariant system in its ground state. Then, if we were to add an additional site to the system and allow it to relax, one would expect that the new site would change to match the rest, while the other sites in the system remain essentially unchanged. In MPS language, let us consider the case in which we already have an infinite MPS with bond dimension $\chi$ representing the ground state of our system. Then adding a site to our system would correspond to adding another tensor in the MPS. The relaxation process could be simulated by minimizing the energy with respect to the new tensor in the environment given by the original MPS. We would then obtain a tensor which looks, mostly, like all of the tensors in our infinite MPS. The idea of the algorithm is to start with a representation of the infinite system in terms of an approximative environment. This environment is then progressively refined by embedding new sites, allowing the sites to relax, and then absorbing them into the environment. Eventually this procedure will converge, thus simulating the environment experienced by a single site in the infinite system in its ground state.

The infinite-system algorithm thus works as follows: one starts from initial environments $L_{H}$ and $R_{H}$ (e.g., random) representing the left and right halves of the (infinite) system with respect to the added tensor $A$ of the TN for $\left\langle\psi\right|H\left|\psi\right\rangle$ (see Fig.2(a)). Then, one iterates the following procedure:

Relaxation: compute the eigenvector $\vec{A}$ corresponding to the minimal eigenvalue of the problem 333We choose $A$ as the center site for the mixed canonical form of the MPS. $\mathcal{H}_{eff}\vec{A}=\lambda\vec{A}$, and reshape it back to a 3-index tensor. The effective Hamiltionian is shown in Fig.2(b). 2. 2.

Absorption (odd step): at an odd iteration step, the optimized tensor is absorved into the left environment $L_{H}$. In detail:

1. (a)

Merge the first bond index and the physical index of $A$ to form a matrix, and compute the singular value decompositon $A=U\Sigma V^{\dagger}$ (see Fig.3(a)). 2. (b)

Undo the index fusion for the left index of $U$ to get back to a 3-index tensor (see Fig.3(a)) and compute $E_{H}$ as defined in Fig.3(b). 3. (c)

Refine the approximation for the left environment $L_{H}$ by contracting $E_{H}$ into it, i.e. $L_{H}\leftarrow L_{H}\cdot E_{H}$, as shown in Fig.3(c). 3. 3.

Absorption (even step): at an even iteration step, the optimized tensor is analogously contracted into the right environment $R_{H}$ (see Fig.4). In detail:

1. (a)

Merge the second bond index and the physical index of $A$ to form a matrix, and compute the singular value decompositon $A=U\Sigma V^{\dagger}$. 2. (b)

Undo the index fusion for the right index of $V^{\dagger}$ to get back to a 3-index tensor and compute the analogue of the tensor $E_{H}$. 3. (c)

Refine the approximation for the right environment $R_{H}$ by contracting $E_{H}$ into it, i.e. $R_{H}\leftarrow E_{H}\cdot R_{H}$.

Since $U$ and $V$ are isometries, the mixed canonical form of the MPS is preserved at every simulation step. To check for convergence it is useful to calculate the desired expectation value after, e.g., each first or second iteration step. For a single-site operator acting on the added site this can be done easily, thanks specially to the mixed canonical form of the MPS. The main computational cost is given by the eigenvalue problem and scales therefore as $O\left(\chi^{3}\right)$.

III.3 Two-site infinite DMRG

If only a single site is added at every iteration step, then the MPS bond dimension $\chi$ is fixed right from the start in the algorithm. However, one may think of situations in which it would be advantageous to increase the bond dimension during the calculation. This can be done by a slight modification of the algorithm, namely, by adding two sites at each iteration, see Fig.5. The two-site infinite DMRG algorithm is then as follows:

Relaxation: compute the eigenvector $\vec{\Theta}$ corresponding to the minimal eigenvalue of the problem $\mathcal{H}_{eff}\vec{\Theta}=\lambda\vec{\Theta}$, where the effective Hamiltionian $\mathcal{H}_{eff}$ and the vector $\vec{\Theta}$ are defined as shown in Fig.5(c) and Fig.5(b), respectively. 2. 2.

Absorption: the optimized tensor is simultaneously contracted into the left environment $L_{H}$ and into the right environment $R_{H}$. In detail:

1. (a)

Compute the singular value decomposition $\Theta=U\Sigma V^{\dagger}$ (see Fig.5(d)) 2. (b)

Undo the index fusion for the left index of $U$ and for the right index of $V^{\dagger}$. 3. (c)

Compute the tensors $E_{HL}$ and $E_{HR}$ as defined in Fig.5(e). 4. (d)

Refine the approximations for the left environment $L_{H}$ and for right environment $R_{H}$ by the contractions $L_{H}\leftarrow L_{H}\cdot E_{HL}$ and $R_{H}\leftarrow E_{HR}\cdot R_{H}$ shown in Fig.5(f).

The crucial point is that, if one adds two sites at a time, then the central matrix becomes a square matrix of increased dimension $d\chi\times d\chi$ as can be seen in Fig.5(b). This allows, in principle, for a truncation of the SVD in the second iteration step, see also Fig.5(d), which allows the bond dimension to grow as the algorithm proceeds. This is also particularly useful if one implements symmetries in the algorithm, since the truncation allows to change the symmetry sectors being kept at every step. In practice, this means that the algorithm can readapt itself to more relevant symmetry sectors, which have more weight in terms of the singular values of $\Theta$, leading to improved accuracy.

IV Gauge-invariant infinite DMRG

Following the ansatz introduced in Ref.gaugeMPS , we start one step before integrating out the gauge field degrees of freedom using the Gauss’ law constraint. This is, we consider the Hamiltonian in Eq.(18), with spin variables on the sites for the staggered fermionic matter field, and angular variables on the links for the bosonic gauge (electric) field. An obvious advantage is that this Hamiltonian is local with at most nearest-neighbour actions, and translationally invariant under shifts by two sites. Furthermore, it is practical for possible generalizations to higher-dimensional systems, since the gauge degrees of freedom can only be integrated out in $(1+1)d$.

IV.1 MPO representation of $H$

In the following, we provide an MPO representation of the Hamiltonian in Eq.(18) to be used in an iDMRG simulation. For convenience, we block site $n$ and link $n$ into a single MPS-site, such that at every MPS-site we have a fermionic and a gauge field degree of freedom. Then, the Hamiltonian can be regarded as the sum of 1-site operator and 2-site operators, i.e,

[TABLE]

where

[TABLE]

and

[TABLE]

The first factor in the tensor product $\otimes$ refers to the fermion degree of freedom, and the second to the gauge field degree of freedom at the MPS site. With $\cdot$ we denote here the tensor product between operators acting on different MPS-sites. The Hamiltonian can be written as an MPO with bond dimension $D=4$ where non-zero coefficients of the tensors are given as in Fig.6.

IV.2 Imposing gauge invariance

We now impose gauge invariance to enforce that our algorithm works directly within the physical subspace of the full Hilbert space. In particular, we are only interested in states $\left|\psi\right\rangle$ that are gauge invariant, i.e.,

[TABLE]

where

[TABLE]

Eq.(41) is nothing but the (discretized) lattice version of the Gauss’ law constraint for the system.

A possibility to impose gauge invariance would be to add a penalty term to the Hamiltonian, so that the gauge-invariant subspace is energetically preferred. For example, one could consider the modified Hamiltonian

[TABLE]

instead of $H$, and then take the ground state sector in the limit $\lambda\rightarrow\infty$. However, by doing this gauge invariance would only be approximately realized, and one would have to extrapolate additionally in parameter $\lambda$.

A safer and more direct option is to implement gauge invariance directly at the level of the tensors in the TN, i.e., consider a TN made of $U(1)$ gauge-symmetric tensors gaugeMPS . This implies that many tensor components in the MPS ansatz must vanish, i.e., only components compatible with gauge symmetry can be different from zero symTN .

For the sake of concreteness, let us assume that we have a finite lattice of $N\in 2\mathbb{N}$ sites. Then a general, i.e., not necessarily gauge invariant, MPS ansatz for the system has the form

[TABLE]

where the matrices $(B_{n}^{s_{n}})_{\alpha\beta}$ correspond to fermionic degrees of freedom, and the matrices $(C_{n}^{p_{n}})_{\alpha\beta}$ to gauge degrees of freedom. We denote the bond dimension with $\chi$, i.e., the bond indices take the values $\alpha,\beta=1,\ldots\chi$.

From Eq.(41) and Eq.(42) we can see that Gauss’ law is basically a prescription to update the electric field $L_{n}$ at the right link of site $n$, namely,

[TABLE]

Therefore, if there is no charge at the site $n$, then $L_{n}$ stays with the value $L_{n-1}$ at the left. At the same time the electric field $L_{n}$ is increased/decreased by one unit if there is a positron/electron 444Recall that an occupied positron or electron state corresponds to $s_{n}=1$ or $s_{n}=-1$, respectively, and that positrons/electrons live on even/odd sites $n$. at site $n$. This “update rule” can be implemented by giving the bond indices a multiple index structure, $\alpha\rightarrow\left(q,\alpha_{q}\right)$, and imposing the following form on the tensors in the bulk:

[TABLE]

If one chooses a vanishing electric field to the left of the first lattice site, i.e. $L_{0}=0$, then the tensors representing the boundaries are gauge invariant if:

[TABLE]

In the above equations, the indices $q$ and $r$ label the electric charge sector, and are sometimes referred to as structural or charge indices. They label the representation of the gauge symmetry group for the index, and run from $1$ to a structural bond dimension $\chi_{c}$. The indices $\alpha_{q}$ and $\beta_{r}$ label the degeneracy subspace within each charge (symmetry) sector, and run from $1$ to a degeneracy bond dimension $\chi_{d}$. Every bulk or boundary tensor which is chosen according to Eq.(46) or Eq.(47), respectively, preserves the gauge symmetry exactly. The variational freedom lies now within the matrices $b_{n,q}^{s_{n}}$ and $c_{n}^{p_{n}}$, and the total MPS bond dimension is given by $\chi=\chi_{c}\cdot\chi_{d}$. The rather lengthy derivation of the result can be found in Ref.gaugeMPS . We also refer the reader to Ref.symTN for details on the implementation of symmetries in TNs and its consequences.

IV.3 Further details

The strategy presented above is very general. For an iDMRG simulation, the MPO is itself gauge-invariant by construction. If the MPS ansatz is also gauge-invariant, then the whole algorithm preserves gauge symmetry at every iteration step, of course provided that the initial conditions for the left and right environments are also gauge-invariant. This initial condition for the environment tensors is very easy to impose.

Since our main goal is to learn from the simulations in $(1+1)d$, we use a coarse-grained version of the one-site iDMRG algorithm presented previously to find a ground state approximation in the thermodynamic limit (let us mention that we also tested a non-coarse-grained version of the two-site iDMRG algorithm, leading to essentially equivalent results). As in the construction of the Hamiltonian, we again block a lattice site and a link into one MPS-site. This leads to an MPS ansatz with a two-site unit cell due to alternating spin-gauge systems for positrons and electrons. The initial tensors are defined according to Eq.(46), but are otherwise chosen randomly (or according to some educated guess) within the variational gauge-invariant subspace. In order to obtain a system that is invariant under translations of one site, we also block neighbouring MPS-sites corresponding to a positron and electron spin-gauge systems together. This procedure is shown in Fig.7.

A list of all the relevant simulation parameters is shown in Table 1.

V Results

V.1 Numerical benchmarks

We computed the chiral condensate for four different values of the fermion mass, $m/g=0,0.125,0.25,0.5$ where in each of the cases we took many points in the interval $x\in\left[10,600\right]$. Such a large interval allows us to extrapolate to the continuum limit, as well as to see the effect of the finite bond dimension as this limit is approached. The parameter $p_{max}\geq\left|p_{n}\right|$ truncates the infinite local Hilbert space of the gauge bosons, and amounts to a maximum bosonic occupation number. Physically it can also be seen as truncation in the gauge $U(1)$ charge. In our calculations we choose $p_{max}=2$, i.e., we truncate the infinite dimensional Hilbert space to five dimensions 555$L_{n}\left|p_{n}\right\rangle=p_{n}\left|p_{n}\right\rangle$ with $p_{n}=0,\pm 1,\pm 2$.. In practice we have seen that this truncation is sufficient for our purposes 666A more detailed analysis and justification of this truncation can be found in Ref.SchwingerTruncation .. Furthermore, we set $N=500$ which corresponds to adding 1000 sites in the physical system due to the two-site coarse-graining. At every simulation step we check that the expectation value of the gauge operator $G_{n}$ defined in Eq.(42) is zero, as required by gauge invariance.

Importantly, in order to get an approximation of the (subtracted) chiral condensate in the continuum, we have to perform a sequence of extrapolations. First, for every $x$ we extrapolate to infinite MPS bond dimensions $\chi_{c}$ and $\chi_{d}$. Second, for the extrapolated curve as a function of $x$, we extrapolate to the continuum limit so that $1/\sqrt{x}\rightarrow 0$.

Let us show an example of the extrapolation in the MPS bond dimensions for $m/g=0.25$. In Fig.8 we show the substracted chiral condensate as a function of $1/\sqrt{x}$ for different bond dimensions $\chi_{c}$ and $\chi_{d}$. One can see that the effect of the truncation becomes stronger as the lattice parameter $a$ becomes smaller, i.e., in the region tending towards the continuum limit. This is an important observation: the closer we are to the continuum limit, the harder the simulation becomes. It may be possible to simulate the lattice system always in an “easy” regime far from the field theory limit, but it is important to remember that in such a case we would not be simulating a field theory, but rather some (interesting but discrete) lattice spin model. In our simulations, the results for the substracted chiral condensate seem to be well converged over the chosen spectrum of bond dimensions.

In practice, for every $x$ we do an extrapolation in the total MPS bond dimension, i.e., $\chi\equiv\chi_{c}\cdot\chi_{d}\rightarrow\infty$. We find that the dependency of the chiral condensate is well described by the fitting function

[TABLE]

where $\Sigma\left(x,\chi=\infty\right)$ is the value extrapolated to infinite bond dimension for inverse coupling $x$. in Fig.9 we show some of these extrapolations for $x=100,250$ and $500$.

Finally, the extrapolation to the continuum limit $a\rightarrow 0$ is taken by considering the regime $x\rightarrow\infty$. Following the procedure in Ref.SchwingerDMRG , we fit the subtracted chiral condensate using the following ansatz:

[TABLE]

where $\Sigma_{sub}$ is the extrapolated continuum value of the subtracted chiral condensate. In Fig.10 we can see that this ansatz describes our data overall very well. This is especially true in the case of larger lattice constants, where the influence of the finite bond dimension is also smaller and the results are therefore easier to converge. Using this fit, we perform a continuum extrapolation for $x\in\left[10,300\right]$, where the convergence of our algorithm is particularly good. The obtained results for the four different fermion masses can be found in Table 2.

V.2 Discussion

As one can see in Table.2 our results are in agreement with the results in Ref.SchwingerDMRG and Ref.gaugeMPS for a fitting region $x\in[0,300]$, which is very well converged. Notice, though, that the approach in Ref.SchwingerDMRG is conceptually very different, since it is based on finite-size DMRG calculations using the non-local Hamiltonian from Eq.(20). In our gauge-invariant iDMRG approach, however, we start from the local Hamiltonian in Eq.(18). We think that this approach is more convenient in order to generalize the calculations to higher-dimensional systems, since it preserves locality explicitly and is therefore more amenable to, e.g., approaches based on infinite Projected Entangled Pair States (iPEPS) iPEPS . This is particularly true, also because in higher dimensions the Gauss’ law cannot be integrated out, and therefore the most natural option is to work with a gauge-invariant 2d PEPS targeting a 2d local Hamiltonian on the lattice, as we shall discuss in Sec.IV.

Let us also stress that in our $(1+1)d$ calculations with iDMRG, the bond dimensions did not need to be too large in order to get decent results. In particular, for $m/g=0.25$ we used $\chi=91$ as the highest total bond dimension, while in the other cases it was $\chi=63$. The extrapolations in Ref.SchwingerDMRG were attained from calculations up to bond dimension $\chi=140$, though via a different algorithm (as mentioned above). Remarkably, relatively small bond dimensions in iDMRG already allows us to provide results which are in quite good agreement with the ones for large bond dimension in Ref.SchwingerDMRG , on top of not having to do any finite-size extrapolation since we work directly in the thermodynamic limit.

For further comparison, in Ref.gaugeMPS , besides working with gauge invariant MPS in the thermodynamic limit, the authors also exploited CT symmetry, i.e., invariance by a one-site translation and charge conjugation. The ground state calculations were done via the so-called time-dependent variational principle (TDVP) tdvp . In this work symmetries were treated in a more sophisticated way by distributing variational freedoms to different charge sectors. In contrast to that, our approach here is simpler, since we just choose gauge invariant initial tensors and then let the algorithm evolve, which naturally preserves gauge invariance. As such, it is remarkable that our simple approach produces results which are also in qualitative agreement with those produced by more sophisticated methods.

Moreover, we remind that here we used the one-site iDMRG algorithm. Despite being more efficient, we know that a two-site iDMRG calculation would bring some extra advantages, e.g., a dynamical increase of the bond dimension, and a dynamical truncation of the gauge-symmetry sectors. Still, we run some checks with a 2-site algorithm but did not obtain much greater accuracy in the regimes explored in this paper. However, it is good to keep in mind that the two-site approach may still be very useful in the more entangled regimes.

VI Prospects for QED in $(2+1)d$

Taking into account what we have learned in the simulation of the $(1+1)d$ case, we would like now to consider the possibility of simulating the lattice version of QED in $(2+1)d$, directly in the thermodynamic limit. This gauge theory is interesting for a number of reasons: it is “closer” to our $(3+1)d$ space-time and also has confinement confQED3 which, unlike in the case of the Schwinger model, appears through a mechanism much more similar to the one in $(3+1)d$ QCD QFT .

Simulating first the Schwinger model has allowed us to learn a number of useful things about how the simulation in $(2+1)d$ should proceed. In particular, for the $(2+1)d$ case one needs to face the following facts:

Gauss’ law cannot be explicitly integrated out. Therefore, the safest choice is to work with a TN of gauge-invariant tensors. 2. 2.

The Jordan-Wigner transformation in $(2+1)d$ introduces non-local strings when mapping some fermionic terms into spins. Therefore, it is more convenient to work directly in fermionic Fock space. 3. 3.

Additionally to the electric field, in $(2+1)d$ there is a also magnetic field term which, in the lattice formulation, corresponds to a plaquette energy term in the Hamiltonian. 4. 4.

Moreover, and as in the Schwinger model, the gauge-boson Hilbert space should be truncated in a maximum occupation number in order to do the simulation (quantum link model) qlm .

Considering the above, and following the intuition built from the simulation of the $(1+1)d$ case, we would therefore need the following ingredients for the $(2+1)d$ simulation:

A TN in 2d as a variational ansatz in the thermodynamic limit. The so-called infinite-PEPS is the most natural option iPEPS . 2. 2.

The ability to simulate fermions in 2d. This has already been achieved, with fermionic implementations of the iPEPS algorithm ftn . 3. 3.

The ability to implement $U(1)$ gauge symmetry in the tensors. This has also been done already in 2d PEPS gtn ; u1peps . 4. 4.

The ability to deal with plaquette interactions. This has also been done in the past for iPEPS, e.g., when simulating the Toric Code model in a field and its generalizations plaqpeps . 5. 5.

Efficient and accurate optimization strategies. Regarding this, important developments in 2d iPEPS methods have been put forward recently optpeps .

We conclude, therefore, that a priori all the necessary ingredients for this simulation are already available. In the following section we would like to be a little bit more specific on how such a simulation could proceed.

VI.1 Lattice formulation

The Hamiltonian of QED in $(2+1)d$ on a lattice can be derived in a similar way as the one in $(1+1)d$ in Sec.II, but taking into account that this time one has two spatial dimensions instead of one. We give here a lattice Hamiltonian that has the correct continuum limit 777See, e.g., Ref.dmrgOld for more details.. On a $2d$ spatial square lattice, the Hamiltonian is given by

[TABLE]

In the equation above, $a$ is again the lattice spacing, $m$ the mass of the fermionic field, and $g$ the coupling between fermonic matter and the gauge boson. Fermionic fields are again staggered, but this time on a 2d square lattice, i.e., along both spatial directions, see Fig.11. The gauge boson variables $\theta_{n,m}$ live on the link between sites $n$ and $m$, and the sum $\langle n,m\rangle$ runs over nearest neighbours. The factor $s(n)$ decides the $+1$ or $-1$ prefactor for the mass term depending on the staggered pattern of the fermionic field: $+1$ for positrons, and $-1$ for electrons. Finally, the term with the cosinus is the curl of the gauge variable around a plaquette, see Fig.11, and corresponds therefore to the magnetic field energy.

In this setting, the Gauss’ law in $(2+1)d$ reads

[TABLE]

where the first equation is for horizontal links, and the second for the vertical. Finally, in order to implement a simulation, it is advisable to truncate again the local dimension of the Hilbert space of the gauge boson, as we did in the $(1+1)d$ case.

VI.2 Variational ansatz: a proposal

As a variational TN ansatz to approximate the ground state of the above Hamiltonian we propose a 2d infinite PEPS with the structure from Fig.12. There are two types of tensors: on the sites, for the staggered fermionic field (positrons and electrons), and on the links, for the bosonic gauge field. The physical indices at the sites are fermionic, as well as the unoriented bond indices. These indices satisfy the fermionic PEPS rules ftn , namely, every time that two of such lines cross, one needs to include a fermionic swap gate in the TN diagram. Additionally, the physical indices at the links are purely bosonic and correspond to the truncated Hilbert space of the gauge variable for the corresponding link. Finally, bosonic bond indices are introduced with an orientation (arrow), which implement the $U(1)$ gauge symmetry in the tensor components.

In terms of equations, the non-zero components are the following for the tensors at the sites:

[TABLE]

where $\pm$ refers to a positron or an electron, tensor $b^{f}_{\alpha_{q},\beta_{r},\gamma_{s},\delta_{t}}$ corresponds to the free parameters, the first delta implements fermionic $\mathbb{Z}_{2}$ parity symmetry, the second delta implements the gauge $U(1)$ symmetry, and $f=0,1$ is the fermionic occupation number. Similarly, for the tensors at the link the non-zero components are given by

[TABLE]

where $c^{b}_{\alpha_{q},\gamma_{s}}$ are the free variational parameters, $b$ is the bosonic physical index, the first delta implements the fermionic parity symmetry for the bond indices, and the last two deltas take into account $U(1)$ gauge symmetry.

As mentioned above, this ansatz can be optimized in the thermodynamic limit to approximate the ground state of the Hamiltonian in Eq.(50). Such an optimization could be done variationally by using techniques recently introduced optpeps , but it could also be optimized by imaginary time evolution with usual iPEPS algorithms iPEPS . In any case, at every step in the algorithm one must carefully take into account (i) gauge invariance, as we did for the $(1+1)d$ case, but now also (ii) fermionic swaps, coming from the crossings of fermionic wires in the TN diagrams. The optimization of this ansatz by imaginary-time evolution is currently work in progress, and its results will be presented in a future publication.

VI.3 Discussion

Let is now discuss briefly several aspects of QED in $(2+1)d$ that may be relevant for our simulation. First, it is possible to consider compact and non-compact formulations of lattice QED, both with the correct continuum limit noncompact . At the level of implementation, the main difference is the way we write the pure-gauge term, and both formulations on the lattice have slightly different behaviours for the scalings of the chiral condensate and the monopole density. In our case, non-compact QED in $(2+1)d$ could also be simulated with essentially the same scheme that we presented: in fact, one only would need to change the specific form of the plaquette gates. It would be interesting, thus, to benchmark both lattice formulations with our numerical approach. Second, there is also the controversy about the dependence of the chiral condensate with the number of flavours $N_{f}$ Nf . Several works have argued in favour of a critical value $N_{f}^{c}$, so that there is no chiral symmetry for $N_{f}>N_{f}^{c}$, though with no agreement on the actual value of $N_{f}^{c}$ and the type of phase transition. This is a problem that, in principle, could be explored also within our approach by including extra fermionic degrees of freedom for the flavours. However, this may involve larger bond dimensions in the ansatz, making the simulations more costly. Third, the interplay between fermions and monopoles is well known in compact QED in $(2+1)d$ mono , where people have studied the possible survival of the monopole plasma in the presence of dynamical fermions, even in regimes where chiral symmetry is restored. This interesting question can also be addressed in principle by our method, studying the monopole density in terms of the number of flavours and the chiral condensate, up to the restrictions mentioned above. And fourth, there is also the issue of the finite-temperature dependence of the chiral condensate for QED in $(2+1)d$, with the presence of a confinement - deconfinement transition conjectured to be of the BKT type fint . Indeed, it should be possible to address this question with mixed-state versions of infinite-PEPS algorithms, which already exist in the literature for finite-temperature and even for dissipation finTTN . In principle one could extend our variational ansatz to a PEPS-operator (PEPO) with the correct symmetries, to do a finite-temperature simulation.

VII Conclusions

In this paper we have simulated the Schwinger model in the thermodynamic limit on a lattice, by using a gauge-invariant version of the iDMRG algorithm. After discussing the details of the theory and the particulars of one-site and two-site iDMRG, we have approximated the ground state and computed the extrapolation to the continuum of the substracted chiral condensate for several values of the coupling, in good agreement with alternative calculations. These results allowed us to build intuition on how a TN simulation of QED in higher dimensional systems should proceed. In particular, we proposed a gauge-invariant variational ansatz for the ground state of QED in $(2+1)d$ in terms of an infinite-PEPS with bosonic and fermionic degrees of freedom, as well as $U(1)$ gauge-invariant tensors. We discussed also that all the ingredients for such a simulation are in principle available in TN methods: $2d$ fermions, $U(1)$ gauge symmetry, plaquette interactions, and accurate optimization schemes. This simulation in $(2+1)d$ is currently work in progress. We hope that this paper will help to clarify, at least qualitatively, the “big picture” towards TN simulations of lattice gauge theories in higher dimensions, with the target of lattice QCD in $(3+1)d$ on the horizon. We also hope that this paper helps to clarify, specially to the lattice gauge theory community, how one can handle the different ingredients of these field theories in the TN language directly in the thermodynamic limit, in order to simulate elusive regimes in quantum Monte Carlo.

Acknowledgements.

We acknowledge discussions with M.-C. Ba$\tilde{{\rm n}}$uls, K. Cichy, I. Cirac, K. Jansen, E. Rico, M. Rizzi and H. Saito, as well as the Donostia International Physics Center (DIPC), where part of this work was written.