Phase Transitions in $Z_{n}$ Gauge Models: Towards Quantum Simulations of the Schwinger-Weyl QED

TL;DR

This paper investigates $ ext{Z}_n$ lattice gauge theories coupled to fermions in 1+1 dimensions, analyzing phase transitions and their relation to quantum electrodynamics, with implications for quantum simulation of gauge theories.

Contribution

It provides a detailed numerical analysis of $ ext{Z}_n$ models for $n=2$ to 8, demonstrating their phase transition behavior and connection to the Schwinger model in the large-$n$ limit.

Findings

All $ ext{Z}_n$ models exhibit Ising universality class phase transitions.

Phase transition parameters approach those of the Schwinger model as $n$ increases.

Models show spontaneous $CP$ symmetry breaking in the absence of background fields.

Abstract

We study the ground-state properties of a class of $\mathbb{Z}_n$ lattice gauge theories in 1 + 1 dimensions, in which the gauge fields are coupled to spinless fermionic matter. These models, stemming from discrete representations of the Weyl commutator for the $\mathrm{U}(1)$ group, preserve the unitary character of the minimal coupling, and have therefore the property of formally approximating lattice quantum electrodynamics in one spatial dimension in the large-$n$ limit. The numerical study of such approximated theories is important to determine their effectiveness in reproducing the main features and phenomenology of the target theory, in view of implementations of cold-atom quantum simulators of QED. In this paper we study the cases $n = 2 \div 8$ by means of a DMRG code that exactly implements Gauss' law. We perform a careful scaling analysis, and show that, in absence of a background field, all $\mathbb{Z}_n$ models exhibit a phase transition which falls in the Ising universality class, with spontaneous symmetry breaking of the $CP$ symmetry. We then perform the large-$n$ limit and find that the asymptotic values of the critical parameters approach the ones obtained for the known phase transition the zero-charge sector of the massive Schwinger model, which occurs at negative mass.

Tables

Table 1. Table 1: Parameters of the numerical fit of the critical mass as a function of t 𝑡 t , according to the formula m c ​ ( t ) = m 0 + α ​ t + β ​ t subscript 𝑚 𝑐 𝑡 subscript 𝑚 0 𝛼 𝑡 𝛽 𝑡 m_{c}(t)=m_{0}+\alpha\sqrt{t}+\beta\,t , for the various ℤ n subscript ℤ 𝑛 \mathbb{Z}_{n} -models with n = 2 ÷ 8 𝑛 2 8 n=2\div 8 .

$n$	$m_0^{\left(n\right)}$	${\alpha }_n$	${\beta }_n$
2	$0.004\pm 0.001$	$\left(8\pm 5\right)\cdot {10}^{-6}$	$0.0149\pm 0.0003$
3	$-1.0472\pm 0.0001$	$-0.603\pm 0.001$	$-0.02\pm 0.01$
4	$\left(-3\pm 1\right)\cdot {10}^{-7}$	$0.626\pm 0.005$	$0.0290\pm 0.0006$
5	$-0.628\pm 0.001$	$-0.494\pm 0.004$	$-0.015\pm 0.001$
6	$\left(-7.2\pm 0.1\right)\cdot {10}^{-6}$	$0.543\pm 0.005$	$0.026\pm 0.001$
7	$-0.448\pm 0.001$	$-0.435\pm 0.003$	$0.004\pm 0.001$
8	$\left(1.8\pm 0.1\right)\cdot {10}^{-7}$	$0.503\pm 0.004$	$0.022\pm 0.001$

Table 2. Table 2: ℤ 3 subscript ℤ 3 \mathbb{Z}_{3} -model. Critical values m c ​ ( t ) subscript 𝑚 𝑐 𝑡 m_{c}(t) for different values of t 𝑡 t . The error is always 0.025 0.025 0.025 .

$t$	$0.000$	$0.523$	$1.047$	$1.571$	$2.094$	$6.283$	$10.472$	$20.944$	$41.888$	$62.832$	$104.720$
$m_c$	$-1.047$	$-1.340$	$-1.571$	$-1.770$	$-1.948$	$-2.927$	$-3.596$	$-4.767$	$-6.329$	$-7.449$	$-9.115$

Equations

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

$$\{\psi (t,x),\psi (t,x^{\prime })\}$$

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

$$[E(t,x),A(t,x^{\prime })]=\mathrm{i}\delta (x-x^{\prime }).$$

$$H=\int \mathrm{d}x\left\{{\psi }^{†}{\gamma }^0\left[-{\gamma }^1(\mathrm{i}{\partial }_1+gA)+m\right]\psi +\frac{E^2}{2}\right\},$$

$$G(x)={\partial }_{1}E(x)-g{\psi }^{†}(x)\psi (x),$$

$$G(x)\approx 0,$$

$$H=$$

$$+m\sum _x(-1{)}^x{\psi }_x^{†}{\psi }_x+\frac{g^{2}a}{2}\sum _x{E}_{x,x+1}^2,$$

$$V_{x,x+1}(\xi )U_{x,x+1}(\eta )={\mathrm{e}}^{\mathrm{i}\eta \xi }{U}_{x,x+1}(\eta )V_{x,x+1}(\xi ),$$

$$h_t$$

$$+\frac{2m}{g^{2}a}\sum _x(-1{)}^x{\psi }_x^{†}{\psi }_x+\sum _x{E}_{x,x+1}^2,$$

$$G_x\equiv {\psi }_x^{†}{\psi }_x+\frac{1}{2}[(-1{)}^x-1]-(E_{x,x+1}-E_{x-1,x})\approx 0.$$

