Duality transformations and the entanglement entropy of gauge theories

The study of entanglement in gauge theories is expected to provide insights into many fundamental phenomena, including confinement. However, calculations of quantities related to entanglement in gauge theories are limited by ambiguities that stem from the non-factorizability of the Hilbert space. In this work we study lattice gauge theories that admit a dual description in terms of spin models, for which the replica trick and R\'enyi entropies are well defined. In the first part of this work, we explicitly perform the duality transformation in a replica geometry, deriving the structure of a replica space for a gauge theory. Then, in the second part, we calculate, by means of Monte Carlo simulations, the entropic c-function of the $\Z_2$ gauge theory in three spacetime dimensions, exploiting its dual description in terms of the three-dimensional Ising model.

In general, defining entanglement measures for multipartite and mixed states involves some difficulties [31], but in the case of bipartite entanglement for pure states there exist well established measures, among which the entanglement entropy is one of the most studied.Entanglement entropy can be thought of as the amount of entropy associated with the degrees of freedom in a spatial region A of a quantum system, and is usually defined in terms of the von Neumann entropy S = − Tr(ρ A ln ρ A ) of the reduced density matrix ρ A , obtained by tracing the density matrix ρ of the system over the degrees of freedom of the complement of A (that we denote as B): ρ A = Tr B ρ.It should be noted that this definition involves non-trivial subtleties [32]: in particular, it relies on the hypothesis that the Hilbert space of the theory factorizes into a direct product of the Hilbert spaces defined for the subsystem A and for its complement, H = H A ⊗H B , which, strictly speaking, is not satisfied for a quantum field theory in continuum, nor for gauge theories-not even when they are regularized on a lattice.When the factorization of the Hilbert space holds, the entanglement entropy can be expressed as the n → 1 limit of the Rényi entropy of order n associated to the subsystem A: In quantum field theories the entanglement entropy is affected by ultraviolet divergences: for highly excited states the divergence is proportional to the spatial volume of the system [33], whereas ground states exhibit the so-called "area law" [34], namely the ultraviolet divergence is proportional to the area of the entangling surface ∂A separating A and B. Therefore, a proper regularization is required, in order to identify the subleading contributions, which encode important physical information about the "resolution scale" at which the degrees of freedom of the system are probed [18,19,[35][36][37][38][39], the geometry of the system [40][41][42][43][44] and the nature of the state itself [45].For a system defined in D = d + 1 spacetime dimensions, it is then convenient to consider a geometry in which the entangling surface does not depend on the size of the subsystem A; a common choice is the slab geometry, in which A is maximally extended in all space-like directions except for the one separating A and B, where it has a finite size l.One can then define the ultraviolet-finite entropic (Rényi) c-function as [46][47][48] C where |∂A| denotes the area of the entangling surface.It has been proven in a large number of cases that the entropic c-function is monotonically decreasing along the renormalization group flow [36,[48][49][50], providing a measure of the effective number of degrees of freedom of a theory [51].Calculations of such highly non-local quantities are typically challenging, both analytically and numerically.A tool which is commonly used is the replica trick [2,52], that allows one to express the trace of a power of the reduced density matrix appearing in eq.(1.1) in terms of ratios of partition functions, where Z n denotes the product of the partition functions of n independent copies of the system, while Z n is the partition function of the theory defined on a Riemann surface, obtained connecting different replicas together through a cut corresponding to the subsystem A at fixed time.The manifold thus defined has a conical singularity corresponding to the entangling surface ∂A.This method, originally designed for analytical calculations, has also been implemented in lattice simulations .On a Euclidean lattice, the replica geometry is directly implemented by means of different boundary conditions along the Euclidean-time direction for the subsystem A and for its complement B, namely the period of A is n times the period of B.
Another approach that is widely used in literature to study the entanglement entropy of theories that admit a holographic dual [77][78][79] is to use the Ryu-Takayanagi formula [12][13][14][15]46], which relates the entanglement entropy of a theory defined on the boundary to minimal-area surfaces extending in the bulk.This approach led to important predictions for the behavior of entanglement in strongly coupled gauge theories.In particular, in ref. [80] it was conjectured that entanglement can be a probe to study the phenomenon of confinement (see also refs.[81][82][83][84]); this possibility has been investigated numerically, by means of Monte Carlo lattice simulations, in SU(N ) gauge theories [53-55, 59, 60, 67, 69].
A non-trivial aspect in the study of entanglement entropy in gauge theories, however, is the ambiguity of the definition itself in the presence of a local symmetry.Gauß's law and consequently the non-local constraint defining physical states makes the physical Hilbert space of a gauge theory non-factorizable.This problem has been extensively discussed in the literature, and different definitions of entanglement entropy have been proposed [54,[85][86][87][88][89][90][91][92][93][94][95][96][97][98][99][100], among which two approaches are most commonly used: one can either extend the Hilbert space of states, including gauge non-invariant ones to ensure factorizability; or work at the level of the algebra of observables defined over a specific region, specifying which generators have to be included in region A and which ones in its complement.It turns out that for some choices of the algebras the two approaches are equivalent [95], and lead to the following general result: the trace which defines the entanglement entropy splits into the sum of different contributions, labeled by the eigenvalues of the operators at the boundary between the system A and its complement, which in general belong to the center of the algebra of observables of A. As a result, in Abelian gauge theories defined on the lattice the entanglement entropy admits the following decomposition where k labels different sectors, determined by the eigenvalues of the boundary operators, p (k) denotes the probability for the sector k, while S(ρ A ) is the Shannon entropy of the reduced density matrix restricted to the sector k.For theories based on a non-Abelian gauge group, eq.(1.4) includes a further term, which is related to the dimension of the representation of the electric flux at the boundary.
From eq. (1.4) one might worry that the entropy thus defined is not really quantifying the entanglement between two subsystems, since only the second term is distillable, while the first one is related to the classical probability distribution of the eigenvalues of the boundary operators, which depends on the specific (and somehow arbitrary) choice of operators that are included in A. While this ambiguity is indeed present on the lattice, in refs.[87,101,102] it was shown that, in the continuum limit, suitably regularized quantities, such as the mutual information or relative entropies, as well as the universal information encoded in the entanglement entropy, do not get contributions from the classical part of eq.(1.4) (which depends on high-momentum modes localized close to the boundary), and, therefore, are independent from the choice of boundary operators.
These subtleties in the definition of the entanglement entropy imply that the replica trick itself presents some ambiguities for gauge theories; in particular, working with the lattice regularization in the Feynman path-integral approach, the problem translates in how to treat the entangling surface where different replicas are joined together.
In this work we address one facet of this problem, focusing on the construction of a replica space for gauge theories through a duality transformation.For gauge theories defined in three spacetime dimensions, the latter can be carried out through a Kramers-Wannier transform [103,104], which maps the gauge theory to a spin system.It is well known that such dualization is not straightforward when the gauge theory is defined on a space with a boundary and there exist operators (including, in particular, those defining Gauß's law) defined on the boundary.However, it was shown in different works [97,[100][101][102] that, even though the entanglement entropy is not preserved by the duality transformation, its universal part, encoded in the entropic c-function, is exactly mapped in the continuum limit.Note that the investigation of quantum lattice models by means of duality was recently addressed (in the Hamiltonian formalism) in ref. [105], albeit with a focus on aspects different from the ones of the present work.
In the first part of this work, we explicitly perform the duality transformation starting from a spin model defined on a lattice in three spacetime dimensions with an arbitrary positive integer number of replicas, to derive the corresponding geometry for the dual gauge theory.As we discuss in what follows, the physical results we obtain from our construction are purely geometric, namely they do not depend on the degrees of freedom of the theory (nor on the details of how they are assigned to A or B); rather, the duality transformation gives us a well-defined prescription on how to implement a conical singularity on the lattice.The logical scheme of how the direct and the dual formulation of a theory are connected with each other, in the Hamiltonian and in the Lagrangian formalisms, is summarized in figure 1 (which also includes the references to previous literature, in which the various aspects discussed herein were studied).
Then, in the second part of this work, we make use of the Kramers-Wannier duality in a complementary way: by means of Monte Carlo simulations, we calculate the entropic c-function of the second Rényi entropy of the three-dimensional Ising model, which is dual to the Z 2 gauge theory.This allows one to study the entanglement content of the ground state of a strongly coupled gauge theory and to test the predictions that were derived from the gauge/gravity correspondence in ref. [80] in a particularly simple (albeit non-trivial) system.The results we obtain are consistent with such predictions and with similar recent studies in low-dimensional gauge theories [106].
The structure of the article is the following.In section 2 we review the generalities about the duality transformation, and discuss its formulation in the presence of replicas; we mostly focus on the case of the Ising model, considering both the case of two and three spacetime dimensions, but we also discuss the extension of our results to other gauge theories, including U(1) and Z N gauge theories (and briefly comment on the generalization to higher dimensions and to non-Abelian gauge theories).Section 3 illustrates the connection between the Hamiltonian and the Lagrangian formulations in the presence of replicas, with a detailed discussion of the (1 + 1)-dimensional Ising model case study, while section 4 is devoted to a Monte Carlo study of the entropic c-function of the three-dimensional Z 2 gauge theory.Finally, in section 5 we discuss our findings and comment on possible Casini and Huerta [101]; Radičević [97]; Moitra, Soni, and Trivedi [102] Lin and Radičević [100];

Figure 1:
Starting from a spin model in the Hamiltonian formalism, one can derive the replica geometry of the dual theory following two different paths, which yield equivalent physical results.One consists in discretizing the replica space first, and then performing the duality transformation in the Lagrangian formulation: this is worked out in this article in section 2. Alternatively, one can transform the spin system at the Hamiltonian level (effects of the duality transformation on universal information encoded in the entanglement entropy were discussed in refs.[97,101,102]), introduce the replicas as described in ref. [100] and discretize the resulting system.This approach is discussed in this work in section 3. generalizations thereof.The appendix contains a discussion of the rôle of boundary conditions in the duality transformation (section A), a study of the dual replica space for other two-dimensional spin models (section B), a derivation of the spin theory dual to the threedimensional U(1) gauge theory (section C), and some technical details about our numerical simulations (section D).

Duality transformation
In this section we discuss the duality transformation for a Euclidean lattice including replicas, with a specific focus on the Ising model in two and in three dimensions.We first review the generalities about the Kramers-Wannier duality in subsection 2.1, then in subsection 2.2 we present the detailed derivation of the duality transformation for a replica geometry, also commenting on the physical and unphysical contributions to the Rényi entropies defined in the direct and in the dual formulation.Finally, we briefly discuss the generalization to other Abelian models in subsection 2.3 and the case of higher dimensions and non-Abelian theories in subsection 2.4.

Review of the duality transformation
Let us first review the Kramers-Wannier duality transformation, following ref.[107]; for simplicity, we discuss the case of the Ising model defined on an isotropic square (in D = 2 dimensions) or cubic (for D = 3) lattice with |Λ| sites.In what follows, we consider a lattice of arbitrarily large extent in all directions, so that the choice of boundary conditions is immaterial; the case of a finite-extent lattice, with either periodic or antiperiodic boundary conditions in each of the main directions, will be discussed in the appendix A.
Starting from the partition function of the Ising model where the variables σ i ∈ {1, −1} are defined on the sites of the lattice, and µ ranges over the main directions of the lattice, by defining the function (which equals cosh β for k = 0, and sinh β for k = 1), and introducing integer-valued variables k i,µ , associated with the oriented bonds between nearest-neighbor sites on the lattice, eq.(2.1) can be recast in the form By rearranging the product of σ i terms to isolate each of them, one obtains where i k i,µ denotes the sum over the k i,µ on all links touching site i, and δ 2 is the Kronecker delta modulo 2. The constraint expressed by the latter can be enforced by means of a (discretized) Bianchi identity, by introducing a dual lattice, whose sites are displaced by half a lattice spacing in each direction with respect to the sites of the original lattice.To this purpose, let us first define In D = 2, assigning s i ∈ {1, −1} variables to the sites of the dual lattice, the constraint in eq.(2.4) is then satisfied by setting Note, however, that the mapping between configurations of the s variables and configurations of the k variables is surjective, but not injective: configurations of the s variables that differ by an overall sign flip correspond to the same configuration of the k variables.
Taking into account this overcounting by a factor 2, and using eq.(2.6) in eq.(2.2), the partition function takes the form with µ, ν and λ all different. (2.8) With manipulations analogous to those for the D = 2 case, and denoting the product of four U i,µ variables around a plaquette as U □ , one obtains . (2.9) As in D = 2, the mapping between configurations of the U variables and configurations of the k variables is not a bijective one.In this case, each configuration of the k's corresponds to 2 Ng configurations of the U 's, that differ from each other only by gauge transformations (i.e., by flips of the signs of the U 's on all of the links touching a site): here N g denotes the number of different U variables that can be set to a chosen value (say, 1) purely by gauge transformations, i.e., the number of links in a maximal tree on the lattice; note that N g is strictly smaller than |Λ|.Equation (2.9) shows that, up to an overall normalization factor, the partition function of the D = 3 Ising model evaluated at coupling β equals the partition function of the three-dimensional Z 2 gauge theory, evaluated at coupling β ⋆ .

Duality in the presence of replicas
Let us now discuss the generalization of the duality in the presence of a set of replicas.Note that our presentation in this subsection is strictly a geometric one, while section 3 contains a more detailed discussion about how a more algebraic approach can be mapped to the construction that we present here, and how the various degrees of freedom of the theory are encoded in this setup.The geometry we are interested in is the one of n replicas connected through a cut, and let us first consider the two-dimensional case, as sketched in figure 2. In a replica geometry, the expression (2.4) still holds, since its derivation is only based on the nearest-neighbor nature of the interaction and is independent from the geometry of the lattice.The difference in the duality transformation is in the representation of the k i,µ variables that satisfy the constraint; clearly, the introduction of a cut does not affect the duality transformation for links not crossed by the cut, since the transformation only involves nearest-neighbor spins in the dual lattice and the links in the direct lattice separating them.Moreover, the duality transformation remains unchanged also for links crossed by the cut that are far from the cut boundary. 1The only difference arises at the edges of the cut; in particular, in two dimensions the dual lattice has sites exactly at (the lattice analogue of) the branch points of the replica space.Consider, for example, the cartoon in figure 3, where two of the n replicas are depicted: the equation expressing the constraint associated with the site i of the direct lattice, in the replica labeled 1, is Proceeding as in subsection 2.1, and denoting the site of the dual lattice in the lower left corner of the plaquette surrounding i as j, one gets and similar relations on the other replicas.Clearly, the constraint is satisfied if and only if the spins defined on the branch-point sites of all replicas coincide: s meaning that the sites lying on (the lattice counterpart of) the conical singularity of the replica manifold belong to all replicas, hence on a square lattice they have 4n nearest neighbors.Thus, the partition function of the system with n replicas can be written as where Z ⋆ n (β ⋆ ) denotes the partition function of a system in which the spins defined on the branch points have 4n nearest neighbors.In fact, the enlarged number of nearest neighbors close to a branch point has been pointed out before [108].Note that for a finite lattice with specific boundary conditions eq. (2.11) is expected to get contributions from different topological sectors, as discussed in detail in section A of the appendix.
Using eq.(2.7) and eq.( 2.11), one can derive an explicit relation between the Rényi entropies in the direct and in the dual formulation of the model: which leads to therefore the Rényi entropies evaluated in the direct and in the dual formulation are not equal: note, however, that their difference can be traced back to the factor 1 2 arising from the change of variables from k to s and therefore it is related to the lattice degrees of freedom we are using to describe the model, which is clearly an unphysical piece of information.By contrast, the entropic c-function, which is a quantity encoding physical information, is unaffected by the additive constant ln 2 and we get where l denotes the length of the cut.
It is straightforward to carry out the dualization in the presence of replicas for the three-dimensional Ising model, too.In this case (see the sketch in figure 4) the "cut" is actually a surface, and the conical singularity is a line where the cut ends.As for the two-dimensional model, the way the duality maps degrees of freedom in the direct and in the dual formulation is not affected by the presence of conical singularities, for spins far Duality transformation in a three-dimensional replica space.The solid line is the entangling surface along the links of the dual lattice, the red region is the cut, the red dot the intersection between the link between the sites i and i + 0 and the cut.
enough from them.Consider, for example, the site i in the direct lattice shown in figure 4; in the first replica, the Kronecker-delta constraint is of the form Repeating the construction carried out in subsection 2.1, one obtains j+ 0,2 U (1) Inspecting the previous equations we can notice that gauge fields U (1) j+ 0,2 and U (2) j+ 0,2 , the ones lying on the conical singularity, appear only once in eq.(2.15).Therefore the constraint is satisfied if j+ 0,2 .This is true for all gauge fields on the entangling surface.Therefore in the dual replica space the gauge fields lying on the conical singularity are common to all replicas and they appear in the action in 4n different plaquette terms, 4 for each replica.Note that this construction can be performed independently for each link of the entangling surface, and is independent from the shape of the latter.
The presence of gauge fields shared among different replicas leads to some constraints on the gauge transformations that leave the physical properties of the system invariant.In particular, there exist plaquette terms that can be schematically written as U edge U (r) staple , where U edge is the gauge field on the conical singularity and U (r) staple is a staple in replica r forming a plaquette together with U edge .The presence of such terms require that gauge transformations performed on the sites along the entangling surface to act on all replicas.
More precisely, a gauge transform on a site on the entangling surface acts upon (4n + 2) gauge fields: the 4n links orthogonal to the conical singularity attached to that point, and the two links along the entangling surface shared by all replicas.This restriction on the possible gauge transformations in a replica space was also noted in the context of SU(N ) gauge theories [98,109].
An interesting implication of the geometry of the entangling surface is that it is possible to define Wilson loops belonging to different replicas, passing through sites on the conical singularity.These loops are closed in the full replica space, but looking at the single replica they appear to be open, with both endpoints on the entangling surface.
We can now write down the expression of the duality transformation in replica space: for a connected entangling surface made of |∂A| lattice sites, we have where Z Z 2 n is the partition function of the gauge model in the geometry that we discussed.Note that the factor 2 (n−1)(|∂A|−1) is a correction related to the number of gauge transforms which can be performed in this geometry.As for the two-dimensional case, also here for a sufficiently large system one can neglect the rôle of boundary conditions and the non-trivial topological implications they entail.
Combining eq.(2.9) with eq. ( 2.16), one can derive the expression for the Rényi entropy of the dual gauge theory (2.17) The difference between the Rényi entropies in the gauge and in the spin theories is the presence of an ultraviolet-divergent term proportional to the length of the entangling surface, therefore contributing to the area law.This term makes the Rényi entropy of the gauge theory larger than the one of the Ising model, as was already pointed out in the literature [110].The origin of this difference stems from the constraint on the number of possible gauge transforms that can be defined in a system with replicas, and results into a multiplicative constant affecting the partition function for the dual gauge theory, which obviously does not carry any physical information.Indeed, one can also note that the term |∂A| − 1 depends on the gauge choice at the entangling surface, as was also discussed in refs.[87,92,99].Like in the D = 2 case, physical information can be obtained from a suitably defined entropic c-function, in which the area-law term cancels out; in D = 3 if the length of the entangling surface does not depend on the size of the region A, as in the slab geometry, the entropic c-function can be simply defined according to eq. (1.2), and one gets Even though the first term appearing in eq.(2.17) is unphysical, it is worth inspecting it further.In particular, in the high-temperature limit of the spin model, β → 0 (which corresponds to the weak-coupling limit β ⋆ → ∞ of the gauge theory), the Rényi entropy reads which is the result obtained with different definitions of entanglement entropy in Hamiltonian lattice gauge theories [54,85,87,88,92,94,95,97,100]; in particular the − ln 2 term appearing on the right-hand side of eq.(2.19) corresponds to the topological entanglement entropy [3,111].
It is also worth noting that, in the opposite limit (β → ∞, or β ⋆ → 0) the result one would naïvely obtain is not the same as in the Hamiltonian formalism: (which is different from the expected value S Z 2 n (β ⋆ = 0) = 0): this mismatch is simply due to the fact that in the ordered, low-temperature phase of the Ising model it is no longer justified to assume that the topologically inequivalent sectors induced by the different possible boundary conditions can be neglected in the duality transformation.We present a more detailed discussion about this point in section A of the appendix.
Before concluding this subsection, we would like to point out a possible generalization of the duality transformation for the three-dimensional Ising model in the presence of replicas.
So far, we assumed the entangling surface to lie on the lattice dual to the Ising-model lattice, that is, to be placed along the links of the dual gauge theory.However, it is also interesting to analyze a geometry in which the entangling surface lies on the links of the spin model; in this case there are spin degrees of freedom on the conical singularity.From the discussion of the D = 2 case, we know that such degrees of freedom are shared among all replicas, and have a larger number of nearest neighbors, which here is (4n + 2), of which 4n in the directions orthogonal to the entangling surface and 2 along the conical singularity.
As before, let us focus on the duality transformation close to the entangling surface.Consider a link in the direct lattice connecting two spins on the conical singularity, which is shared by all replicas: the dual plaquette orthogonal to that link has an edge, oriented in the Euclidean-time direction, that crosses the cut, therefore the plaquette, before closing on itself, has to loop around the conical singularity n times.The link variable in the direct lattice is therefore represented as where U □×n denotes the closed loop of length 4n that winds around the conical singularity.Note that the presence of two such loops of length 4n is necessary to satisfy the Kronecker-delta constraint around a point of the direct lattice with an enhanced number of nearest neighbors: for the 4n links touching a site on the conical singularity one can make the usual duality transformation according to eq. (2.8), but then there are only two link variables to be adjusted, such that in the sum i k i,µ all gauge variables appear twice.As a consequence, the two k variables along the conical singularity must be built from a larger number of gauge variables.Note that this geometry was already discussed in ref. [93], whose authors carried out a strong-coupling calculation of the Rényi entropy in this "central-plaquette" geometry, where the central plaquettes are the ones encircling the conical singularity, as sketched in figure 5.
In this case, the relation between the spin and the gauge partition functions reads where Z n and Z Z 2 n are the partition functions for the spin system and for the gauge model with the geometries discussed in the present section, that are different from the ones appearing in eq.(2.16).In this central-plaquette geometry, the relation between the Rényi entropies associated with the spin model and with the gauge theory reads which is different with respect to eq. (2.17), but, again, the difference consists of unphysical terms and eq.(2.18) still holds.

U(1) and Z N models
The construction described above can be readily generalized to other Abelian gauge theories defined in three spacetime dimensions.In particular, here we will discuss the case of the U(1) gauge theory (and briefly comment on the generalization to Z N gauge theories, too).The U(1) gauge theory is defined by the action where , and each link variable U i,µ is an element of the U(1) group.We can equivalently write U i,µ = exp(iθ i,µ ) and directly work with the algebravalued variables θ i,µ , having defined the discretized exterior derivative of the θ variables as θ □ ≡ θ i,µ + θ i+μ,ν − θ i+ν,µ − θ i,ν .The theory admits a dual description in terms of a non-compact, integervalued spin model, with degrees of freedom {ϕ}; an explicit derivation of the duality is reported in section C of the appendix.The partition function of the spin model reads where I ν (β) is a modified Bessel function of the first kind.Our purpose is to use the duality transformation to derive the (geometric) structure of the replica space of a gauge theory, starting from a spin system.It is therefore useful to reverse the derivation presented in section C of the appendix to obtain the U(1) model as the theory dual to the spin system with the partition function in eq.(2.26), before extending it to a replica geometry.Using a well-known property of modified Bessel functions of the first kind the partition function in eq. ( 2.26) can be rewritten as Now k i,µ is a continuous and periodic link variable and the argument of the Kronecker delta involves the sum over the six links touching site i.The constraint is satisfied if one takes where □ is the orthogonal plaquette in the dual lattice.
In the presence of replicas, the duality transformation can be generalized along the lines that we discussed for the three-dimensional Ising model; the geometry we consider is the same as in figure 4, but with ϕ variables instead of σ ones.The constraint associated to the site i of the left replica is again expressed by eq.(2.15).Note that, with the cut placed as in figure 4, all links along the 0 direction (both in the direct and in the dual lattice) intersect the cut, and therefore are shared between different replicas.As a consequence, a plaquette in the dual lattice that lies along the 0 direction has some links in replica 1 and some in replica 2. Imposing eq.(2.31) as a solution of the constraint in this geometry, we obtain the following equations Using these equations, the constraint is identically satisfied if and only if (2.32) (and similar relations for other replicas), i.e., if the gauge variables lying on the entangling surface are the same in the different replicas.Also in this case we obtain the same picture as in the Z 2 gauge theory: gauge fields on the conical singularity belongs to all replicas, and they belong to 4n different plaquettes.
The construction is the same for the case of Z N gauge theories, too.Also, it is easy to see that, one can also recover the central-plaquette geometry for Z N and U(1) gauge theories starting from a spin model in which some spins lie on the conical singularity.
We conclude that, in all of the Abelian gauge theories we analyzed, the geometry of the replica lattice is the same, regardless of the details about the degrees of freedom.In particular, it is found that gauge fields on the entangling surface belong to a larger number of plaquettes.

Higher dimensions and non-Abelian theories
We conclude this section with some comments on generalizations of the previous analysis to other theories.
The generalization to higher dimensions is straightforward.In spacetime dimension D the entangling surface is a codimension-2 manifold.Any s-simplex lying on the entangling surface has a larger coordination number, i.e., a larger number of (s + 1)-simplices touching it.For example, a site on the conical singularity is connected to (4n+2D−4) sites: 2(D−2) on the entangling surface and 4n in the two transverse directions.A link on the entangling surface belongs to (4n + 2D − 6) different plaquettes, and so on.
The central-plaquette geometry can be generalized to higher dimensions, too.Consider an s-simplex with at its boundaries a set of (s − 1)-simplices, with s > 1.If the simplex encloses the conical singularity, the size of its boundary is enlarged, since it has to wind around the singularity n times before closing on itself.
An important and less simple question is, whether it is possible to extend this construction to non-Abelian gauge theories.Duality transformations for non-Abelian gauge theories are known, but they are considerably more involved [112][113][114][115][116][117][118][119], and it is not straightforward how to apply them to derive a geometric construction for the dual replica space.
Note however that, from a purely geometric point of view, the key features that we discussed in this section for the construction of the dual replica space (namely, the enlarged number of nearest neighbors and the central-plaquette setup) can be directly implemented-and, in fact, have been implemented [93]-to study Rényi entropies in arbitrary gauge theories.Our analysis shows that, at least for Abelian gauge theories, such geometries are justified by the duality construction.While we have not proven that these geometries are the correct ones also for non-Abelian theories, the observation that the geometry of the replica space appears not to be dependent on the gauge group suggests that this may be the case.

Connection with the Hamiltonian formalism
In section 2 we applied the duality transformation to derive the geometry of the replica space for Abelian gauge theories, working in the Lagrangian formalism.Our derivation followed the path depicted in figure 1, from the lower left corner to the upper right corner, passing through the upper left corner: one starts from a quantum lattice model in the Hamiltonian formalism, where there exists a well-defined notion of Hilbert space, and chooses a factorization H = H A ⊗ H B .For the spin models considered here, such as the Ising model, the factorization can be performed without ambiguities.At this stage one can discretize the path integral that appears in the replica trick, namely where r labels the different replicas and σ A r σ B r denotes a complete set of states (where we made the bipartition of the system explicit).In the simple case of the (1+1)-dimensional Ising model, the Hamiltonian in each replica can be written as the sum of three terms, H A , H B , and h: the first two respectively describe the interactions between degrees of freedom inside A only and in B only, while h represents the interactions between the degrees of freedom in A and those in B: (3.2) σ z and σ x are Pauli matrices and, in the boundary term σ z a σ z b , a and b are nearest neighbors.From eq. (3.1) one can obtain the partition function of a statistical-mechanical system by dividing the (compact) time direction in N τ intervals of length a τ = 1/(T N τ ), Trotterizing the exponential and inserting a resolution of the identity at each discrete time step.If we take, for definiteness, the σ A r σ B r to be eigenstates of σ z , the hopping term σ z i σ z j leads to a nearest-neighbor interaction in the spatial direction, while the σ x operators generate interactions in the temporal direction.The replica structure shown in figure 2 emerges from terms such as where σ A r σ B r ; τ is a state inserted at Euclidean time τ .The previous term leads to an interaction in the time direction between two spins in the same replica, if they belong to B, and between two different replicas, if they belong to A. From the discrete replica space one can then obtain the dual geometry as described in section 2.
In this section we discuss how one can derive the same geometric picture by reversing the order of the steps, i.e., by first performing the duality transformation in the Hamiltonian description, and then discretizing the resulting path integral.For simplicity from now on we will discuss the case of the (1 + 1)-dimensional Ising model, but the discussion is easily generalizable to other spin models and higher dimensions.
For this purpose, instead of working with states in a Hilbert space, it is more convenient to work with algebras acting on that Hilbert space [87].Given a Hilbert space H and an algebra of observables A, a bipartition of the system is now defined by identifying a subalgebra A A of operators acting on the subsystem A. The reduced density matrix is then defined as the unique operator ρ A ∈ A A satisfying for all O ∈ A A .In general ρ A can be written as a linear combination of the generators of the algebra.For the (1 + 1)-dimensional Ising model on a chain of N sites, the maximal algebra one can define is generated by the two Pauli operators σ x and σ z on each site The choice of a subalgebra which leads to eq. (3.1) is the following: let A be the set of sites with 1 ≤ i ≤ M < N , then which is the maximal subalgebra we can define on the subchain from site 1 to site M .One can then act on the subalgebra defined in eq.(3.6) with the Kramers-Wannier duality [97]: the subchain of M sites maps to a subchain of (M +1) sites whose edges lie on the entangling surface, and the algebra A A maps to the non-maximal algebra (in this section we label dual sites with half-integers) where the σ x generators at the edges are removed.As a consequence, this algebra has now a center, i.e., a set of generators commuting with all the other generators of the algebra.Any operator belonging to an algebra with a center must be block-diagonal in a basis that diagonalizes the center, therefore the Hilbert space splits into superselection sectors, labeled by the eigenvalues of the generators belonging to the center.Therefore, states with reduced density matrices in this algebra must be eigenstates of the generators of the center, in this case σ z 1/2 and σ z M +1/2 .The path-integral representation of Tr ρ n A in the presence of superselection sectors was analyzed in ref. [100].The key difference with respect to eq. (3.1) is that one has to keep different superselection sectors separated in Tr ρ n A .One can first compute the path integral in a fixed sector, then sum the result over different sectors.Importantly, the sector is the same for all replicas.For the algebra defined in eq.(3.7), this leads to In this expression the eigenvalues of σ z 1/2 and σ z M +1/2 in the summand are fixed to be s and s ′ respectively, and they are the same for all replicas.The h term in the Hamiltonian, describing interactions between pairs of degrees of freedom of which one is in A, the other in B, can be written as The difference with respect to the previous case is in the presence of terms such as that is, for all replicas spins at sites i = 1 2 and i = M + 1 2 for τ = 0 have the same value in all configurations that contribute to the partition function.Note that this happens only for τ = 0, since at all other times resolutions of the identity have been inserted, leading to independent sums over all eigenvalues of σ z on each site.
Given that the spins on the entangling surface have the same value in all configurations, they are effectively the same spin: this is how a geometry including two spins with an enlarged number of nearest neighbors emerges directly in the Hamiltonian formalism.
The derivation of the replica space in the presence of superselection sectors discussed in this section is completely general, and can be applied to other spin systems, as well as to higher-dimensional theories.In higher dimensions the explicit derivation of the discretized path integral starting from the Hamiltonian of the theory is less straightforward, and we do not address the calculation here.However, it would be interesting to study how different choices of algebras for gauge theories map to different replica spaces.In particular, the main choices one can make are two: the electric-center algebra and the magnetic-center algebra [87]; it seems natural to expect that these two choices correspond to the two different geometries we found in section 2, namely the central-plaquette geometry and the one with the enlarged number of nearest neighbors, while combinations of the two choices of the center might map into combinations of the previous two geometries.
As a final comment, we stress that an algebraic choice in the Hamiltonian formalism is mapped to a purely geometric one in the Lagrangian formalism.In particular, different choices of the algebra A A correspond to different geometries of the entangling surface.In both pictures, the differences between the various choices seem to be related to the ultraviolet physics of the model, and this is consistent with the expectation that in the continuum limit the universal part of the entanglement entropy, as well as properly regularized physical quantities such as the mutual information, are independent from the choice of the algebra.

Monte Carlo results for the Z 2 gauge theory in three dimensions
In this section we present the results of a Monte Carlo study of the entropic c-function of the ground state of the three-dimensional Z 2 gauge theory.An early study of the entanglement entropy in this model by means of Monte Carlo simulations was reported in ref. [54], while more recently this problem has been addressed in various studies based on tensor networks [120][121][122] and exact diagonalization [123].
As discussed in refs.[2,52], the replica trick allows one to express Rényi entropies in terms of ratios of partition functions, which are quantities that today can be computed to high precision by means of Monte Carlo simulations on the lattice.This is true also for the associated entropic c-functions, which can be written as In this work we make use of a highly efficient algorithm, based on Jarzynski's theorem [124,125], that we introduced in ref. [75], generalizing previous work done in ref. [68].Jarzynski's theorem is an exact equality between ratios of partition functions defined by equilibrium Boltzmann distributions and averages over non-equilibrium processes.Specifically, consider two probability distributions p 0 = exp(−S 0 )/Z 0 and p 1 = exp(−S 1 )/Z 1 which can be continuously connected by tuning a parameter λ, such that S(λ = 0) = S 0 and S(λ = 1) = S 1 ; Jarzynski's theorem then states that where W 0→1 is the (generalized, dimensionless) work performed to drive a system, initially at equilibrium, from λ = 0 to λ = 1, and the average is performed over an ensemble of different trajectories that the system follows in its out-of-equilibrium evolution.It is important to stress that the equality holds regardless of the specific details of how the system is driven off equilibrium: the equality holds irrespectively of the way in which λ is varied, as long as this evolution is the same for all trajectories, and that the system is genuinely out of equilibrium for every λ > 0. In recent years, this theorem has been harnessed in highprecision lattice simulations, to study a variety of different quantities [126][127][128][129][130][131], including in computations of Rényi entropies through quantum Monte Carlo simulations [70,[132][133][134][135][136][137][138].
To obtain high-precision results for the Z 2 gauge theory we exploit its duality to the three-dimensional Ising model-using, in particular, the fact that the entropic c-function computed in the two theories is the same, as shown by eq.(2.18).The geometry we simulated is a slab geometry, with the entangling surface lying along the links of the dual lattice.
A major advantage of directly simulating the Ising model is the possibility to use cluster algorithms, which essentially allow one to perform simulations arbitrarily close to critical points, without suffering from critical-slowing-down problems.In our simulations we used the Swendsen-Wang cluster algorithm [139] in a parallelized version [140][141][142] adapted for GPU-accelerated supercomputers.
The non-equilibrium protocol we used to calculate the entropic c-function is the same as the one we introduced in ref. [75].First, note that the ratio of partition functions appearing in eq.(4.1) can be approximated on the lattice in terms of Z n (l)/Z n (l + a), with a denoting the lattice spacing (i.e., the shortest length that can be resolved on the lattice).We can then design a non-equilibrium evolution connecting the two probability distributions by varying the couplings at the edge of the entangling surface: starting from a slab of length l + a, we reduce the couplings of spins in replica r near the entangling surface with spins in replica r + 1 and, at the same time, we increase the coupling between spins in the same replica, decreasing the length of the slab by one lattice spacing.
We performed simulations in a geometry of two replicas to compute the entropic cfunction associated to the second Rényi entropy.Simulating different values of β and for various values of the lattice volume allowed us to perform both thermodynamic and continuum extrapolations.To the best of our knowledge, this is the first study in which such extrapolations are performed for a (2 + 1)-dimensional gauge theory.For the scale setting we followed ref. [143], where the dependence of the critical value of the Ising coupling β c on the lattice size in the Euclidean-time direction in units of the lattice spacing, denoted as N τ , was determined with high precision for a wide range of values.The critical deconfining temperature of the dual gauge theory can then be obtained as aT c = 1/N τ,c , so that one can then express all of the other dimensionful quantities in the appropriate units of T c .To make sure that the system is not affected by undesired finite-temperature corrections, we performed our simulations on lattices whose extent in the Euclidean-time direction is at least ten times larger than the critical one.Further details on the scale setting and on the different extrapolations can be found in section D of the appendix.
Figure 6 shows a summary of our final results.The plotted data correspond to the function where C CFT   2   was estimated in ref. [144].The length of the cut is measured in units of two different mass scales, namely, the critical temperature T c and the mass m g of the lightest glueball-like state in the theory, taken from ref. [145].
The behavior of the entropic c-function of the ground state is known only for few strongly interacting theories.Some predictions for gauge theories can be derived from the holographic correspondence, using the Ryu-Takayanagi formula [12,13] to compute the entanglement entropy in terms of the area of minimal surfaces.The authors of ref. [80] calculated the entanglement entropy for confining backgrounds in non-Abelian gauge theories similar to quantum chromodynamics (QCD), and in a slab geometry at different length scales, finding the following picture: at small length scales (where, due to asymptotic freedom, the theory is invariant under conformal symmetry), the entropic c-function is constant and proportional to the number of degrees of freedom-which, for SU(N ) gauge theories in the large-N limit, scales like the square of the number of color charges N .As the length of the slab is increased, the theory eventually hadronizes, and, since the number of hadronic degrees of freedom is independent from N , the entropic c-function is then expected to exhibit a sharp transition, from a regime in which its parametric dependence on N is O N 2 to one in which it scales as O N 0 .The length scale at which this transition occurs is controlled by a mass scale associated to the boundary theory.For the QCD-like theories considered in ref. [80], all of the typical energy scales are parametrically equivalent, and so is the scale associated to the transition of the entanglement entropy.
In order to gain insight on the behavior of the entropic c-function in different regimes, we fitted it with different models.At large length scales, the theory can be modeled in terms of a gas of non-interacting glueballs [146,147].On the other hand, for (1 + 1)-dimensional free scalar and Dirac fields of mass m it was found that, for n > 1 and lm ≫ 1/(n − 1), the entropic c-function admits the expansion [148,149] In particular, for n = 2 the expansion is exact at leading order.Following ref. [80], one can think of a free scalar field in 2 + 1 dimensions as a collection of infinite two-dimensional scalar fields with transverse momentum k and mass m(k) = √ m 2 + k 2 .Summing over the entropic c-functions of each mode leads to Therefore we fitted our data with the Ansatz The function fits the data well from large values of lm g down to lm g ≃ 0.84.The results of the fit are A = 0.33(3), α = 0.360 (19), with χ 2 red = 0.82 (where χ 2 red denotes the χ 2 divided by the number of degrees of freedom of the fit).Note that the parameter α appearing in eq.(4.6) controls the exponential decay of the entropic c-function at large distances, and it reduces to 1 for a free scalar field.For the Z 2 gauge theory we obtain a value of α which is smaller than in the free-field case, albeit roughly of the same order.A discrepancy with respect to the free case is expected, since the system we are considering is interacting.
Consider now the opposite regime (small lm g ); the recent ref. [106] reported a study of the two-mode logarithmic negativity in the Schwinger model, which was studied as a measure of entanglement at different length scales.It was found that, at distances shorter than the inverse of the mass of the lightest state in the continuum, the decay of the twomode negativity is polynomial.Motivated by this observation, we fitted our data for small lm g with a power-law function The function fits well our results from small values of lm g up to lm g ≃ 1.26.The results of the fit are B = 0.360( 9), c = 0.48(2), with χ 2 red = 1.02.Note that a deviation from the power-law decay is expected at sufficiently short length scales, since C 2 (0) is finite.At very short distances, however, the simulations become more expensive in terms of computing resources, since one should get very close to the critical point of the theory and very large lattices are necessary to extrapolate to infinite volume.Despite these difficulties, looking at the data we obtained from finite lattices and for different values of the lattice spacing, as shown in figure 6, one observes clear convergence towards a finite value for lm g → 0, even though some data sets appear to tend to a value slightly different from 1.This is likely a systematic uncertainty due to the finiteness of the lattice spacing and/or of the length of the slab, which become particularly relevant effects in the lm g → 0 limit.

Conclusions
In this work we discussed the replica space and Rényi entropies of Abelian gauge theories; in particular, we exploited duality transformations relating them to spin models, and discussed the geometric features of this mapping.The advantages of this approach are twofold: on the one hand, the duality transformation can be explicitly performed in a replica space, leading to a derivation of the replica geometry which is valid for different gauge theories.On the other hand, the duality mapping allows one to compute entanglement measures for strongly interacting gauge theories by simulating spin models, which can be efficiently studied with Monte Carlo methods.
In the first part of this work we analytically performed the duality transformation for two-and three-dimensional spin systems with Abelian global symmetry group.The picture emerging from our study is a purely geometric one, in the sense that it does not depend on the specific degrees of freedom of the theory one is considering.Fields lying on the conical singularity belong to all replicas and the local action gets contributions from interactions with fields in all replicas; in a spin model, this translates into an enlarged number of nearest neighbors for the spins on the conical singularity, whereas in a gauge theory the gauge field therein contributes to a larger number of plaquettes.Note that this geometry realizes in a natural way a number of features of the replica space of a gauge theory that have already been discussed in past works, such as the fact that gauge transformations on the entangling surface must be the same on all replicas [98,109] or the presence of Wilson loops that appear to be open, if one looks at a single replica, but close in a different one [54].The existence of a larger number of neighbors in the dual formulation of the theory has also been discussed recently in refs.[150,151].
If the conical singularity is placed in the dual lattice, plaquettes encircling the entangling surface have to loop around it n times before closing onto themselves.This is the central-plaquette geometry discussed in ref. [93].In the present work we focused mostly on the first geometry, since the Ising model we simulated is dual to a gauge theory with gauge fields on the conical singularity, without central plaquettes.It would be interesting to investigate how the results change if one uses a different geometry, and to test whether the same continuum limit is recovered.
In order to make contact with the Hamiltonian formalism, in section 3 we derived the replica geometry for the (1 + 1)-dimensional Ising model in presence of superselection sectors, finding the same geometry of section 2, with spins on the entangling surface having an enlarged number of nearest neighbors.The derivation is completely general and in principle it might be extended to gauge theories with magnetic-and electric-center algebra.
In the second part of our work, we numerically exploited the duality transformation, simulating the three-dimensional Ising model to compute the entropic c-function associated with the second Rényi entropy of the Z 2 gauge theory.Despite its simplicity, this theory displays many phenomena that also characterize SU(N ) gauge theories, such as a mass gap and confinement.Combining the efficient algorithm [75] based on Jarzynski's theorem [124,125] and GPU parallelization, we were able to extrapolate our results to the infinite-volume and continuum limits, and to compare our final results with predictions based on the holographic correspondence [80].
As expected by semi-empirical arguments, the entropic c-function decays exponentially for large lengths of the slabs.At large length scales, the ground state of the gauge theory can be described as a gas of weakly interacting bosons, and indeed the exponential decay seems to be parametrically controlled by the mass of the lightest glueball, up to order-one numerical factors.
At short length scales, the entropic c-function scales as a power of the length of the slab.The exponent of the power-law decay we obtained fitting our data is compatible with 1  2 and it might be related to the scaling dimension of some operator in the conformal field theory describing the three-dimensional Ising model at criticality.Note, however that such exponent is likely to be connected in a non-trivial way to the scaling dimension of the twist fields [2,152].The transition from power-law to exponential scaling is observed to occur around a length scale controlled by the mass of the lightest glueball, that is the mass gap of the theory.This behavior resembles the transition that was found in many holographic models dual to a confining background [13,80,83] and also in Monte Carlo studies of SU(N ) gauge theories [53-55, 67, 69].Moreover, we note that similar results were recently obtained from a study of the two-mode entanglement negativity in the Schwinger model [106].
In this work we studied the simplest, non-trivial gauge theory which admits a dual spin description.A natural generalization is to study the entropic c-function of other Abelian gauge theories by exploiting the duality transformation.In this respect, a particularly interesting theory is the U(1) gauge theory in 2 + 1 dimensions.Moreover, as we previously stated, the same geometries we derived for Abelian theories can in principle be used to study Rényi entropies in non-Abelian gauge theories.We leave these directions of research for future work.Consider now a lattice with periodic boundary conditions in all directions.Using the fact that the partition function can be written in terms of the dual variables as Note that for τ ⋆ = 0 one recovers eq.(2.7).Let us now inspect the term (s i s i,µ ) τ ⋆ i,µ : if the direct link is frustrated, the dual variables orthogonal to that link appear in the dual partition function as operator insertions.Consider a model with a set of frustrations that is gauge-equivalent to a configuration with no frustrations at all; as previously noted, such sets are dual to closed paths, therefore all operator insertions appear twice in eq.(A.4) and one obtains again eq.(2.7).
It is now clear that a set of frustrations that is inequivalent to Z(β) is dual to a set of open paths in the dual lattice.This leads to the well-known result namely a generic spin correlator in the D = 2 Ising model can be expressed in terms of the ratio of the partition function with frustrations piercing paths that connect the various insertions, over the partition function without frustrations.
From the invariance of the partition function under insertion of a closed loop of frustrations which is contractible to a point, it follows that the constraint appearing in eq.(2.4) admits additional solutions, beside the one discussed in section 2, if the topology of the system is such that there exist non-contractible closed loops.In the latter case, the constraint will have different solutions, related to the inequivalent partition functions with the insertions of topological lines of frustrations along non-contractible paths.In general, if we start from the Ising model without frustrations, the general duality transformation reads For example, in the geometry of a torus where Z ij (β) is the partition function with periodic (p) or antiperiodic (a) boundary conditions in directions i and j; antiperiodic boundary conditions are equivalent to the insertion of frustrations along a non-contractible path. 2ote that topology affects the calculation of Rényi entropies, because the contributions from the partition functions with the different types of boundary conditions have to be taken into account.However, in some regions of the parameters space the contribution from antiperiodic boundary conditions is suppressed.In particular, antiperiodic boundary conditions generate phase boundaries, which are exponentially suppressed in the ordered phase.Therefore, up to numerical constants This fact can be seen also at the level of the Rényi entropies.Consider eq.(2.13): in the limit β → 0 it gives the expected result: indeed at β = 0 the ground state of the Ising model is in a product state, while for β ⋆ → ∞ the entropy takes a contribution from the two degenerate vacua.However, in the opposite limit, we do not recover the same result, rather: The ln 4 term can be related to the existence of the four different topological sectors appearing in eq.(A.7).

Ising model in three dimensions
The previous discussion can be easily generalized to three dimensions by looking for other solutions of the Kronecker-delta constraint in the case of the Z 2 gauge theory.With the same logic as before, consider the following gauge transformation of the Z 2 gauge theory: given any link of the lattice, we frustrate all four plaquettes sharing that link; it is easy to see that such transformation is a symmetry of the partition function.Moreover, the set of frustrated plaquettes is pierced by a closed path in the dual lattice.New solutions to the constraint imposed by the Kronecker delta appearing in eq.(2.4) can be obtained by choosing any set of closed paths in the lattice dual to the one of the Z 2 gauge theory and introducing antiferromagnetic couplings for all plaquettes pierced by those paths.
If the paths are all contractible, then the resulting partition function is equivalent to the partition function of the Z 2 gauge theory, therefore in a lattice with trivial topology there are no further solutions.If the topology is non-trivial, on the other hand, the general solution of the Kronecker-delta constraint corresponds to the sum of partition functions with insertions of frustrations along all of the possible types of non-contractible paths.For example, in the geometry of a three-dimensional torus, the partition function dual to the Ising model reads where Z Z 2 ijk is the partition function of the gauge theory with (−) or without (+) the insertion of a topological line of frustrations along a non-contractible path in each of the three directions.
The same derivation can be carried out if one starts from the Z 2 gauge theory to derive the Ising model as the dual theory.Like in the two-dimensional case, the three-dimensional Ising model is invariant if one frustrates all links touching one site.The difference from the two-dimensional case is that in three dimensions the set of frustrations is dual to a closed surface.Different solutions to the Kronecker-delta constraint are therefore associated with the partition functions of systems with frustrations that are dual to a set of closed surfaces.As a consequence, for a lattice with non-trivial topology one has to sum over all possible partition functions with insertions of non-contractible surfaces.In the geometry of a threedimensional torus, one then recovers the well-known result (A.12) As in the two-dimensional case, partition functions with antiperiodic boundary conditions are suppressed in the confining phase of the gauge theory (β ⋆ < β ⋆ c ), which corresponds to the ordered phase of the spin model (β > β c ). Therefore This is the reason why, in our Monte Carlo study of the entropic c-function of the Z 2 gauge theory, we neglected contributions from topologically non-trivial sectors.In addition, the previous arguments explain why in section 2 we found that the strong-coupling result expressed by eq.(2.19) is consistent with the prediction that can be derived in the Hamiltonian formalism, while the weak-coupling one is not.In the weak-coupling regime, the result is spoiled by the effects of non-trivial topology, which are no longer suppressed.

B Duality in two-dimensional Z N and U(1) models
In this section we discuss the duality transformation in systems with Z N and U(1) symmetry in a replica space, as a natural generalization of our analysis for the Ising model.A relevant question we want to address here is whether the duality on a replica manifold is affected by the nature of the degrees of freedom of the theory, or just depends on geometric aspects, namely the presence of branch-point singularities.Therefore, for these models we will not keep track of the numerical prefactors in the partition function and we will just focus on the geometry of the dual lattice.We focus on a specific class of models, the so-called clock models, defined by the partition function where q ∈ {0, 1, . . ., N − 1} and ∆ µ q i ≡ q i − q i+μ .The Boltzmann factor can be Fouriertransformed as where the C k (β) are the coefficients of the expansion.The partition function then becomes In the second line we used the discretized version of integration by parts, while in the last expression we used the Fourier decomposition of the Kronecker delta enforcing the constraint ∆ µ k i,µ = 0 mod N .Let us start by considering the simplest case of an infinite lattice without any replica structure.
It is easier to solve the constraint in terms of another set of variables k, which are equal to k mod N and satisfy ∆ µ ki,µ = 0 directly.Since there is a unique change of variable from k to k (up to a Z N transformation), the representation is the same in terms of k.Using the k variables, the constraint can be solved by introducing a set of integer-valued spins ϕ on the sites of the dual lattice and imposing ki,µ = ϵ µν ∆ ν ϕ.This expression can be manipulated further; for N = 2, 3, or 4, the model is self-dual.Consider now a system of n replicas.It is clear that also in this case the duality is only affected by the presence of branch points, therefore we focus our attention on the spins close to one edge of the cut.The geometry is the same as in figure 3, but with the variables s replaced by ϕ.We can solve the constraint at site i imposing eq.(B.6), k (1) i,0 = ϕ which satisfy ∆ µ ki (1) ,µ = 0 only if ϕ (1) j+ 0 = ϕ (2) j+ 0. We see again that the dual spin defined on the branch point is shared by all replicas and has a larger number of nearest neighbors.Note that this geometric picture does not depend on the specific choice of the symmetry group.
The case of a model with U(1) symmetry is straightforward, since the partition function is equivalent to eq. (B.1) with continuous variables, and the duality transformation is performed following the same steps, leading again to a geometry where the spin on the branch point singularity has a larger number of nearest neighbors.

C Derivation of the model dual to U(1) gauge theory
In this section we derive the expression of the partition function of the spin theory dual to the U(1) gauge theory in three dimensions, eq.

D Details of the simulations
In order to convert our results from lattice units to physical units, we set the scale using the study of the Ising model in three dimensions reported in ref. [143], where the values of β at which the deconfinement transition occurs were computed numerically, for a range of temporal extents up to N τ = 16.To perform simulations even closer to the critical point and access finer lattices, we extrapolated the results of ref. [143], using eq.( 22) in the same paper, to larger N τ values, as shown in table 1.Note that here and in the following we are implicitly trading the a → 0 limit for the equivalent N τ → ∞ limit.For additional information on our simulation, see table 2. At fixed values of the lattice spacing and for fixed length of the slab we performed an infinite-spatial-volume extrapolation, keeping N τ fixed to at least 10N τ,c and N τ larger than the spatial extent (in units of the lattice spacing) N s .For values of β close to the critical point, finite-volume effects are non-negligible, and an extrapolation is necessary, whereas for larger values of (β − β c ) we were able to simulate physical volumes that are compatible with the infinite-volume limit within their uncertainties.This allowed us to bypass the need of extrapolation to the thermodynamic limit for simulations of large slabs, where the signal is poor and large statistics is required.
The fitting function we used for the infinite-volume extrapolation, shown in figure 8, is f (L; c, A, M ) = c + A exp(−M L), (D.1) where L = aN s .The fit parameter c is the infinite-volume result, while M is a mass parameter which we imposed to be the same for all the extrapolations, by performing a global fit.The value we obtained is M/T c = 1.31(2) with a global χ 2 red = 1.87.The fit function in eq.(D.1) is empirical, since no study on the scaling of the entropic c-function at finite volumes exists in literature.A better understanding of the scaling properties of such non-local quantities would be important, to have better control on the extrapolations from results at smaller volumes.
The continuum extrapolation was performed in two different ways, as shown in figure 9.The first one consists in approximating the derivative with the backward difference, which, in principle, has O a 2 discretization errors.However, with this approach the different points we use for the continuum extrapolations are evaluated at (l + a/2)m g , which tends to lm g as a goes to zero but introducing further discretization errors which are linear in a.We took the absolute difference between the continuum extrapolations carried out with these two methods as an estimate of (a source of) the systematic error on our results.The final result reported is the value of the mid-point approximation with an uncertainty obtained by summing in quadrature the statistical and the systematic errors.

Figure 2 :
Figure2: Two-dimensional replica space.The vertical direction is the Euclidean-time direction, while the horizontal one is the spatial one.The red wavy line represents the cut; green links connect spins in different replicas, therefore the periodicity of the lattice in the subsystem A is n times the periodicity in the subsystem B.

Figure 3 :
Figure 3: Structure of the dual lattice close to a branch point (red dot).

2 Figure 5 :
Figure 5: Central-plaquette geometry.The blue loop crosses the cut (pink surface) and therefore has length 4n, where n is the number of replicas.

Figure 6 :
Figure 6: Monte Carlo results for the entropic c-function of the second Rényi entropy.Inset: entropic c-function at small values of lm g .Blue symbols denote the results extrapolated to the thermodynamic and continuum limits, while the other symbols correspond to results obtained from finite lattices.

Figure 7 :
Figure 7: Frustrations in the direct lattice (wavy links) and the corresponding closed paths in the dual lattice (dashed lines).

Figure 8 :
Figure 8: Thermodynamic limit for different values of the lattice spacing and length of the slab.

Table 1 :
Values of β associated to each N τ,c .This approximation leads to discretization errors that are linear in the lattice spacing.The second method is to use the mid-point approximation,

Table 2 :
Details of our simulations.First column: length of the slab in units of the inverse of the critical temperature.Second column: different critical temporal extents (i.e., lattice spacings) used for the continuum extrapolation; for the values marked with a star, we calculated the entropic c-function for two different lengths of the cut and interpolated linearly to the corresponding lT c .Third column: ranges of spatial extents used for the thermodynamic extrapolation.