$$P\,:\,\left\{\begin{array}[]{ll}\psi_{x}\rightarrow\psi_{-x},&\psi_{x}^{\dagger}\rightarrow\psi_{-x}^{\dagger},\\ E_{x,x+1}\rightarrow E_{-(x+1),-x},&U_{x,x+1}\rightarrow U_{-(x+1),-x}^{\dagger},\end{array}\right.$$

$$C\,:\,\left\{\begin{array}[]{ll}\psi_{x}\rightarrow(-1)^{x+1}\psi_{x+1}^{\dagger},&\psi_{x}^{\dagger}\rightarrow(-1)^{x+1}\psi_{x+1},\\ E_{x,x+1}\rightarrow-E_{x+1,x+2},&U_{x,x+1}\rightarrow U_{x+1,x+2}^{\dagger}.\end{array}\right.$$

$$U_{x,x+1}^{\ell }{V}_{x,x+1}^k={\mathrm{e}}^{\mathrm{i}\frac{2\pi }{n}k\ell }{V}_{x,x+1}^k{U}_{x,x+1}^{\ell }\text{mbox}withk,\ell \in {\mathbb{Z}}_n.$$

$$V|v_k⟩={\mathrm{e}}^{-\mathrm{i}2\pi k/n}|v_k⟩.$$

$$U|v_k⟩=|v_{k+1}⟩\text{mbox}fork<n-1,U|v_{n-1}⟩=|v_0⟩.$$

$$e_k=\sqrt{\frac{2\pi }{n}}\left(k-\frac{n-1}{2}+\varphi \right).$$

$$h_c=$$

$$+\frac{2m_c(t)}{\sqrt{t}{g}^2\tilde{a}}\sum _x(-1{)}^x{\psi }_x^{†}{\psi }_x+\sum _x{E}_{x,x+1}^2,$$

$$m_c(t)=\alpha \sqrt{t},$$

$$h_t^{(n)}=$$

$$+\frac{2m}{g_n^{2}a}\sum _x(-1{)}^x{\psi }_x^{†}{\psi }_x+\sum _x{\tilde{E}}_{x,x+1}^2,$$

$$g_n=g\sqrt{2\pi /n}.$$

$$m_c(t)={\alpha }_n\sqrt{t}.$$

$$\alpha =\underset{n\to \infty }{\lim }{\alpha }_n\sqrt{n/2\pi }.$$

$$\Sigma =\frac{1}{N}\sum _x⟨E_{x,x+1}⟩$$

$$\displaystyle{\color[rgb]{1,0,0}h^{(3)}_{t=\frac{2\pi}{3}}}=$$

$$+\frac{3}{2\pi }m\sum _x(-1{)}^x{\psi }_x^{†}{\psi }_x+\sum _x{\tilde{E}}_{x,x+1}^2.$$

$$\varSigma=N^{-\frac{\beta}{\nu}}\lambda\big{(}N^{\frac{1}{\nu}}(m-m_{c})\big{)}$$

$$S_L\left(\frac{L}{2}\right)=\frac{c}{6}{\log }_2(L)+s_0,$$

$$\Delta ={\epsilon }_1(N)-{\epsilon }_0=\frac{\pi v_s{x}_s}{N^2},$$

Full text

I Introduction

The advent of quantum technologies has opened unprecedented perspectives for the simulation of difficult problems, whose solution depends exponentially on input size. In particular, there is the intriguing possibility to explore lattice gauge theories rothe1992lattice ; montvay1997quantum , such as quantum electrodynamics (QED) in 1D. This represents a complex and historical problem, that was investigated by a number of outstanding physicists in the 70’s wilsonlgt ; kogut1975hamiltonian ; susskind1977lattice ; kogut1979introduction .

The possibility of using quantum simulators to tackle these problems has been made possible by recent developments in low-temperature physics and atomic control techniques in optical lattices. Many proposals have been put forward in the literature to use ultra-cold atomic optical lattice systems to simulate many body models bdz ; qsim1 ; qsim2 ; qsim3 ; qsim4 , as well as Abelian and non-Abelian lattice gauge theories simul1 ; simul4 ; simul3 ; simul5 ; simul6 ; simul7 ; simul9 ; simul11 ; zoharreview and in particular quantum electrodynamics in 1D, that appears as a realizable option in the not-too-distant future simul10 ; KCB ; NEFMPP ; simul12 .

The key idea is that cold-atom quantum simulators make possible the implementation of matter fields in presence of artificially designed gauge fields by suitably identifying the gauge degrees of freedom with the internal (for example spin) states of the atom. The first experiments with fermions in presence of such “synthetic” fields have already been proposed and performed, offering very promising perspectives cmr ; ytterbium ; mpi_sun ; fallani ; lcd . Also, a first experiment reproducing 1D QED with few qubits has been reported martinez2016 .

In general, the numerical study of approximated lattice-gauge theories is of paramount importance to determine their effectiveness in reproducing the main features and phenomenology of the target theory. Novel quantum-inspired numerical techniques, such as DMRG- and MPS-based algorithms DMRG1 ; MPS , fully exploit the entanglement of the states that contribute to the dynamics, and are able to reduce the computational cost by suitably tailoring the relevant (effective) subspaces in the Hilbert space. From the theoretical point of view, an approach based on quantum simulators paves the way towards a number of problems that were traditionally difficult to analyze, such as the investigation of possible phase transitions, non-perturbative phenomena and dynamical aspects simul2 ; simul6 ; BCJC ; rico2016 ; montangero2015 ; schwinger_mps ; buyens_prx ; buyens .

One of the main problems one has to face in encoding a lattice gauge-theory model in a cold-atom system is the fact that the number of states of the gauge field, sitting on the links between lattice sites, must be finite. In the literature, this has been accomplished in different ways. In simul2 ; rico2016 , $U(1)$-gauge fields have been replaced by spin variables, which allow for finite dimensional (but non-unitary) representations of the canonical commutation relations. For one-dimensional QED this amounts to consider the so called quantum link model qlm1 ; qlm2 ; qlm3 ; qlm4wiese . Otherwise, one can try to discretize the gauge group, with the advantage of preserving the commutators, and guaranteeing at the same time the unitary character of the minimal coupling structure. This has been done for example in KCB ; buyens , where the gauge group $U(1)$ was replaced by the discrete group $\mathbb{Z}$, which however admits only infinite dimensional representations. In such a case, the Hilbert space describing the gauge degrees of freedom must then be truncated to perform numerical simulations buyens .

In this paper we will follow a different strategy, that arises from the fact NEFMPP that the Weyl group commutation relations arising form the canonical commutators of the U(1) gauge fields admit discrete and finite dimensional representations, thus yielding a discrete and finite implementation of the gauge group commutator at the Hamiltonian level. This is important for two reasons. First, numerical simulations do not introduce any further approximation except for standard finite size effects, that can be dealt with by looking at scaling properties. Second, one can think of experimentally simulating the class of models derived by such strategy, and described in the following in an exact way NEFMPP , by means of trapped cold atoms with synthetic degrees of freedom that satisfy periodic boundary conditions aoc ; cmr .

We will consider $\mathbb{Z}_{n}$ Abelian gauge theories, with varying $n$, in which the gauge field is coupled with a spinless matter field in one spatial dimension. These models have the twofold advantage of a careful control over finite-$n$ effects, and of providing a controlled approximation of lattice U$(1)$ quantum electrodynamics in the large-$n$ limit. We will focus on the ground state properties for different values of $n$ and other parameters of the model. In particular, we will find that, with no background field (corresponding, in the $n\to\infty$ limit, to the $CP$-invariant point of the Schwinger model), the system undergoes a quantum phase transition towards regions of the parameter space in which the parity and charge conjugation symmetries are spontaneously broken. The transition will be characterized in terms of finite-size scaling and critical exponents, which will fall into the Ising universality class. Such quantum phase transition does not survive the introduction of a background field, since in this case the excitations become always massive. We will numerically study the $\mathbb{Z}_{n}$-models for $n=2$ to $n=8$ and extrapolate the results to the large $n$-limit, showing that the corresponding known phase transition of the U$(1)$-model is recovered. The features of this class of transitions and its asymptotic properties will be carefully scrutinized.

The article is organized as follows. In Section II, we introduce the massive Schwinger model in $1+1$ dimensions and discuss the paths to discretization of space and gauge degrees of freedom. Section III includes the definition of $\mathbb{Z}_{n}$ gauge models and the presentation of their general features and scaling properties. In Section V we study in detail the case $n=3$, characterizing its ground state properties and the quantum phase transition at a negative critical mass, in absence of background field. Section VI is devoted to a presentation of the results obtained in all cases $n=2\div 8$, $n\neq 3$, focusing on the different phenomenology of the even and odd cases, while the details of the numerical results of all these cases are given in the Appendix. In Section VII we summarize our results and recover the U$(1)$ model in the limit of large $n$. In section VIII, we comment on a possible implementation of the proposed class of models in a cold atomic platform. We finally draw our conclusions in Section IX.

II Discretization of one-dimensional QED

Quantum Electrodynamics in one spatial dimension is a U$(1)$ gauge theory, describing the interaction of a charged particle (“electron”), represented by a spinor field $\psi(t,x)$, and the electromagnetic field $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, associated to the potential $A_{\mu}$, with $\mu,\nu=0,1$. The classical Lagrangian density is determined by the minimal coupling prescription:

[TABLE]

where $m$ and $g$ are the electron mass and charge, respectively, and $\{\gamma^{\mu},\gamma^{\nu}\}=2\eta^{\mu\nu}$ with $\eta=\mathrm{diag}(1,-1)$. The properties of the theory are strongly characterized by the absence of transverse degrees of freedom: the electron, described by a two-component spinor, is spinless, and the only independent component of the electromagnetic tensor is the electric field $E=F_{01}$. While quantization of the spinor field is determined by the canonical equal-time anticommutators

[TABLE]

a gauge choice is necessary to quantize the electromagnetic potential. In the canonical gauge, the temporal component $A_{0}$ is set to zero, while the spatial component $A:=A_{1}$ is taken as the conjugate variable to $E$:

[TABLE]

This choice, leading to the Hamiltonian

[TABLE]

does not allow one to enforce Gauss’ law $G(x)=0$, with

[TABLE]

as an operator constraint. However, since $[G(x),G(x^{\prime})]=0$ and $[G(x),H]=0$ due to (4), it is possible to select the physical subspace of states $|{\psi}\rangle$ for which $G(x)|{\psi}\rangle=0$, which will be denoted by

[TABLE]

at all space points.

In the following, we will consider two kinds of discretization, towards classical and quantum simulations of the model. The first one is spatial discretization: the continuum model will be replaced by an approximation on a linear lattice of points with spacing $a$, making the continuous space variable $x\in\mathbb{R}$ discrete: $x\in\mathbb{Z}$. The second one is the approximation of the gauge group U$(1)$ with a finite group, which is essential if one wants to work with a finite number of local degrees of freedom in the gauge variables. This can be done essentially in two ways, based on the generalization of the commutation relation $[E,U]=\eta U$, where $E$ and $A$ are two conjugated operators ($[E,A]=\mathrm{i}$), $U=\mathrm{e}^{-\mathrm{i}\eta A}$ is the gauge comparator and $\eta\in\mathbb{R}$ is a constant with the same dimensions as $E$. One option is to focus on the preservation of the above commutator. This is the approach taken, for example, in quantum link models qlm1 ; qlm2 ; qlm3 ; qlm4wiese , in which the operators $E$ and $U$ are replaced with spin variables. Another option is to require that the group commutator $\mathrm{e}^{\mathrm{i}\xi E}\mathrm{e}^{-\mathrm{i}\eta A}=\mathrm{e}^{\mathrm{i}\eta\xi}\mathrm{e}^{-\mathrm{i}\eta A}\mathrm{e}^{\mathrm{i}\xi E}$, which is equivalent to the previous one in the U$(1)$ case, be satisfied by unitary operators for discrete values of $\eta$ and $\xi$ NEFMPP . We will follow the latter strategy, that entails the reduction of gauge invariance to a finite group $\mathbb{Z}_{n}$. A pictorial representation of the gauge degrees of freedom is shown in Fig. 1.

As emphasized in the Introduction, we thus obtain an exact finite implementation of the gauge group commutator at the Hamiltonian level. What is yet to be understood is the proper scaling of the dynamical variables in order that the $\mathbb{Z}_{n}$ model correctly reproduce a lattice U$(1)$ model, with a continuous electric field. Clearly, $n$ must go to $+\infty$. However, proper scaling will be necessary for $E=E^{(n)}$. As we shall see, this scaling is dictated by Schwinger himself schwinger2001quantum .

III Lattice $\mathbb{Z}_{n}$-QED model

The Schwinger model (5) described in the previous section can be discretized on a one-dimensional lattice. For convenience, we shall first redefine, the vector potential $A\to A/g$ and the electric field $E\to gE$. With this transformation, that leaves the commutation relation (4) unchanged, the charge is absorbed in the minimal coupling, to reappear in the energy density of the free electric field. Correspondingly, we redefine the parameters in the one parameters groups $U(\eta),V(\xi)$ by $\eta\rightarrow\eta/g,\xi\rightarrow g\xi$.

The lattice Hamiltonian reads montvay1997quantum ; qlm4wiese ; melnikov

[TABLE]

with $x$ labelling the sites of a one-dimensional lattice of spacing $a$. Here:

i) fermionic matter is represented by one-component spinor creation/annihilation operators $\psi_{x}^{\dagger},\psi_{x}$, defined on each site $x$, so that $\sum_{x}\psi_{x}^{\dagger}\psi_{x}$ represents the total number of fermions in the system;

ii) a staggered mass $(-1)^{x}m$ is introduced, so that the positive- and negative-mass components of the Dirac spinor live respectively at even and odd lattice sites, avoiding in this way the fermion-doubling problem kogut1975hamiltonian ; susskind1977lattice ; rothe1992lattice ; montvay1997quantum ;

iii) gauge fields are defined on the links $(x,x+1)$ of the lattice through the pair of variables $E_{x,x+1}$ (electric field) and $A_{x,x+1}$ (vector potential) or equivalently the comparators $U_{x,x+1}(\eta)=\mathrm{e}^{-\mathrm{i}\eta A_{x,x+1}}$ and $V_{x,x+1}(\xi)=\mathrm{e}^{+\mathrm{i}\xi E_{x,x+1}}$, commuting at different sites, and satisfying the group canonical commutation relations:

[TABLE]

for $\xi,\eta\in\mathbb{R}$, which are equivalent exponentiated versions of the (algebra) canonical commutation relations $[E_{x,x+1},A_{x,x+1}]=\mathrm{i}$.

Notice that in (III), the coupling constants $m$ and $g$ are the same as in the continuum, both with the dimensions of energy, while all the field operators are dimensionless.

In analogy to what is done in fermionic lattice models that are usually used to describe condensed matter systems and for reasons that will become clear in the next section, in the following we will actually use a slightly modified Hamiltonian by introducing an additional dimensionless parameters $t$ in front of the first, kinetic, term so to write the dimensionless Hamiltonian

[TABLE]

proportional to the Hamiltonian density. We can also see from this expression that the coefficient $g^{2}a/2$ fixes the scale of the mass, while the coefficient $g^{2}a^{2}$ the one of the parameter $t$. Thus, the numerical simulations described in Sections V and VI will be performed by setting $g^{2}a/2=1$ and $g^{2}a^{2}=1$. The standard lattice Schwinger model is then recovered for $t=1$.

In addition, the theory must respect Gauss’ law, encoding the U$(1)$ gauge symmetry of the model, that reads

[TABLE]

The Hamiltonian (III) is invariant under $C$ and $P$ symmetries that for staggered fermions read (assuming an infinite lattice or setting $-x\equiv 2L-x$ for a finite lattice with an even number $N=2L$ of sites, labeled from [math] to $2L-1$):

[TABLE]

[TABLE]

The Hamiltonian (III) can be approximated via a discrete Abelian $\mathbb{Z}_{n}$-gauge model NEFMPP , that can be obtained from a finite dimensional representation of the two-parameter projective unitary Weyl weyl1950theory group $\left\{\mathrm{e}^{\mathrm{i}(\xi E_{x,x+1}-\eta A_{x,x+1})}\right\}_{\xi,\eta\in\mathbb{R}}$.

For the two particular cases $(\xi,\eta)=(0,\sqrt{2\pi/n})$ and $(\xi,\eta)=(\sqrt{2\pi/n},0)$ one gets the two operators $U_{x,x+1}=\mathrm{e}^{-\mathrm{i}\sqrt{\frac{2\pi}{n}}A_{x,x+1}}$ and $V_{x,x+1}=\mathrm{e}^{\mathrm{i}\sqrt{\frac{2\pi}{n}}E_{x,x+1}}$, that satisfy the commutation relations

[TABLE]

which is a discrete $\mathbb{Z}_{n}$ version of (9). This representation can be implemented by considering an $n$-dimensional Hilbert space ${\cal H}_{n}$ defined on each link, and choosing an orthonormal basis $\{|v_{k}\rangle\}_{0\leq k\leq n-1}$. Dropping the link index, we consider the diagonal operator $V$ acting as

[TABLE]

The operator $U$ is instead defined as that operator that performs a cyclic permutation of the basis states:

[TABLE]

Some simple algebra shows that these operators do indeed satisfy the Schwinger-Weyl commutation relations (14). Let us remark that this representation exactly implements the unitarity of both operators.

Thus the dynamics of the $\mathbb{Z}_{n}$-model is determined by the Hamiltonian (III), where the discrete version of the electric field $E_{x,x+1}$ is given by the Hermitian operator that is diagonal in the $\{|v_{k}\rangle\}$ basis, with eigenvalues

[TABLE]

In all cases, the eigenvalues of the electric field are symmetric around zero, with a maximum value $E_{\mathrm{max}}=\sqrt{2\pi/n}(n-1+\phi)/2$. We notice that, for $\phi=0$, it is possible to have zero electric field only if $n$ is odd. A value $\phi\neq 0$ corresponds to adding a background field that can be obtained by placing charges at the boundaries of the chain, thus yielding different charge sectors, that are known to be super-selected. It is indeed known that this model displays $\theta$-vacua coleman2 , which can be related to the axial anomaly manton via the spectral flow of the Hamiltonian operator (when imposing periodic boundary conditions) or to unusual twisted boundary conditions for the fermionic field wipf . In this case the $P$ and $C$ symmetries are explicitly broken. In the case of even $n$, the minimum eigenvalues (17) are doubly degenerate for $\phi=0$. As a consequence, in the strong-coupling limit, in which the $U$-dependent terms in (III) are neglected, the energetic cost of creating a fermion-antifermion pair from the vacuum vanishes: this feature is typical of theories with $\theta=\pi$ coleman2 .

IV Scaling properties of the Hamiltonian

Before starting to numerically investigate the Hamiltonian, some comments are in order to establish the correctness of Eq. (III) to suitability represent a quantum simulator for 1-dimensional QED. Being in particular interested in its critical properties, our analysis needs to contain a careful check of the scaling properties of the discretized Hamiltonian as we change the different parameters that appear in it, including the $a\rightarrow 0$ limit (continuum limit), the $N\rightarrow\infty$ limit (infinite volume limit), the $n\rightarrow\infty$ limit ($U(1)$-limit). Close to a critical point, at which physical constants and observables are functionally related by universal laws, it is very hard to control these different cases independently, both from an analytical and a numerical point of view. However, we can resort to well-known techniques based on a finite-size scaling analysis guided by universal scaling properties. We have chosen to perform this study in two steps: first, we consider a particular $\mathbb{Z}_{n}$-model, by keeping $n$ fixed, and perform a finite-size scaling in the dimension of the spatial lattice; second, we let $n$ increase and analyze the large-$n$ limit.

To this end, we first notice that the parameter $t$ can be used to understand critical properties of the model. Let us suppose that the system undergoes a phase transition for a critical value of the mass, $m_{c}(t)$, which may depend on $t$. At this particular point the dimensonless Hamiltonian (III) should be scale invariant. The coefficient $t$ can be absorbed in a re-scaling of the lattice spacing, $a\rightarrow\tilde{a}=a/\sqrt{t}$, and:

[TABLE]

In a mean-field approach, in which possible anomalous dimensions of the field are neglected, the coefficient in front of the second addend must be independent of $t$. In other words, the critical value of the mass scales like

[TABLE]

where $\alpha\equiv m_{c}(t=1)$. We will examine accurately how the critical value of the mass depends on $t$ in the numerical simulations of the next sections, where we will see that its behaviour does not significantly deviates from the one predicted here. Therefore we will obtain the continuum limit critical mass $m_{c}$ by setting: $m_{c}=m_{c}(t=1)=\alpha$.

Incidentally, let us remark that the limit $t\rightarrow\infty$ is not equivalent to the limit $a\rightarrow 0$, since the coefficient $t$ weights the kinetic term differently with respect to the mass and the electric energy terms; in particular, in the $t=0$ case we recover a classical limit which can be exactly solved, while in the large $t\rightarrow\infty$ limit we deal with a pure kinetic Hamiltonian which cannot display any phase transition.

Second, we want to study the large-$n$ limit. It is important to notice that, as shown in NEFMPP , the scaling of the eigenvalues of the electric field with $n$ as given in Eq. (17) is fixed by requiring that the U$(1)$-limit is recovered when $n\rightarrow+\infty$. Also, recalling that two consecutive values of the electric field differ by $\sqrt{2\pi/n}$, it is convenient to collect such a factor and work with the dimensionless Hamiltonian

[TABLE]

where now $\tilde{E}_{x,x+1}$ has eigenvalues $(k-(n-1)/2+\phi)$ with unit spacing (hence, independent of $n$), and

[TABLE]

In the same spirit as before, we can conclude that now (with a slight abuse of notation)

[TABLE]

Comparing the Hamiltonian density (IV) to (IV) in the limit $n\to\infty$, we can conclude that the asymptotic value of $\alpha_{n}g/g_{n}$ must approach the coefficient $\alpha$ appearing in Eq. (19), namely

[TABLE]

In the following, to perform numerical calculations, we will consider the Hamiltonian (IV) defined on a lattice of size $N=2L$ with open boundary conditions. We will work in the sector with one fermion for each “physical site”, i.e. with $N_{\textrm{part}}=N/2=L$ particles. Also, as explained above, we will set $g^{2}a/2=1$ (which sets the units of energy) and $g^{2}a^{2}=1$ (which sets the units of $t$).

We will first present the $\mathbb{Z}_{3}$-model, in order to illustrate all the details of our treatment. We will then discuss the general $\mathbb{Z}_{n}$-model, for both odd and even $n$. As we will see, these two cases need to be considered separately.

V Lattice $\mathbb{Z}_{3}$-QED model

V.1 Hilbert space and gauge-invariant subspace

As mentioned in Sec. III, in the Schwinger model each “physical fermion” is represented by a pair of staggered fermions sitting in nearby sites, with even/odd sites occupied by positive/negative mass particles. Thus the vacuum state (Dirac sea) is obtained by leaving the even sites empty and occupying the odd ones. The presence/absence of a fermion in an even/odd site is interpreted as the presence of a quark/anti-quark, while a meson is a configuration made up of a quark and an anti-quark. This is shown in Fig. 2(a). On each link $(x,x+1)$, the electric field can only assume one of the three values $E=\sqrt{2\pi/3}(k-1)$, with $k\in\mathbb{Z}_{3}=\{0,1,2\}$, which will be represented as an arrow pointing left, an un-oriented segment, and an arrow pointing left, respectively, as shown in Fig. 2(b).

Thus the total Hilbert space associated with an even site, together with its two adjacent links, is $2\times 3\times 3=18$ dimensional. But Gauss’ law forces the physical states to belong to an invariant subspace which is constructed out of those states for which the value of the electric field $k_{r}$ at the right link is either $k_{r}=k_{l}$, if the site is empty, or $k_{r}=k_{l}+1(\mathrm{mod}\,3)$ if the site between links is occupied by a (positive mass) fermion. The situation is similar for odd sites, for which Gauss’ law constrains physical states to have either $k_{r}=k_{l}-1(\mathrm{mod}\,3)$, if the site between links is empty, or $k_{r}=k_{l}$ if the site is occupied by a (negative mass) fermion. This is displayed in Fig. 3(a). Notice that we have $2\times 3=6$ independent configurations for each site. The gauge invariant states of a “physical site” are obtained by gluing together an even and an odd site that share a common value for the electric field in the link between them, obtaining $2\times 6=12$ possible configurations, as shown in Fig. 3(b).

It is easy to see that, for a chain with $N$ sites (with open boundary conditions), the dimension of the gauge-invariant subspace is $2^{N}\times 3$. Some notable examples of gauge-invariant states in a chain are shown in Fig. 4.

The ground state of the Hamiltonian (III) will be given by the completely filled Dirac sea (see Fig. 5(a)) for large positive $m$, while for large negative $m$ the system will tend to choose between the two states shown in Fig. 5(b), where mesons/antimesons have formed. Notice that the Dirac sea is invariant under both parity and charge conjugation, while $P$ and $C$ map the mesonic and antimesonic states into each other.

These two cases are clearly distinguished by the mean value of the electric field operator

[TABLE]

that we will use as a kind of order parameter, since it vanishes for the Dirac sea and takes the values $\pm\pi/3$ for the mesonic/antimesonic states. An abrupt change of this quantity signals the existence of a possible phase transition, whose existence has to be confirmed by looking at the appropriate scaling of observables and thermodynamical quantities.

V.2 Numerical investigation of the critical point with no background field

We perform numerical calculations by means of a finite-size DMRG code DMRG2 in which gauge invariance is exactly enforced. This is obtained by using a unit cell given by a pair of nearby (even and odd) sites, whose local Hilbert space is the span of the 12 gauge-invariant states described in Fig. 3. This is implemented at each step of the algorithm, with a twofold advantage: decreasing the computing time by working in a restricted space, and avoiding transitions out of the gauge-invariant subspace. We work with up to $N=80$ sites ($L=40$ pairs), while keeping $1000$ DMRG states at most. These values are large enough to ensure stability of our findings and small errors.

We first numerically study the Hamiltonian (IV) at the $CP$-invariant point, i.e. in absence of a background field: $\tilde{E}_{x,x+1}\in\{-1,0,+1\}$. We start by choosing $t=2\pi/3$ so as to work with the operator

[TABLE]

Notice that, here and in the following sections, mass is expressed in units of $g^{2}a/2$. The behavior of the observable $\varSigma$ as a function of $m$ is displayed in Fig. 6(a) for different system sizes, ranging from $L=12$ to $L=40$. We see that, as expected, $\Sigma$ essentially vanishes at large positive $m$ and tends to the value $\sqrt{2\pi/3}/2\simeq 0.724$ for large negative $m$.

In Fig. 6(b) we zoom on the central region, showing a steeper transition as the system size increases. This strongly suggests that we are in presence of a phase transition, at a critical value of the mass which corresponds to the point where all curves intersect. We can estimate this value if we make a hypothesis about the nature of the phase transition: indeed, if we know the critical exponents, we can calculate $m_{c}$ by using the fact that $\varSigma$ should scale with the system size $N$ according to the finite-size scaling formula Cardy

[TABLE]

where $\lambda$ is a universal function. By taking into account suggestions from the continuum limit coleman2 and the symmetries of the model, we anticipate that the phase transition is of the Ising-type, so that $\beta=1/8$ and $\nu=1$. A fit of the data yields then $m_{c}=-1.948\pm 0.025$, where the error has been estimated as the semi-interval between the numerical points. We now have to look at the numerical curves given by $N^{\frac{\beta}{\nu}}\varSigma$ versus $N^{\frac{1}{\nu}}(m-m_{c})$, for different $N$, which should all collapse onto the same universal curve $\lambda(x)$. This behaviour is clearly seen in Fig. 7.

Once we have an estimate for the critical mass, we can confirm the validity of our hypothesis by calculating other quantities. Figure 8 displays the entanglement entropy of a subsystem of size $L/2$ at the critical point, which—according to conformal field theory cardy_calabrese —should scale logarithmically with the system size according to the law

[TABLE]

where $s_{0}$ is a constant (which can depend on the boundary conditions and other details of the model) while $c$ is the central charge. The fit yields $c=0.51\pm 0.01$, in perfect agreement with the central charge of the Ising model, $c=1/2$.

Additional information can be obtained by looking at the scaling of excited states with the size of the system: Figures 9(a) and 9(b) display the behaviour of the gaps $\Delta$ and $\Gamma$ of the first two excited states, which—again according to conformal field theory henkel —should obey

[TABLE]

$\varepsilon_{0}$ being the ground state energy density, $v_{s}$ a speed, and $x_{s}$ the surface critical exponent, which is equal to 2 for the Ising model with open boundary conditions. We numerically find

[TABLE]

Plugging this result back into Eq. (28) we can also estimate the speed $v_{s}$, obtaining

[TABLE]

(a number very close to $\pi/2$). We remark that surface exponents are found for states that can be obtained from the ground state by changing from periodic to anti-periodic boundary conditions. The system here shows spontaneous symmetry breaking to one of the two degenerate polarized states and thus the first excited state can be represented as a kink-like solution which interpolates between these two degenerate polarized minima, in agreement with what is found in the continuum coleman2 . These results are also fully compatible with recent buyens ; shimizu_kuramashi and less recent byrnes2002 numerical results. Interestingly, similar conclusions hold also at finite temperature sachs_wipf .

For $t=0$, the system undergoes a (first order) phase transition between the Dirac sea, depicted in Fig. 5(a), with an energy per pair $E_{\textrm{Dirac}}/L=-m$, and the mesonic state in Fig. 5(b), with energy per pair $E_{\textrm{meson}}/L=m+2\pi/3$. The critical value of the mass $m_{0}^{(3)}=-\pi/3\simeq-1.047$ is obtained at $E_{\mathrm{meson}}=E_{\mathrm{Dirac}}$.

In order to test the validity of Eq. (23), we have repeated this procedure for several values of $t$, checking that the Ising transition is always present and calculating numerically $m_{c}(t)$. Our numerical findings for $m_{c}(t)$ as a function of $t$, as well as other useful information, are given in Appendix A.

A numerical fit of the form

[TABLE]

yields the values

[TABLE]

Let us notice that, as expected from the predicted behaviour (23), the coefficient of the linear term is much smaller than the one of the square-root term, thus yielding a negligible correction, at least for not too large values of $t$. This will also be apparent in Fig.  15 (green points and green fitting curve).

V.3 Numerical investigation of the critical point in presence of a background field

It is known coleman2 that the Schwinger model should exhibit a phase transition only at the $CP$-invariant point. In order to check if this is the case also in our model, we have scrutinized the effects of a constant background field. We present here just one representative example, by considering an electric field

[TABLE]

in the Hamiltonian (V.2).

As shown in Fig. 10(a), the observable $\varSigma$ still shows a very sharp transition between a negative and a positive value. But we are now in presence of a cross-over, rather than a phase transition, as it can be inferred by performing a scaling analysis. Indeed, the function $\lambda$ in Eq. (26) changes for different system size $N$ and does not have a universal character, as one can infer from Fig. 10(b). Also, the entanglement entropy $S_{L}(l)$ does not scale with the size $l$ of the interval, as predicted by conformal field theory cardy_calabrese , but is rather constant, except for some small edge effects, also at the crossing point $m^{*}=-0.325$ (see Fig. 11(a)). These results suggest that, in presence of a background field, the gap never closes, as the numerics confirms (see Fig. 11(b), red dots).

Finally, we also checked the case of a background electric field which is halfway between two integer values:

[TABLE]

Also in this case the model is gapped for any value of the mass, as shown in Fig. 11(b) (green dots). At a first sight, this situations looks very similar to the case of even $n$ with no background field, when the possible spectrum of the electric field does not include zero, being still invariant under a sign change. However, the two cases are very different, since, as we will discuss in the next section, the $\mathbb{Z}_{n}$-model with even $n$ and no background field, which is $CP$-invariant, still exhibits a phase transition.

VI Lattice $\mathbb{Z}_{n}$-QED model for other values of $n$

The analysis performed in the previous section for the $\mathbb{Z}_{3}$-model can be repeated for all values of $n$. In the following we will consider only the case with no background field. The dimension of the gauge-invariant Hilbert subspace for a chain with $N$ sites is $2^{N}\times n$, increasing only linearly with $n$, since the electric field can now take the $n$ values $-\sqrt{\frac{2\pi}{n}}\frac{n-1}{2},\cdots,+\sqrt{\frac{2\pi}{n}}\frac{n-1}{2}$. Similarly to what was done for the $\mathbb{Z}_{3}$ case, we rescale the electric field as: $E_{x,x+1}=\sqrt{2\pi/n}\tilde{E}_{x,x+1}$ and study the Hamiltonian (IV), with $g^{2}a^{2}=g^{2}a/2=1$. One must consider odd and even $n$ separately.

VI.1 Odd $n$

As for the case $n=3$ presented in the previous section, if $n$ is odd we anticipate a phase transition from a phase where the ground state is the Dirac sea for large positive $m$, to a phase in which the ground state is a meson/antimeson state for large negative $m$. At $t=0$, there is a first-order phase transition between these two states, at a critical mass $m_{0}^{(n)}$ that can be easily found by comparing the energy of these two states, given respectivley by $E_{\mathrm{Dirac}}/L=-m$ and $E_{\mathrm{meson}}/L=m+2\pi/n$, thus yielding the critical value

[TABLE]

Fot $t\neq 0$ we must resort to our DMRG code and perform an analysis identical to the one presented for the $\mathbb{Z}_{3}$-model, assuming again that the phase transition falls in the Ising universality class. As an example, in Fig. 12 we show the behaviour of the function $\lambda$ of Eq. (26) for different system size in the $\mathbb{Z}_{5}$-model and for $t=2\pi/5$, proving its universality in this case as well. We have performed an exhaustive analysis of the $\mathbb{Z}_{5}$- and $\mathbb{Z}_{7}$-models, obtaining the value of $m_{c}(t)$ as function of $t$ in both cases, as summarized in the Appendix and in Fig. 15.

We can now fit these data with the formula (32) to get an estimate of the coefficients $\alpha_{n}$, $\beta_{n}$ with $n=5$ and $7$. The numerical results are summarized in Table 1 and show an excellent agreement with the theoretically predicted value $m_{0}^{(n)}=-\pi/n$. Also, as for the $n=3$ case, the coefficient of the linear term is much smaller than the one of the square-root term.

VI.2 Even $n$

Even-$n$ models are different from odd-$n$ ones since the electric field cannot take the value zero, still being $CP$-invariant. This means that we are working in a different super-selection sector corresponding to a different total charge at the boundary.

Let us first consider the case $n=2$. This is a very small (in fact, the smallest non-trivial) value, and one may expect some peculiarities, due to the fact that the scaling arguments in Eqs. (IV)-(19) do not apply. Actually, since the electric field Hamiltonian becomes trivial for $\phi=0$, the presence of a phase transition only depends on the ratio of the coefficients of the hopping term and the mass term. Hence, the critical mass would be linear in $t$. The gauge-invariant Hilbert subspace for a pair of sites is $8$-dimensional and a basis is shown in Fig. 13(a). The electric field can assume the two values $-\sqrt{\pi}/2,+\sqrt{\pi}/2$. The transition is from a phase where the ground state is the uniformly polarized vacuum, for large positive $m$, to a ground state in which the electric field has alternating signs on links, for large negative $m$. These states are shown in Fig. 13(b).

For $t=0$ the energy per pair of these two states can be calculated exactly:

[TABLE]

while

[TABLE]

Thus a first order phase transitions occurs at $m_{0}^{(2)}=0$. For $t\neq 0$ we look for the phase transition by numerically calculating the observable $\varSigma$ as function of $m$ and performing a finite-size scaling of the universal function that describes the order parameter. Fig. 14(a) displays $\varSigma$ for $t=2\pi/2=\pi$, from which we can calculate the critical value $m_{c}=0.016\pm 0.025$, while Fig. 14(b) shows the corresponding universal function.

As for the odd-$n$ case, we can numerically evaluate the critical value of the mass for different values of $t$ and get the fit of the function $m_{c}(t)$ according to Eq. (32)

[TABLE]

From these values, we can see that both coefficients are very small, the dominant one being associated with the linear term. This is indeed one issue of the $\mathbb{Z}_{2}$-model that, as we will presently see, is not shared by higher $n$-models.

Taking into account this peculiarity and considering that we want to perform a large-$n$ limit by using at least three different values of $n$ both in the even and odd case, we have performed a similar analysis for the $\mathbb{Z}_{4}$, $\mathbb{Z}_{6}$ and $\mathbb{Z}_{8}$-models, whose results are summarized in the Appendix. By fitting these data, we obtain the coefficients $m_{0}^{(n)}$, $\alpha_{n}$ and $\beta_{n}$ of Eq. (32) as given in Table 1. As for the odd case, we find an excellent agreement of the numerical value for $m_{0}^{(n)}$ with the theoretically predicted value, $m_{0}=0$, and confirm that the dominant term is the one containing $\sqrt{t}$, as expected from Eq. (19).

VII Large-$n$ limit

The values of $m_{c}(t)$ for $n=2\div 8$, are summarized in Table 1 and plotted in Fig. 15. Leaving aside the peculiar $n=2$ case, these data clearly show that the critical mass actually shows a square-root dependence on $t$

[TABLE]

where, for any $n$, the critical mass at $t=0$ can be calculated analytically (with $g^{2}a=2$) according to the formula

[TABLE]

and vanishes in the large-$n$ limit. The coefficients $\alpha_{n}$ can be read from the third column of Table 1. As we can see from Fig. 16, they obey the scaling

[TABLE]

and $b=0$ within numerical error in both cases. Except for the different sign, which is due to the fact that the even $n$-models do not admit a zero electric field, these two values are the same.

Thus, combining (23) and (46), we conclude that the continuous U$(1)$ theory exhibits a phase transition at the critical mass ($t=1$)

[TABLE]

with the sign depending on the charge sector. This value is in very good agreement with the estimates $m_{c}/g=0.33(2)$, obtained by using a lattice Hamiltonian approach HKCM , and $m_{c}/g=0.3335(2)$, obtained by studying the truncated $\mathbb{Z}$-model (at most at the first five loop levels) byrnes2002 .

VIII Cold-atom simulator

The implementation of the $\mathbb{Z}_{n}$ models is complicated by the presence of the correlated hopping terms, related to elementary processes in which the hopping of a fermion to a nearest-neighboring site is always associated to an action on the link between the sites, which amounts at increasing the electric field in the case of hopping to the left, and decreasing it if the fermion hops to the right. The accuracy of correlated hopping terms is vital for any cold-atomic simulator of the described theories, since it guarantees that, once the system starts in the physical subspace, in which Gauss’ law is satisfied, it will not leave this subspace during the evolution. However, in a quantum simulation, in which matter and gauge fields are encoded in the external and (possibly) internal degrees of freedom of cold atoms, Gauss’ law does not emerge as a natural property. Implementation of correlated hopping and enforcement of Gauss’ law are therefore still open problems. We can identify two possible ways to simulate the gauge variables with cold atoms:

• •

Gauge variables can be encoded in the internal degrees of freedom of single atoms trapped at intermediate positions between each couple of adjacent sites. Hopping of a fermion induces transition towards different states according to the hopping direction. This realization require a fine tuning of atomic transitions, as care must be taken in ensuring that all the allowed transition amplitudes between states with given fermion occupation numbers and electric field are equal. Moreover, the pure-gauge term requires that the energy levels of the intermediate atoms at $a\to\infty$ are quadratically spaced.

• •

Gauge variables can be encoded in an external, transverse degree of freedom. A possible interesting implementation arises from the possibility of trapping cold atoms in circular lattices, obtained by interaction with Laguerre-Gauss laser modes aoc . The scheme is represented in Figure 17, where the red spots represent the bottoms of potential wells in which the fermions are trapped, while the blue ones host one particle per link (statistics is immaterial), which can hop through neighboring sites of each circle (identified with eigenstates of the electric field), but cannot hop towards other links due to a large energy barrier. The equal amplitude of hopping between sites on the circle arises from a natural circular symmetry, and the pure-gauge term can be implemented by adding an external potential that properly varies along one of the transverse directions.

In both cases, the Gauss law could be implemented either by tailoring the transition amplitudes in order to enhance correlated hopping and suppress the forbidden terms, or by adding an energy or noise penalty to the states that violate Gauss’ law simul2 ; simul11 ; zohar2011 . In the latter case, the desired interaction Hamiltonian can emerge as a higher-order effective dynamics NEFMPP .

IX Conclusions

We have investigated discrete $\mathbb{Z}_{n}$ models, that approximate QED in one dimension (Schwinger model). In these models the electric field can take a finite number of values, and one important common feature is the preservation of the unitarity of the comparator. Thus, we have put the large-$n$ limit on a firm mathematical ground, adding novel rigorous results in the field of quantum simulations of gauge fields, that may soon find experimental verifications in cold-atomic systems.

In particular, we have unveiled the presence of phase transitions, whose features depend in an interesting way on whether $n$ is even or odd. Although the details of these transitions depend on $n$, their universality class, as well as some of their main features, are $n$-independent, so that by looking at the large $n$ limit, in which $\mathbb{Z}_{n}\to U(1)$, one can establish the presence of a phase transition for one-dimensional lattice QED, and extract crucial information.

A possible implementation of $\mathbb{Z}_{n}$ models on a cold-atom simulator, discussed in NEFMPP and reviewed in Section VIII, relies on the identification of the discrete values taken by the electric field with some suitable additional degrees of freedom of the simulator aoc . A realization appears realistic and would be important to elucidate some important features of one-dimensional QED. Clearly, $1+1$-dimensional models have to be considered as toy-models with respect to the more realistic $3+1$-dimensional ones, but the possibility of using quantum simulators for the investigation of collective and non-perturbative features of gauge theories would enable us to shed new light on old problems, and provide new insights on crucial but still unsolved questions.

The dynamics of lattice U$(1)$ gauge theories has been recently experimentally investigated in a few-qubit trapped-ion quantum computer martinez2016 . Therefore, among the computational perspectives, it is worth mentioning the study of real time dynamical phenomena, such as string breaking marcos2013 ; montangero2015 ; buyens_prx , and the time evolution of localized excitations in the different phases identified in this work.

Acknowledgments

The authors are grateful to Fabio Ortolani for his precious help with the DMRG code. EE, PF, GM and SP are partially supported by INFN through the project “QUANTUM”. FVP is supported by INFN through the project “PICS”. PF is partially supported by the Italian National Group of Mathematical Physics (GNFM-INdAM).

Appendix A Additional information on the phase transition

We give here additional details on the phase transition, as well as numerical figures, for different $\mathbb{Z}_{n}$-models.

Let us start with $n=3$. For $t=0$, the transition is sharp for every system size, as can be seen in Fig. 18. As explained at the end of Sec. V.2, the system undergoes here a first order phase transition between the Dirac sea and the mesonic state shown in Figs.  5(a) and 5(b). Our numerical findings for $m_{c}(t)$ as a function of $t$ are summarized in Table 8.

Our numerical findings for $m_{c}(t)$ for the $\mathbb{Z}_{5}$ and $\mathbb{Z}_{7}$-models are given in Tables 8 and 8, respectively. Those for the $\mathbb{Z}_{2}$, $\mathbb{Z}_{4}$, $\mathbb{Z}_{6}$ and $\mathbb{Z}_{8}$-models are reported in Tables 8, 8, 8, and 8, respectively.

All the values given in the Tables are plotted in Fig. 15, to yield the fit in Eq. (44).