Disentangling the gravity dual of Yang-Mills theory

A construction of a gravity dual to a physical gauge theory requires confronting data. We establish a proof-of-concept for precision holography, i.e., the explicit reconstruction of the dual background metric functions directly from the entanglement entropy (EE) of strip subregions that we extract from pure glue Yang-Mills theory discretized on a lattice. Our main focus is on a three-dimensional Euclidean SU(2) theory in the deconfining phase. Holographic EE suggests, and we find evidence for, that the scaling of the thermal entropy with temperature is to power 7/3 and that it approaches smoothly the critical point, consistent with black hole thermodynamics. In addition, we provide frugal results on the potential between quenched quarks by the computation of the Polyakov loop correlators on the lattice. Holographic arguments pique curiosity in the substratum of Debye screening at strong coupling.


Introduction and summary
Our current understanding of fundamental particle physics rests on non-Abelian gauge field theories. Simplest examples among these are pure Yang-Mills (YM) theories, the studies of which constantly contribute to our ever-improving understanding of nonperturbative aspects of Standard Model of particle physics. The lucrative simplicity of YM theories void of fermionic matter is deceptive, however. The prime example is the inability to currently extract a satisfying narrative of the deconfinement phase transition. Attempts to extrapolate obtained lessons to situations in the presence of fundamental quarks described within the theory of strong interactions, quantum chromodynamics (QCD), falter. A description for the deconfinement phase transition is inherently nonperturbative, which means that standard diagrammatic methods based on weak coupling expansions may not be optimal starting points. On the other hand, lattice approach hits the challenge of simulating fermions in the fundamental representation of the gauge group. New approaches to shed light on this immensely important and highly convoluted issue are certainly welcome.
Here we do not directly address the phenomenon of deconfinement phase transition. Rather, we wish to point out another outstanding issue in pure YM theories, the lack of understanding of the behavior of entanglement entropy (EE), S EE . Entanglement entropy may though play a role in our quest to understand the deconfinement phase transition [1,2] as it is a probe of the underlying energy scale or the finite correlation length [3], an intrinsic dynamical one in the cases of interest to us in this work. In somewhat similar fashion to the case of deconfinement phase transition, computations of EE in non-Abelian gauge field theories are highly complicated. Only existing results directly from the field theory for slab subregions have been obtained at zero temperature from lattice simulations of (3+1)dimensional pure YM theories with SU(N c ), N c = 2, 3, 4 [4][5][6][7][8][9] as well as at around critical temperature in SU(3) [7]. 1 In addition, note that the (1+1)-dimensional case can be solved exactly since the Migdal-Kadanoff decimation approach [17,18] is exact. This leads to trivial entanglement entropy [5,19], i.e., in the sense of not depending on the spatial extent of the subregion; in Sec. 2 we provide a complementary exposition for this fact. In the context of gauge/gravity duality, the Ryu-Takayanagi proposal [20] allows extracting the entanglement entropy for spatial subregions in the large-N c limit in a variety of supersymmetric field theories in diverse dimensions at large coupling.
In this paper, we will report on extracting the lattice derivatives of the entanglement entropy in pure SU(N c ) Yang-Mills theories using a powerful new method. Our treatment is general, but we gear our attention to Euclidean three-and four-dimensional cases and consider entanglement entropies for single strip and slab geometries, respectively. Similarly to bulk thermodynamic quantities, as the entanglement entropy itself is not an observable, we will measure the derivative of the entanglement entropy with respect to the width of the strip (slab) for a given : ∂S EE /∂ . More precisely, we utilize the replica method [21][22][23], which on the lattice can be schematically summarized as follows. If we are interested in region A, we construct a reduced density matrix by tracing over the degrees of freedom in the complement, region B, ρ A = Tr B ρ, where ρ is the density matrix of the whole system, A ∪ B. The entanglement entropy,  associated with region A can then be written with the help of s replicas, Here the first equality is an identity and the second one can be understood through the connection to path integrals. Z is the partition function of the original system, where the division to A and B is only imaginary, fields are periodic in 1/T , and hence Z does not depend on (nor on s). The partition function Z( , s) is for the field configurations which have more complicated boundary conditions: for spatial coordinate values in the B region the fields are 1/T periodic and for the spatial coordinate in the region A, the fields are s/T -periodic; see Fig. 1 for a sketch. As alluded to before, we are not interested in S EE , however, but its derivative with respect to . To this end, by acting with this derivative on (1.2) we get where we replaced for the free energy F ( , s) = − log(Z( , s)). In the latter step we approximated the spatial derivative with a discrete derivative on the lattice, with lattice spacing a, and for the derivative with respect to replica number s we only consider the dominant contribution. Higher replica numbers give rise to contributions whose sizes are expected to be small [9]; in free field theories the dependence of EE and the second Rényi entropy can even be shown to be the same [22,24]. The right-hand-side of (1.3) is a free energy difference, a finite quantity that we extract from the lattice. This will be detailed in Sec. 2 and a brief introduction to lattice methods can be found in Appendix A. Notice that there are no issues of gauge invariance as long as one chooses not to pierce the cut through the sides of plaquettes [25]. Therefore, this lattice formulation does not suffer from the muddle associated with the tensor decomposition of the physical Hilbert space that one promptly runs across in gauge theories when formulating the tracing over the complement regions [26,27].
Our results for extracting precise data for entanglement entropies are based on a novel computational method. The advantage of our method is that it enables us to outperform previous simulations, based on a clever interpolation ansatz for the free energies [26] needed in (1.3) that we improve upon. This advancement mitigates the severe signal-to-noise ratio, thus allowing more accurate simulations with still moderate computational resources. We can push down statistical error margin considerably relative to previous works in the field. Our method is, in principle, applicable for any number of colors (and dimensions), and so although our main focus in this paper is three-dimensional field theory at non-zero temperature, we will provide original results for the zero temperature SU(5) in four dimensions elsewhere; see also [28] for the associated entanglement c-function that we have amused on already.
Armed with the success in four dimensions, we continue the discussion in three-dimensional case. Although we are limiting our focus on SU(2) theory, since the simulations are costly, we will consider a novel scenario and consider entanglement entropies at nonzero temperature. In addition to for the first time to extract entanglement entropies in three dimensions, there are interesting features that we get to test. First of all, contrary to four dimensions, there is no UV fixed point, so one does not get to exploit the underlying conformal symmetry to claim success of the comparison to holographic result. Interestingly, we will find support that the gravity dual to supersymmetric Yang-Mills theory at strong coupling and in the limit of large-N c elucidates many features that would be otherwise unexpected from field theory point of view.
So, what is the gravity dual of the Euclidean three-dimensional pure Yang-Mills theory that we consider? In an optimal scenario we should really be considering a supersymmetric Yang-Mills theory on the lattice, e.g., as in [29]. Since the results that we obtain are pretty encouraging, we believe that our work is a first step towards this direction and similar methods that we construct could be considered in vastly more demanding cases with fermions. We will review salient details of the holographic approach in Sec. 3.
The key observation made in [30] is that the derivative (1.3) can be directly and very effectively used to reconstruct selected components of the metric in the dual geometry to an accuracy given by the statistical error in the measurements of (1.3). This boils down to a few facts, that ∂S EE /∂ is independent of a regularization scheme and that it can be written without an explicit integral that would otherwise be present in the evaluation of the entanglement functional S EE . The derivative is instead written in terms of components of the background metric, evaluated on the position of the RT surface. This quantity can be understood as the conjugate momentum that is conserved under translations of the position of the strip (slab) [31] and it becomes a very simple expression when one chooses to evaluate ∂S EE /∂ at the tip of the RT surface. In particular, in the limit → ∞ one expects to probe the IR physics of the ambient field theory. On general grounds in the deconfining phase, the finite part S EE / should pick only the contribution from the degrees of freedom on the black hole horizon and indeed the expression ∂S EE /∂ | →∞ reduces to the Bekenstein-Hawking entropy S BH ∝ N 2 c T 7/3 /λ 1/3 , where λ = g 2 Y M N c indicates the 't Hooft coupling [32,33]. Thereby, per holographic dictionary, extracting the RHS of (1.3) will match onto the thermal entropy of Yang-Mills theory at strong coupling and large-N c when → ∞ is taken. We will show numerical evidence for this from our lattice simulations in Sec. 5. That the derivative of the entanglement entropy of some region A becomes in the limit → ∞ proportional to the thermal entropy, S th , can be understood also directly from the definition of the entanglement entropy, Eq. (1.1), by noting that in the limiting case, where A represents the whole system, and its complement, B, shrinks to zero volume, one has ρ A → ρ, and therefore: lim B→∅ S EE (A) = tr(ρ log(ρ)) = S th,A . (1.4) Now, as the thermal entropy is an extensive quantity, S th,A grows linearly with the volume of A, which in turn depends linearly on . One would therefore expect that with s th,A being the density of thermal entropy in A, and ∂ vol(A) = |A ⊥ | is simply the area of the cross section of A that is perpendicular to the direction in which grows. We will discuss this in more detail in Sec. 4. The relation Eq. (1.5) could also be of interest for future thermodynamic studies, as it allows for a non-perturbative determination of the thermal entropy of a system without having to fix integration constants. At intermediate energy scales (or strip widths) we find evidence for the scaling S EE ∼ −4/3 that is expected from (supersymmetric) YM theory at zero temperature in three dimensions [34]. We present the results that we have obtained from lattice studies in Sec. 5. We also show that the corresponding geometry can be reconstructed using this data. Although the reconstructed metric might give further information than the entropy density on some field theory quantities, and possibly even make useful predictions, we defer these investigations in the future.
The rest of the paper is organized as follows. Sec. 6 contains our conclusions and a discussion on future outgrowths of our work. Appendix A contains a review on selected topics in lattice gauge theory for ease of reference. We also supplement the article with an Appendix B where we give a peek to other novel results that we extracted from the lattice in the case of 3d YM at finite temperature. We present data for the Polyakov loop and extract the potential V (L) between quarks and anti-quarks, separated by a distance L. Among other things, in the intermediate energy scales, we show evidence for the scaling V (L) ∼ L −2/3 [35] stemming from the same D2-brane background in the UV as we utilized for the entanglement entropy. At the deep IR, i.e., for very large separation, the D2-brane background suggests that the real part of the potential ReV ∼ L −10/3 , while the imaginary part would grow as ImV ∼ L. Our results are consistent with the former behavior though we cannot dissect it from the expectation of the Debye screening behavior below T c : ReV ∼ e −m D L /L 1/2 where m D would be the Debye mass [36]. 2 To this end, a systematic scan in the parameters is required, in particular to analytically continue the Euclidean correlators to real time [38,39] to see how the imaginary part of the potential behaves.

Entanglement entropy on the lattice
In this section we describe in some more detail how the aforementioned replica method [21][22][23] can be applied to determine entanglement measures in SU(N c ) Yang-Mills theories using non-perturbative lattice Monte Carlo methods. For the reader who is unfamiliar with the lattice formulation of Yang-Mills theory and Monte Carlo techniques, we introduce in Appendix A the basic concepts.
The first implementation of the replica method for SU(N c ) lattice gauge theory was introduced almost 15 years ago [4] and used ever since [6][7][8][9]. We will here briefly review the working principle of the original method and sketch the considerations that lead to our new approach.
In Appendix C we provide some more details on our update algorithm from which also an argument arises for why in (1 + 1) dimensions the entanglement entropy must be independent of the width, , of the slab region A.

Established method
The replica method is based on a path-integral representation of matrix elements, ψ 1 |ρ|ψ 2 , of the density matrix ρ, which can be thought of as the (Euclidean) time-evolution operator from some initial time x d 1 to some final time x d 2 , and ψ 1 and ψ 2 represent instantaneous states of the system at these times. For an SU(N c ) gauge theory on a d-dimensional finite lattice of temporal extent N t , the initial and final states ψ 1 and ψ 2 can be thought of as defining values (up to gauge transformations) for the gauge links that are within or touch the time slices x d = 0 and x d = N t . Explicitly, where S G is the Wilson gauge action (A.18) and U ψ are made up from link-variables (A.10). The normalization factor Z can be obtained by noting that taking the trace over ρ means to identify the two boundaries at x d = 0 and x d = N t with each other and sums over all possible values. From the normalization condition one then finds that is the lattice partition function for the given system in the case of periodic temporal boundary conditions. Note that in this case, a finite temporal extent, N t , of the lattice corresponds to a finite temperature T = 1/(a N t ).
By splitting the lattice system into two parts, A and B, as depicted in Fig. 2, matrix elements for the reduced density matrix of part A can then be represented as where only the temporal boundary states, ψ B,1 and ψ B,2 for part B get identified with each other and summed over all possible values, indicated in the drawing after the last equality sign by assigning the label r B to both temporal boundaries. The temporal boundaries for part A are specified by the boundary states ψ A,1 , ψ A,2 . Figure 2: Illustration in (2 + 1) dimensions of how we divide the system with fixed temporal boundary states, ψ 1 and ψ 2 , into two parts, A and B. After the partitioning, we have ψ 1 = ψ B,1 ⊗ ψ A,1 and ψ 2 = ψ B,2 ⊗ ψ A,2 , where ψ B,1 , ψ B,2 represent the boundary states for part B and ψ A,1 , ψ A,2 the boundary states for part A.
With these definitions, it is now straightforward to obtain an expression for the so-called Rényi entropy of order s ∈ N in terms of a ratio of two lattice partitions that represent tr(ρ s A ): Here Z c ( , s, N t , N s ) is the partition function for a system with the topology depicted on the left-hand side of Fig. 3. The example shows the case of s = 2 in (2+1) dimensions with region A having width = 2. The right hand side of Fig. 3 shows the same lattice but with a different topology, so that it decomposes into a product of two copies of the system described by the normal partition function Z(N t , N s ). Using the expression for entanglement entropy given in Eq. (1.2), one finds with these definitions: where on the last line we approximated with ∆s = 1, and used that for s = 1 one has F c ( , 1, N t , N s ) = F (N t , N s ). With this approximation, the entanglement entropy turns out to be just the Rényi entropy of order two: Figure 3: Illustration of two different topologies for a (2 + 1) dimensional lattice of size V = N 2 s × s N t with s = 2. The left-hand side shows the topology for the system described by the partition function Z c ( , s, N t , N s ) in Eq. (2.5), and the right-hand side the system described by Z s (N t , N s ), i.e., the product of s copies of the system described by the ordinary lattice partition function Z(N t , N s ).
As the entanglement entropy has a -independent UV-divergence, it is convenient to look at the derivative of S EE with respect to instead of S EE itself. On the lattice, this corresponds again to a discrete lattice derivative, Apart from getting rid of UV divergences, Eq. (2.9) has from lattice Monte Carlo perspective also the advantage over Eq. (2.6), that the change that has to be done to the system to increase the width of region A from to ( + 1) is in general much smaller than the required change to get from width = 0 to . However, also the change from to ( + 1) represents a highly non-local change that affects a large number of degrees of freedom. As a consequence there is a bad overlap problem, meaning that link-variable configurations that contribute to the partition function of the system play essentially no role for the partition function of the ( + 1) system, and vice versa.
To overcome the overlap problem, the authors of [4] proposed to define a one-parameter family of interpolating partition functions, Z * (α), with α ∈ [0, 1], so that where we leave the dependency of Z * on N t and N s implicit in order to avoid too lengthy expressions. The interpolating partition function is chosen to be where S G, [U ] and S G, +1 [U ] represent for a given configuration of link-variables, U = {U x,ν } x,ν , the two values of the gauge action where the subsystem A has width and ( + 1), respectively. The derivative of the entanglement entropy with respect to from Eq. (2.9) is then obtained as where the integration over α is performed by interpolating lattice results for the expectation values for a sufficiently dense set of α-values in the integration domain [0, 1]. An example for a possible outcome of this procedure is shown in Fig. 4, using data for a SU(3) gauge theory on a (3 + 1)-dimensional lattice with N s = 16, N t = 16, s = 2 and = 2 at β g = 6.0, which has been extracted from [6]. The left-hand panel of Fig. 4 shows the interpolating function for Eq. (2.13) as function of α and the right-hand panel shows the corresponding running free energy differences, where ∆F (1) corresponds to the value of Eq. (2.12) we are interested in. The example depicted in Fig. 4 illustrates the problem that arises when using the just described method for interpolating between the partition functions (2.10). During the interpolation process, the free energy has to pass through a huge maximum. This is presumably due to the fact that the configurations that contribute to the interpolating partition function from Eq. (2.11) at 0 < α < 1 are simultaneously subject to both actions, the one corresponding to a width of region A, and the one for region A having width ( + 1). As the configuration distributions sampled with either of two actions separately have very little overlap (overlap problem), the number of configurations that can give a significant contribution to the interpolation partition function, Eq. (2.11), must be highly reduced when α ∼ 1/2, resulting in a significantly increased free energy. This is problematic because the values of Eq. (2.13) for different α values are determined with Monte Carlo methods and therefore come with a statistical error. Upon integration over α, the positive and negative values of Eq. (2.13) that occur as α runs through the interval [0, 1] lead to large cancelations, so that that total free energy difference, ∆F (1), is again relatively small, the errors in the measured values of Eq. (2.13) cannot cancel but simply accumulate, giving rise to a bad signal-to-noise ratio.

Improved method
To avoid the formation of huge free energy barriers and resulting bad signal-to-noise ratios, we can use an alternative interpolation method. To this end, let us denote by C the set of all plaquettes in a V = s N t × N d−1 s lattice, which are not affected by a change of temporal boundary conditions. These are all plaquettes except for the temporal ones that touch from below the time-slices with x d = r · N t for r = 1, . . . , s (cf. Fig. 5).
A generalized partition function can then be written as: wherex = x 1 , . . . , x d−1 ∈ Z d−1 labels the spatial positions on the lattice and U νd (x) and x y A B Figure 5: The figure illustrates at the example of a (2+1)-dimensional lattice of size N 2 s × s N t with s = 2, how a change of temporal boundary conditions, associated to a change of the width of region A from = 2 to = 3, affects the way in which plaquettes are computed. The plaquette P 1 in the top-left panel touches from below the dashed line r B 1 , meaning that the top-edge of P 1 does not consist of the link variable that connects x 1 and x 2 , but rather the link variable that connects x 1 and x 2 . Similarly, the upper edge of plaquette P 2 does not consist of the link between x 1 and x 2 but of the link between x 1 and x 2 . In the top-right panel, the same temporal plaquettes are marked in red and belong now to region A: this means that the link that closes the plaquette P 1 along its top-edge is now indeed the link variable that connects x 1 and x 2 and the link that closes P 2 along its top-edge is the one that connects x 1 to x 2 . The two lower panels show for = 2 (left) and = 3 (right) the spatial links over which the temporal plaquettes are subject to the boundary conditions of region A (red) and region B (blue).
U (1) νd (x) are the two different values for the temporal plaquette depending on whether the plaquette is subject to the boundary conditions from outside or inside region A (see Fig. 5: (0) corresponds to blue, (1) to red). Which of these two cases applies is for each spatial link controlled by the value of a corresponding discrete variable the generalized partition function (2.15) reduced to Z c ( , s, N t , N s ) as introduced in the previous section.
To describe how the interpolation between Z c ( , s, N t , N s ) and Z c ( + 1, s, N t , N s ) can be carried, using the generalized partition function Eq. (2.15), let us denote by ..,j be the subset of the first j elements of K. Let us then further define: and abbreviate Z i = Z(β, s, N t , N s , {n i }) and F i = − log(Z i ). The expression for the derivative of the entanglement entropy (2.9) then becomes: This free energy difference can be measured from a single simulation, using a form of multi-canonical method [40] for discrete sets of states, known as Wang-Landau (WL) sampling [41]. With this method one samples the modified partition function,  Possible choices for ordering the spatial boundary link variables in K are illustrated in Figs. 6 and 7 for the (2+1)-dimensional and the (3+1)-dimensional case, respectively. These examples exploit the fact that gauge invariance ensures that a change of temporal boundary conditions over spatial links for which at least one end is always (before and after the update) either completely in region A or completely in region B, does not change the free energy (cf. Appendix C). This is the reason why the free energy graphs in Figs. 6 and 7 are flat at small i, where i is used to enumerate the boundary link variables. To determine the total free energy difference in Eq. (2.23), one therefore has to consider only the boundary states for i > N d−2 s . Note also, that with this choice of ordering of the spatial links in K, the change in free energy as function of i is piecewise linear. This will be particularly useful for future studies on larger systems, as it means that one does not have to measure the change in free energy for all N K boundary states, but only has to determine the slope of ∆F i as a function of i in each of the piecewise linear interval.  After all spatial links perpendicular to the boundary of region A have been added, one continues to add spatial so as to form complete cubes that fill one row after the other. Note that the sketches on the right-hand side show a much smaller lattice than the one used to obtain the data in the left-hand plots. The sketches are meant to illustrate why the free-energy difference as a function of n changes slope for a certain points. The total free energy difference for the here discussed example is the same as in Fig. 4 but is obtained without having to overcome a equally huge free energy barrier.

Holographic approach
Entanglement entropy is a difficult quantity to compute in general QFTs and even if the theory admits a lattice prescription the computation proceeds nontrivially as we have seen in the previous section. However, if a QFT admits a holographic description, its computation becomes easy using the Ryu-Takayanagi (RT) formula [20,42,43]. The RT-formula states that the entanglement entropy is given by the area of a certain minimal bulk surface in the holographic dual of the QFT. In this section we will review salient details on the computation of the derivative of entanglement entropy. In particular, we will show how the derivative can be used to solve the inverse problem, to infer the dual background metric, which can only be obtained up to statistical error inherent in the data of the derivative.

Derivative of the entanglement entropy
In a holographic theory the subsystems A and B are regions on the conformal boundary of the bulk dual theory. The QFT is said to live on this conformal boundary while the dual gravity theory lives in the bulk enclosed by the boundary. The boundary QFT and bulk gravity descriptions are dual in the sense that they encode the same physical information. The two descriptions are equivalent but some quantities might be much easier to compute on one side of the correspondence than on the other. A great example is the entanglement entropy which is given by where Σ A is the minimal area surface associated to the boundary region A and G (10) Nc is the ten-dimensional bulk Newton constant which can be written explicitly in terms of gauge theory quantities if the dual theory is known. We consider the dual to be Yang-Mills theory and imagine 1/G (10) Nc ∝ N 2 c . Here, we are assuming that the external bulk geometry is (d + 1)-dimensional, i.e., that the boundary QFT is d-dimensional. The 8-dimensional bulk Figure 8: The holographic bulk surface used in computing S EE (A) using the RT-prescription. The QFT lives on the boundary B d = A ∪ B, and the dual gravity theory lives in the bulk M d+1 (times the internal geometry). The bulk codimension-2 surface Σ A is anchored to the boundary of A so that ∂Σ A = ∂A. The bulk metric in M d+1 is such that it causes Σ A to hang in the bulk when its area is being minimized.
surface Σ A has to be anchored to the boundary region A whose entanglement entropy we wish to compute, that is, ∂Σ A = ∂A. Additionally, Σ A has to be homologous to A in the bulk meaning that one has to be able to retract Σ A smoothly close to A through the bulk. There is an infinite family of bulk surfaces that satisfy these conditions. The RTformula then instructs to pick the surface with minimal area as measured by the bulk metric tensor. The area of this minimal surface then gives S EE (A). Note that since the region A is defined on a boundary time slice, the bulk surface Σ A is of codimension two. Also note that we are working with an asymptotically AdS d+1 bulk theory and therefore the area of Σ A is divergent because it has to reach all the way to the boundary of the spacetime. This divergence dominates the entanglement entropy, and this implies that the holographic formula follows the area law. This is the RT-prescription for static cases which is the case of interest to us in the present work, but there are extensions of the prescription applicable to more general situations [44].
In this paper, we consider slab-like entangling regions A. That is, regions defined by where is the width of the slab. We will always orient the slab orthogonal to the first spatial boundary dimension x 1 . The slab is assumed to be large in directions other than x 1 . This slab shape combined with the translation invariance of the theory implies that the profile of the bulk surface can be represented by a single function where z is the holographic coordinate orthogonal to the boundary. We are interested in the entanglement of such slabs in (3 + 1) and (2 + 1) dimensions. We will start with the former case. Consider a bulk metric in the Einstein frame of the form where b(z) = 1 − z 4 /z 4 h and R is the radius of curvature. This metric describes a family of asymptotically AdS 5 spacetimes if a(0) = 1. The holographic coordinate z is such that z = 0 corresponds to the conformal boundary. The function b(z) is the usual blackening factor of the AdS planar black hole, which puts our theory to a finite temperature. Further, we will only consider functions a(z) such that a(z h ) = 1, where z h is the black hole horizon position. It follows that the dual QFT is at temperature It is worth pointing out that since we consider a static situation, the entanglement entropy actually does not depend on the metric component g tt at all. We chose to fix the relationship between g tt and g zz so that later we will be able to infer the complete bulk metric. 3 The entanglement entropy is then given by the integral , V is the regularized volume of the slab in directions x 2 and x 3 , and we have simply integrated over the internal space so that its volume is absorbed into the definition of the five-dimensional Newton constant G (10) Nc R 5 Vol(Ω 5 ). The point z = z * is called the turning point of the bulk surface, which is the deepest point in the bulk reached by the minimal surface. The profile of the surface x 1 (z) has to be such that its area is minimized. This happens when x 1 (z) satisfies the Euler-Lagrange equation with the Lagrangian L given by the integrand of (3.5). Since x 1 (z) itself does not appear in the integrand, we immediately find a conserved quantity ∂ x 1 L. We can fix the value of this constant at the turning point z = z * where x 1 (z) → ±∞ and find an equation for the profile which minimizes the area: where the ± refers to the two different branches of Σ A as a function of z. After we now know the profile, we can give explicit integrals for the slab width and the entanglement entropy S EE : On the last line we have separated the UV-divergence from the z = limit of the integral. It is instructive to study the behavior of S EE ( ) for small , that is, in the UV-limit. This limit is easily obtained from the integrals (3.8) and (3.10). The UV-limit is equivalent to assigning a(z) = 1 and b(z) = 1, which corresponds to the zero-temperature N = 4 supersymmetric Yang-Mills theory. In this limit the integrals can be computed and one finds where is the UV-cutoff we used in the integral (3.10). Note that here we displayed the canonical field theory example of N = 4 SU(N c ) super Yang-Mills theory on flat space, dual to AdS 5 × S 5 , wherein we plugged in known expressions for the parameters in terms of field theory quantities using Vol(S 5 ) = π 3 , R 4 = 4πg s α 2 N c and G (10) Nc = 8π 3 α 4 g 2 s . The first term reflects the expected UV-divergence of S EE ( ). On the QFT side, it originates from short-range correlations across the entangling surface. On the gravity side, the divergence has a geometric interpretation, as mentioned already above: since the surface Σ A has to reach all the way to the boundary, its area is infinite because the boundary is infinitely far away. Note that S EE is indeed proportional to V , the volume of directions parallel to the slab: this is the area law of entanglement entropy.
The diverging term is not particularly interesting since it is independent of the slab width . The occurence of this divergence is actually related to the fact that the entanglement entropy itself is not an observable, similarly to say the total free energy of a system, which is typically infinite. In order to bypass this problem, one can focuse on the finite piece in (3.11) or on its change. To this end, it is easy to show that the derivative of S EE ( ) takes the following simple form [30]: where we recall that the width is obtained by inverting z * = z * ( ). The expression (3.12) is particularly nice since all integrals cancel each other out and the right-hand side reduces to a simple function of the metric coefficients. If we were to reconstruct the bulk geometry using lattice measurements of entanglement entropy in (3+1) dimensions, we could use the formula (3.12) very effectively. Notice further that when considering the infinite width limit → ∞, the tip of the RT surface approaches z * → z h , and the only contribution to the finite piece stems from (3.12), which gives the thermal entropy where in the latter expression we again wrote the result for N = 4 SU(N c ) super Yang-Mills theory on R 4 at temperature T .
Having reviewed a few key holographic ideas we will next introduce the (2+1)-dimensional model we utilize. We use the following bulk geometry in the string frame: (3.14) where z p is a length scale playing the role of the radius of curvature and e φ is the dilaton. Notice that the our ansatz for the dilaton is not crucial; using Weyl rescaling we could aim for reconstructing in the Einstein frame instead. We stick to the string frame to allow for easier comparison with results in the string literature. The blackening factor is b(z) where z h is the horizon position. This background for a(z) = 1 can be derived from type IIB string theory by considering a stack of D2-branes [45]. This is a top-down holographic model, so the metric includes a compact internal space; we keep the ten-dimensional Newton constant explicit below. The temperature in this family of geometries, assuming again a(z h ) = 1, is The integrals for the strip width and entanglement entropy can be derived in the same way as in the (3 + 1)-dimensional case. The relevant integrals are Like before, the derivative of S EE ( ) is a simple function, free of integral expressions: where if interested in this relation in terms of ten-dimensional Newton constant one uses Vol(S 6 ) = 16π 3 /15 for the volume of the internal space. We can work out the UV limit of S EE by studying (3.17) and (3.18) in the limit z * z 0 . The result is [34], We note that in order to avoid flowing to the weakly coupled Yang-Mills regime, we need to work with finite UV cutoff , however, we are only interested in terms depending on in this work.
On the other hand, in the large width limit Here one should also be wary that the strict IR limit is probably not trustworthy at face value, but better captured using 3D superconformal field theory via M2-branes [32,42]. Our focus is on the middle term in the bracket in (3.21) (or actually that in (3.19)) that arises solely from black hole horizon and corresponds to the thermal behavior. We believe that it universally describes the thermal entropy of the deconfining phase in the field theory. Indeed, below we show evidence for the S EE ∼ T 7/3 behavior as extracted from lattice computations.

Bulk reconstruction
We now have expressions for the holographic entanglement entropy in both (3 + 1) and (2 + 1) dimensions. In this work, however, we only consider the reconstruction in (2 + 1) dimensions since this is the case for which we have sufficient data from the lattice. The formulas require as input the bulk metric and output entropies. Since we have also computed the derivative of S EE ( ) with lattice simulations, it is interesting to ask whether it is possible to find such a bulk geometry for which the holographic predictions would agree with the lattice results. Schemes like this where the holographic bulk geometry is inferred from QFT quantities is called bulk reconstruction. We will use statistical methods for reconstructing the bulk from lattice data [30]. The lattice data consists of a set of M measurements of dS EE /d and corresponding uncertainties, σ, for different widths : The bulk geometry contains only one free function, a(z), which we parametrize as The expressivity of the ansatz is controlled by the number of free parameters, N basis . Furthermore, the ansatz is built so that a(0) = a(z h ) = 1 by construction, which keeps the geometry asymptotically conformally AdS and the horizon at a constant temperature. We note that if we were interested in very low temperatures and the confining phases of the theory, a metric ansatz with a collapsing cycle at the IR end of the geometry (a 'cigar') would seem more appropriate.
Our statistical model then has the parameters Nc ) is the multiplicative constant on the right hand side of (3.19). We will link strip widths i to entropy derivatives by assuming a likelihood We will use weakly informative normal priors for a and c with standard deviation 5 and 1, respectively, around their maximum likelihood estimates. Finally, we sample the posterior distribution of a and c with Hamiltonian Monte Carlo. The parameter distributions can then be used to find the distribution of bulk metrics consistent with the lattice entanglement entropy data. In a sense this is like using the AdS/CFT duality "in reverse". Then, having access to the bulk gravity model, we could use the AdS/CFT dictionary in the usual direction to derive other quantities of interest, for example, quark-antiquark potentials.

Entanglement and thermal entropy
In this section we would like to show that which tell us that at sufficiently large width of the entangling region A, the entanglement entropy, S EE ( , T ), of A reduces to the thermal entropy, S th,A (T ), in A. Note that S th,A (T ) is an extensive quantity and grows therefore linearly with the volume of A( ), i.e. linearly with . A more straightforward way to write this is, to show that the derivative S EE with respect to , becomes -independent for large : where |A ⊥ | is the -independent area of a cross section of the entangling region A perpendicular to , and s th,A is the thermal entropy density for region A.
As mentioned in the introduction, in terms of density matrices, one can see that the relation in Eq. (4.1) should hold, by noting that the entanglement entropy in region A is defined by: where ρ A (T ) = tr B (ρ(T )) is the reduced density matrix for region A, obtained by taking the partial trace of the full density matrix, ρ(T ), with respect to a basis for the complement, B = A c , of region A. In the limiting case, where A represents the whole system, and B is empty, one has ρ A (T ) → ρ(T ) and therefore i.e. the expression for the entanglement entropy reduces to the expression for the thermal entropy if the region A is entangled with nothing. This makes it plausible that Eq. (4.1) holds as soon as A is much bigger than the thermal screening radius, r T ∼ T −1 , as then most of A will not "feel" that B has been traced out; only the part of A that is within a distance r T from the boundary to B can experience entanglement with region B.
In the following, we will provide an additional argument for why Eq. (4.1) resp. Eq. (4.2) should hold, based on a comparison of the replica trick expression for the entanglement entropy and the expression for the thermal entropy in terms of a Euclidean path integral. In order to have well-defined path integral expressions, we will consider directly lattice regularized theories. However, the provided relations should apply also when using other valid regularization schemes.

Thermodynamics in terms of a Euclidean (lattice) field theory partition function
Let us define the lattice free energy, F L (N t , V, N ) as where Z(N t , V, N ) is a canonical, Euclidean lattice partition function, depending on the temporal lattice size, a N t , the spatial lattice volume, a 3 V , and a conserved charge N . We keep the dependency on the lattice spacing, a, implicit as we are interested in the properties of a lattice system at fixed a. If the theory under consideration does not contain any conserved charges, one can drop the dependency on N ; similarly, if there are multiple conserved charges, one will have to extend the formalism to an appropriate set of charges, The pure gauge system we study in this paper does not contain any conserved charges in which are interested, but we keep the N -dependency for the moment, in order to emphasize the relation between F L and the usual Helmholtz free energy F (T, V, N ), which is given by From Eq. (4.6) it follows that the differential of F L in therms of the usual thermodynamic quantities, T = 1/N t (temperature), p (pressure), µ (chemical potential), S (entropy), V (spatial volume) and N (conserved charge number), is given by: where after the second equality sign we have used that the internal energy U is related to F by and that In order to express the thermal entropy S in terms of F L , we can use Eqn. (4.8) and (4.7) to write:

Relation between entanglement and thermal entropy
To demonstrate the relation from Eq. (4.1), we note that in order to determine Eq. (2.9), one does not necessarily have to implement the change → ( + 1) in the width of the region A by changing the temporal boundary conditions as described in Fig. 5 along the current boundary of region A; one could just as well increase by removing a slice of width ∆ from deep inside region B, change the temporal boundary conditions of that slice, and insert it again somewhere in the interior of region A, as illustrated in Fig. 9.
Let us now assume that the width of region A and the width (N s − ) of region B are both sufficiently large, so that the effect of the interfaces between the two regions, located at x = 0 and x = N s − , can be neglected for describing the thermodynamic properties of the system deep inside of either of the two regions. Similarly we assume that the change in the volumes of the two regions, A and B, when → ( + 1) can be neglected. In this case, deep inside region B, the system consists of two copies of a system that can be described by the lattice free energy density whereas deep inside region A, the system consists of just one copy of a system that is described by a lattice free energy density In terms of these lattice free energy densities, the derivative of the entanglement entropy from Eq. (2.9) is given by where N (d−2) s is the spatial area (in lattice units) perpendicular to the x-direction, so that N (d−2) s ∆ with ∆ = 1 is the spatial volume of the slice that is moved from region B to region A. Now, as discussed around Eq. (2.8), the lattice expression for S EE actually corresponds to the second order Rényi entropy, defined in Eq. (2.5). This is a consequence of approximating in Eq. (2.6) the derivative with respect to the replica number, s, in the limit (s → 1) by a discrete forward derivative: 14) The expression on the right-hand side of Eq. (4.13) is therefore a consequence of this forward derivative approximation in (4.14). If one tries to reverse-engineer what the right-hand side of Eq. (4.13) would look like if the discrete derivative approximation in Eq. (4.14) were not necessary, one quickly finds as possible solution that: in which case one would have: with s th,A being the density of thermal entropy, S th,A , of region A, as follows from Eq. (4.10). With Eq. (4.16), one could therefore compute the thermal entropy of a quantum field theory directly from a single lattice Monte Carlo simulation, using the method for determining the derivative of the entanglement entropy, described in Sec. 2.2. This is not possible by using only Eq. (4.10), as this would involve knowledge of lattice free energy, F L , itself. Instead one would have to measure for example ∂ T S(T ), where T = 1/N t , for multiple values of T , that will allow one to numerically approximate the integral in and where the integration constant, S 0 = S(T = 0) = 0, is fixed by the 3rd law of thermodynamics. This latter procedure appears potentially very expensive, as one has to simulate low temperatures, which means large N t and therefore large lattices that are expensive to simulate. The relation Eq. (4.16) might therefore in the future also be useful to study the thermodynamic properties of strongly interacting field theories.

Results
In this section we will discuss results that we have obtained from running simulations on the lattice as described in Sec. 2.2. A comparison between results obtained with our new method and corresponding literature results for the case of four-dimensional Yang-Mills theory has been briefly exposed in [28,46]. In the following we will present our results on the entanglement entropy in (2+1)-dimensional pure SU(2) gauge theory at high temperatures. Along the way we recall how the holographic approach, as described in Sec. 3, guides us, e.g., in picking correct power laws for fitting. We provide an example for how the replica trick, used to determine the entanglement entropy, can be utilized to determine also the thermal entropy, resorting to lim →∞ ∂ S EE ( ) ∝ S thermal , as discussed in Sec. 4. In the final subsection we reconstruct the bulk geometry from the entanglement entropy measurements.

Three-dimensional Yang-Mills theory at high temperature
We are interested in lattice results on the behavior of the derivative of the entanglement entropy with respect to the slab-width, , at high temperatures (T > T c ) and large slab-width ( > T −1 ), where T is the temperature and T c the critical temperature for deconfinement, which serves as physical reference energy scale.
We have performed simulations at five different values of the lattice spacing, a, and various different temperatures. The list of simulated setups is given in Table 1. Note that due to our interest in the high-temperature regime, we are forced to go to rather large values of the inverse gauge coupling, β g . As discussed in [32], this means that the criteria for the applicability of the classical approximation for the bulk theory in which the holographic computations in Sec. 3 are carried out, are satisfied only at relatively large values of the slab-width , compared to the lattice spacing a. The range of validity as function of β g and in lattice units is sketched in Fig. 10 for N c = 2, using the inequalities from [32]: where g 2 c is the continuum theory gauge coupling. In terms of lattice quantities, using the relations β g = 2 N c /g 2 0 , g 2 0 = a g 2 c , and c = a , the inequalities take the form: Our lattice data for the derivative of the entanglement entropy might therefore at short distances not conform with the power law behavior, expected from the discussion in Sec. 3. For our simulations at β g = 16, also the large-plateau value is for the lowest temperatures N t = 8, 10 reached only at values of that are already outside the range of validity, Eq. (5.2).
With this in mind, we allow for certain corrections to the expected behavior discussed in Sec. 3 when performing fits to our lattice data. An example is shown in Fig. 11. From holography (cf. Eq. (3.20)), we would expect our lattice data for dS EE /d as function of to behave in the UV like a power law ∼ −7/3 . As our simulations are carried out at temperatures, T , which are 2-5 times the critical temperature T c , this power law behavior is almost completely screened by the plateau within the range of validity for the holographic predictions from Eq. (5.2).
We will therefore focus on verifying that for T /T c 1 one has  Figure 11: Example of fitting the lattice equivalent of dS EE /d as a function of for SU(2) in (2+1) dimensions at three different temperatures T ∼ 1/N t with N t =8,12,16 for a spatial lattice size of V = 128 2 . As T /T c > 1, the derivative of the entanglement entropy with respect to approaches at large a temperature-dependent plateau-value. The length unit Λ cr corresponds to the inverse critical temperature, i.e., Λ cr = 1/T c . The faint vertical lines represent the lower and upper limit of the holography-validity interval from Eq. (5.2). Note that the data, corresponding to different N tvalues (i.e. different temperatures), does not yet lay on a unique curve at short distances, indicating that our lattice data does not resolve a sufficiently high UV-regime at which the temperature would no longer matter. Furthermore, we note that the short distance data lays outside the holography-validity range and does not follow the expected pure power law ∼ x −7/3 (black line); we had to allow for a shift, x 0 , in order to fit the data (red lines).
(T /T c ) 7/3 for T /T c 1. As a first test of this scaling law, we show in the first three panels of Fig. 12 our lattice results for ∂ S EE ( ), rescaled by (T /T c ) −7/3 and as function of ( T ) for β g = 32 (top), 48 (second from top) and 64 (third from top). For each value of β g the plots consist of data from various temperatures. Thanks to the rescaling of ∂ S EE ( ) by (T /T c ) −7/3 , the plateau data corresponding to different temperatures coincides. This is in particular true for β g = 48, 64, corresponding to lattice spacing values a T c ≈ 0.024, 0.032. For β g = 32, resp. a T c = 0.048 the data does not collapse as nicely anymore, and for the even lower β g -values, resp. larger lattice spacings, it would look even worse. The bottom panel in Fig. 12 shows a linear extrapolation of the rescaled plateau values, which provides a first rough estimate of the continuum value of the proportionality factor, s th,0 in Eq. (5.3).
In order to extract the plateau values for lim →∞ ∂ S EE ( ) more accurately, we can try to fit our lattice data for each value of T /T c separately. As mentioned at the beginning of this section, we have to expect deviations from the ∼ −7/3 power law behavior and we therefore modify the powerlaw by allowing for a shift not just in y-, but also in x-direction.   We therefore fit the ansatz, From these fits, we obtain the plateau-value, b, as function of T /T c = T Λ c and a T c . Here a is the lattice spacing and Λ c = 1/T c is the deconfinement length scale. The ratio T /T c as function of β g , N s , and N t is determined as described in Appendix B.2.3, and a T c = 1/(N t (T /T c )). The so obtained plateau values, b, are plotted in Fig. 13 as function of the reduced temperature τ = T /T c − 1. Superimposed to the data is for each value of β g separately a fit of the form y = a 0 (1 + x) 7/3 1 + a 1 (1 + x) −1 , with (x, y) ∈ {(τ, b(τ )) | ∀τ }. As before, this fit ansatz is motivated by the identification of the plateau value with the thermal entropy density. The latter should, according to holography, for T −1 behave like s thermal ∝ (T /T c ) 7/3 = (1 + τ ) 7/3 , and we include a leading correction term, whose relative magnitude is parametrized by a 1 . The fits are performed on data for which T /T c > 1.9, as lower temperatures would require higher order corrections. The only exception is the case of β g = 16, where T /T c > 1.7 is used as lower bound in order to have a minimum of three data points to fit the two parameters a 0 and a 1 . However, as can be seen, the quality of fit for β g = 16 is poor. This is on the one hand due to the extended lower bound of the fit-range which, as mentioned before, would require higher order corrections, and on the other hand due to the fact that the lattice is too coarse at β g = 16, so that N t = 2 for the highest temperature, T /T c ≈ 5.12. We therefore exclude the β g = 16 data from further analysis. T/T c =4.8 The data from the remaining β g -values is used to perform a continuum extrapolation (β g → ∞) in two different ways (black and cyan bands).
The first extrapolation (black) is obtained using the fitted functions and corresponding error bands to obtain values and errors for b at fixed T /T c but different lattice spacing values, a T c , and perform a fit:  Fig. 13. The second continuum extrapolation in Fig. 13 (cyan) is obtained by considering directly the lattice spacing-dependency of the fit results for the parameters a 0 and a 1 for β g = 24, 32, 48, 64. As shown in Fig. 15, we perform quadratic fits to the data for a 0 and a 1 as functions of a T c . The so obtained fit results for the constants a 0,0 and a 1,0 and their errors are then used to draw the cyan curve and error band in Fig. 13. As can be seen, the two continuum extrapolations are in good agreement.

Reconstructed bulk geometries
In Section 3 we have discussed our statistical model which we use to connect the metric tensor of the holographic dual to entanglement entropy values. We will now use data of the derivative of entanglement entropy with respect to the slab width obtained from the lattice to find which metric fits these observations. More specifically, the data is for SU(2) lattice gauge theory on a V = 128 2 lattice at temperature T /T c = 2.646. Since the the gauge theory is in (2 + 1) dimensions, the appropriate holographic background is the asymptotically D2brane geometry.
We are sampling the posterior distribution of the parameters which define the metric tensor to find metrics which would fit the lattice data well under the assumption that the holographic model represents the gauge theory well. The results are shown in Fig. 16. There is a pronounced bump in the metric towards the UV-region of the geometry. It would be interesting to analyze the ramifications for this feature on some derived quantities. We have also plotted the predictions for dS EE /d using the maximum likelihood holographic geometry and find good agreement with the lattice gauge theory data.

Discussion and outlook
We demonstrated that holography has entered a new phase, precision science. Gravity duals of realistic field theories can be reconstructed from underlying data that is always bound to be imprecise, a fact which subsequently feeds in likelihood estimates on any holographic computations. Systematically reducing the statistical uncertainties on the extracted quantities from the lattice leads to an increasingly accurate metric of the dual geometry. The current work is a proof-of-concept and we did not yet investigate much which operators drive the renormalization group flow away from the UV or what predictions the bulk reconstruction yields, other than the computation of the thermal entropy per usual identification with the Bekenstein-Hawking entropy of the dual black hole. Given the formidable route to establishing the dual quantum geometry at finite numbers of colors, we advocate that there seems to be a (naive) way to quantify the trustworthiness of classical dual geometry. The holographic entanglement c-function has a sudden jump across the energy scale marked by the deconfinement-confinement transition. This can be associated with the degrees of freedom scaling as O(N 2 c ) and O(N 0 c ) in the respective phases. This jump is likely an artifact of the N c → ∞ limit, and is less drastic for finite values of N c , as first illustrated in [8] for N c = 3 and later in [9] for N c = 2, 3, 4. Nevertheless, if one sticks to the deconfining phase, the reconstructed dual metrics behave smoothly as a function of N c , giving one in principle means to interpolate the emerging classical geometries. There are of course many complications along the way to making this precise, especially if the entanglement entropy is being used as a tool to reconstruct the metric, since the RT formula is expected to be modified in the presence of quantum corrections [47]. It is clearly interesting to study further the entanglement entropies for higher rank Yang-Mills theories. We hope to report on entanglement c-functions in (3+1) dimensions for N c > 4 in the future.
Other possible avenues for further research include the shape dependence on the entangling region. In principle, parts of this information is already encoded in our link updates and one could not only extract, say, the corner contribution [48], but globally extend the en-tangling region to spherical shapes. In addition to this, aspects of multipartite entanglement measures are clearly interesting targets. Simplest extension of the present work being to that of two strips (slabs) and questions related with mutual information, entanglement negativity [49], or purification [50]. In particular, it would be interesting to analyze if in the last case scenario one finds non-monotonous behavior as a function of the energy scale as suggested by the holographic dual quantity [I] known as the entanglement wedge cross section [51].
Entanglement measures are by far not the only data with which one can reconstruct the bulk geometry; work in this direction include [30,[52][53][54][55][56]. Other quantities that have been investigated in the literature include operators [57,58], four-point correlators [59], differential entropy [60,61], fidelity [62], complexity [63], light-cone cuts [64,65], chiral condensate [66,67], conductivity [68], hadron spectra [69], and magnetization [70] with varying degree of assumptions about the dual classical bulk gravity. To all of these cases one could incorporate the Hamiltonian Monte Carlo approach discussed in Sec. 5.2 to also infer the maximum likelihood estimates of the reconstructed duals. In addition to these examples, one can also use the Wilson/Polyakov loops to infer the bulk metric [30,71], though the fundamental string breaking at finite temperature limits the most straightforward application and only part of the geometry is visible to the boundary. In an upcoming work [72] the authors overcome this limitation and reconstruct the full geometry. Besides the murky status of Debye screening at strong coupling as revealed in Appendix B, the ability to use Polyakov loop data also for the reconstruction will further motivate putting extra effort in a more systematic study of quark-anti-quark potentials already in the pure Yang-Mills theory in any dimension. Having simultaneous access to bulk metric both in the string and the Einstein frame by reconstruction from qq potential and entanglement entropy data, respectively, would enable us to make a direct comparison with the beta function of SU(N c ) Yang-Mills theory at strong coupling.
Finally, to address the deconfinement phase transition and lower temperatures in general in the bulk gravity side, one seeks for an alternative ansatz for the metric other than those (3.3), (3.14) where a black hole is assumed to reside at the bottom of the geometry.
Acknowledgments We would like to thank Jacopo Ghiglieri, Marco Panero, and Javier G. Subils for useful discussions. N. J. is supported in part by the Academy of Finland grant. no 1322307, and A. S. and K. R. by the Academy of Finland grants 319066, 320123 and 345070. T. R. is supported in part by the Swiss National Science Foundation (SNSF) through the grants no. 200021_175761 and TMPFP2_210064. The authors wish to acknowledge CSC -IT Center for Science, Finland, for generous computational resources.

Appendices A Lattice Yang-Mills theory
The lattice serves as non-perturbative regulator for the path-integral quantization of field theories. By restricting the matter and/or gauge fields to a lattice, one introduces a cutoff for UV-divergences, and also renders the number of degrees of freedom countable, so that the path-integral measure can be defined in an explicit way. For gauge theories with com-pact gauge group (e.g., SU(N c )) lattice regularization has the additional advantage that it allows for evaluation of the path-integral without gauge fixing, and therefore without the introduction of Faddeev-Popov ghosts. In this appendix we briefly discuss how the lattice regularization for the path-integral of an SU(N c ) gauge theory is performed [73], how vacuum expectation values of observables are determined with Monte Carlo methods, and how continuum physics is extracted from such lattice results.

A.1 Gauge fields on the lattice
Consider a SU(N c ) gauge theory in d dimensions in which some matter field φ couples to the fundamental representation of SU(N c ). In the continuum, the gauge field can then be thought of as a choice of connection, with Lie algebra valued connection coefficients, A µ ∈ su(N c ). These encode how the basis of the local internal spaces, C Nc , (in which the field φ takes values) at infinitesimally neighboring space time points are related. More precisely, the infinitesimal transformation where > 0 is an infinitesimal number and n a d-dimensional unit vector, can be thought of as mapping φ(x + n) from the internal space C Nc x+ n at the infinitesimally displaced point x + n back to the internal space C Nc x at the point x. After this procedure, the image of φ(x + n) and φ(x) are both complex vectors in the same vector space, C Nc x , and the covariant derivative of φ(x) in the direction of n µ is obtained as: which holds for arbitrary n ∈ S d−1 .
The Euclidean Yang-Mills gauge action, is simply the integrated squared norm of the corresponding intrinsic curvature of the connection, also known as field strength tensor, In order to achieve that the covariant derivative of the field φ(x) transforms covariantly under a local gauge transformation Ω(x) ∈ SU(N c ), i.e., where φ = Ω φ, and D µ = ∂ µ + i A µ , one finds by writing both sides of Eq. (A.6) in terms of φ (x), that: and therefore: implying that that the connection coefficients need to transform like under the action of Ω(x). The norm of the curvature of the connection is, however, invariant under gauge transformations, so that the Yang-Mills gauge action is gauge-invariant as required for a gauge theory. Note, that we have adopted the convention of absorbing a factor of g c , i.e., the gauge coupling, into the definition of A µ . The covariant derivative in Eq. (A.1) does therefore not explicitly depend on the gauge coupling g c . Instead, the gauge coupling appears explicitly as a factor g −2 c in front of the gauge action in Eq. (A.4). When transitioning from continuous R d to a d-dimensional Euclidean lattice with lattice spacing a, all space time points become separated by a finite distance a. In order to describe the relation between the internal spaces C Nc at neighboring lattice sites, the infinitesimal transformations from Eq. (A.2) need to be replaced by finite group transformations, meaning that the gauge field would now more conveniently be represented by group-instead of algebra-valued configuration variables. As we anyway can no longer resolve variations of the gauge field below the lattice scale, a, we can for example think of the gauge field being given by the SU(N c )-valued variables where x ∈ Z d is now a d-vector of integers, labeling the lattice sites, and µ refers to the unit vector in µ-direction, so that x + µ is the nearest neighboring site of x in µ-direction.
The link-variable in Eq. (A.10) is a parallel transporter that maps the internal space at site x + µ to the internal space at site x, so that the lattice-equivalent of the covariant derivative, D µ φ(x) of the field φ(x), discussed above, can be written as the finite differencẽ The action for the lattice gauge field {U µ (x)} x∈Z d ,µ∈{1,...,d} can be written as which is known as Wilson gauge action [73], with β g = 2 Nc g 2 0 being the inverse bare lattice gauge coupling, with the dimensionless g 2 0 = g 2 c a 4−d in terms of the continuum theory gauge coupling g 2 c , and U µν (x) are the plaquette variables .
( and Taylor expansion of the gauge field, 16) one can show that the plaquette variable from Eq. (A.13) is given by where the substitution g 2 0 = g 2 c a 4−d has been used, as well as in the sense of a d-dimensional Riemann sum (assuming Riemann-integrability).
As mentioned in the introduction to this section, for gauge theories with a compact gauge group G, (e.g., G = SU(N c )), the lattice regularization has the additional advantage that it allows for an evaluation of the corresponding path-integral without requiring gauge fixing. This is thanks to the existence of a so-called Haar integration measure [74], dU , on compact Lie groups, which is 1. bi-invariant, meaning that: ∀ V ∈ G : dU = d(U V ) = d(V U ) , 2. and normalized, so that: If we now restrict ourselves to the study of a SU(N c ) Yang-Mills gauge theory on a d-dimensional, periodic, Euclidean lattice of finite extent V = N t · N d−1 s , where N s is the number of lattice sites in each of the d − 1 spatial directions and N t the number of lattice sites in temporal direction, the lattice-regularized path-integral can be written as a partition function, which is well-defined and finite for all finite N t , N s . Note that in Euclidean field theory, the extent of the periodic time-direction plays the role of inverse temperature. This remains the case also on the lattice, where the physical temperature T is given by T = 1/(a N t ).

The vacuum expectation value (VEV) of an observable O[U ] can be determined as
The so obtained VEVs are, however, results for the system on a finite lattice, i.e., in the presence of UV and IR cutoffs, and not for the physical continuum theory. In Sec. A.3 we will briefly discuss how continuum physics is extracted from such lattice results.

A.2 Monte Carlo sampling
The path-integral for an SU(N c ) Yang-Mills theory on a finite lattice, given in Eq. (A.20), is well defined and in principle computable. However, due to the large number of integration variables, a direct evaluation with, e.g., quadrature methods is in practice impossible. Consider for example a SU(N c ) gauge field on a d-dimensional, hypercubic lattice with N s = N t = 16. This system is described by d × N t × N d−1 s link-variables, each of which is parameterized in terms of (N 2 c − 1) parameters. The partition function is therefore an integral over (N 2 degrees of freedom, which for the example of SU(N c = 3) in d = 4 dimensions amounts to the number 2 097 152. In order to evaluate such high-dimensional integrals efficiently one has to resort to Monte Carlo methods [75], in which the integration over all possible configurations is approximated by a sum over a finite, representative subset of N S configurations U n , where representative means, that the U n are distributed in accordance with the probability density function Formally we thus get: (A.24)

A.2.1 Markov Process
The representative set U of configurations U n needed for the Monte Carlo method can be generated by a discrete time Markov process: Starting from an arbitrary initial configuration U 0 , new configurations U i are generated subsequently by modifying the previous one. There are many ways to define a Markov process that produces the correct distribution of configurations U n . For simplicity we will discuss here only a very basic implementation using local updates: 1. a single link variable U µ (k) is selected from the current configuration U i .

a new random value
3. the new configuration is accepted with transition probability T (U i+1 |U i ). If it is not accepted, the old configuration is replicated: The index i that counts the number of attempted local updates that have taken place since the beginning of the Markov process, in order to reach the configuration U i , is often referred to as Monte Carlo time. In order to produce configurations with the correct distribution, the transition probability has, in addition to the obvious conditions, also to fulfill a balance condition, which is typically chosen to be the so-called detailed balance condition, 27) which requires that the probability for being in a configuration U and undergoing a transition to a configuration U must be the same as the probability for being in the configuration U and undergoing a transition to the configuration U . The above conditions must be satisfied for all possible configurations U , in order to ensure that the update algorithm is ergodic, meaning that it is such that there is a finite probability for the Markov process to connect any pair of admissible configurations. A set {U s n } n∈{0,...,N S −1} of configurations, that consists of configurations generated by a Markov process, and which are separated by a constant amount s of Monte Carlo time, is called a Markov chain.
A possible choice for the transition probability, constructed to fulfill the detailed balance condition, is given by: A Markov process with this choice of transition probability is known as Metropolis algorithm [76]. A more efficient algorithm is given by the so-called heatbath algorithm [77,78], where the new candidate values, U µ , for the link variables during the local updates are directly generated with the correct probability, so that the Metropolis acceptance test in Eq. (A.28) is unnecessary. As computer time for simulations needs to be finite, also the Markov chains obtained during such simulations will always be finite. This means, that if we start from a random initial configuration U 0 , we will in general have to discard a certain number, N therm , of configurations from the beginning of the Markov chain before we can start using the sampled configurations. This is known as thermalization and is required because the chosen initial configuration is in general very different from a typical equilibrium configuration, so that the initial configuration and its N therm descendants are too unlikely to be present in a finite, representative subset of the full distribution of possible configurations.

A.2.2 Error Estimate
The error in a Monte Carlo estimate for the expectation value of an observable obtained from N S measurements is of order O 1/ √ N S , provided the N S measurements are uncorrelated. For local update algorithms like Metropolis, this is obviously not the case as subsequent configurations differ at most by the change of an individual configuration variable ("at most" as it depends on whether the attempted update is accepted or not). One therefore usually takes measurements only after N sweep (attempted) local updates, where N sweep is set to the number of degrees of freedom in the system and defines a so-called sweep. However, configurations separated by single sweeps are usually still correlated and if these correlations are strong, one should reduce further the frequency by which measurements are taken, in order to avoid wasting computer time with redundant measurements. In any case, possible correlations between subsequent measurements have to be taken into account when computing errors in estimates of expectation values of observables.
Integrated auto-correlation time: Let us denote the true expectation value of an observable O by: with O i being the value of the observable O evaluated on the configuration U i . While this true expectation value (A.29), obtained in the limit of an infinitely long Markov chain, would be an ordinary number, the corresponding result obtained from a finite Markov chain of length N S ,Ō is just an estimator which is itself a random variable which fluctuates around O . 4 As Markov chains generated by Monte Carlo are always finite, we have to know how wellŌ approximates O , i.e., how large its statistical error is: which, by inserting (A.30), reads, To get an estimate for (A.33), one could repeat the simulation M times and approximate is the i-th measurement obtained from the m-th simulation. In this case (A.33) could be approximated by but running the same simulation many times is of course highly unpractical and we would therefore like to estimate (A.33) from the data of just a single simulation.
and where after the last equality sign, we have introduced the integrated autocorrelation time for the observable O: of which (A.37) is the sum, decays exponentially with increasing j and at some point j re becomes smaller than its own statistical error Γ O (j) . From this point on, (A.38) cannot decay further but instead fluctuates around zero with an amplitude ∼ Γ O (j) . If we define the running integrated auto-correlation time, this quantity will behave like a random walker for k > j re , as illustrated in Fig. 17. It is therefore necessary to truncate the summation over j in the definition of the estimator for the integrated auto-correlation time (A.37) and use instead where k max is a constant which is large enough so that (A.40) is a good approximation of true integrated auto-correlation time but small enough so that the auto-correlation function is not yet dominated by statistical noise. The statistical error in (A.40) is then approximately given by [79]: and a convenient choice for k max is to require that in which case the error can also be written as with N ef f O as defined below (A.37). Jack-knife method for error estimation: To accurately determine the integrated autocorrelation time τ int O with the above method, one either has to store a lot of redundant data (redundant because it is highly correlated) which can cause problems due to limited storage, or one has to implement an on-the-fly computation of τ int O while the Markov process is evolving, and this, for every observable that one wants to measure. However, if one is only interested in an accurate estimate of the expected error in the Monte Carlo estimate (A.30) for O , it is much simpler to use a re-sampling method like Jack-knife to do so. If we consider again a set S = {O 1 , . . . , O N S } of N S measurements, this method works as follows: 2. define the i-th Jack-knife set to be given by S i = S\b i , 3. determine the M Jack-knife meansŌ (j) = 1

O ∀j and their average valuē
4. and compute the Jack-knife error as: A priori we do not know the integrated auto-correlation time τ int O in step1, but we can measure the error O as a function of decreasing M , and as soon as the error stops to increase, The advantage of using the Jack-knife definition for the statistical error is, that each of the Jack-knife means, {Ō (i) } i=1,...,M , is computed from all data but one bin, whereas the binwise means, {Õ (i) } i=1,...,M , are computed from the much smaller amount of data contained in single bins. The Jack-knife definition is therefore much more stable if the bin-width b is small. This makes the Jack-knife method also particularly useful to estimate statistical errors in functions of expectation values, as for example the parameter errors in non-linear fits to simulation output, which can easily become instable if the data is noisy.

A.3 Extrapolation to continuum
As mentioned at the beginning, the lattice serves as a regulator for ultra-violet divergences by introducing a large momentum cutoff. But as is the case with any such cutoff: if it is removed without renormalizing the theory, the divergences reappear. In perturbation theory, the renormalization consists of a cutoff dependent redefinition of the fields and couplings in the theory, so that the bare couplings become functions of the cutoff. In lattice field theory, the dependency is usually considered to be in the opposite direction: as the lattice spacing a (which is the cutoff) is absorbed into the definition of the lattice fields in order to render them dimensionless, a is no-longer a free parameter of the theory and can only be changed indirectly as a function of the bare parameters. The theory therefore kind-of chooses for itself the preferred cutoff-scale, depending on the choice of the bare parameters.
This obviously implies that the continuum limit cannot be taken by directly changing a, but rather one has to find and approach a region in parameter-space where the system prefers to reduce the lattice spacing. Such regions can be found by monitoring for example the mass spectrum of the theory as a function of the free parameters (i.e., the coupling constants): for dimensional reasons every lattice mass m is of the form m = a m c , where a is the lattice spacing and m c a dimensionful mass parameter of the continuum theory. If m c is a physical mass, it should become independent of a in the limit (a → 0). Then the lattice mass m has to vanish in this limit and the corresponding correlation length in lattice units, ξ = 1/m, will diverge (cf. Fig. 18).
A point in parameter space where this happens is called a critical point as it corresponds to a second-order phase transition-point of the statistical system described by the Euclidean lattice action under consideration.
In order to determine how the lattice spacing, a, depends on the dimensionless lattice gauge coupling, u = g 2 0 , of the Wilson action (A.12), we require that an arbitrary dimensionless physical observable P (u(a), a), which could for example be a mass ratio, becomes independent of a if a is sufficiently small. This yields the Callan-Symanzik equation: is the lattice version of the renormalization group beta function. For small u we can equate (A.46) to the perturbative expansion of the beta function to get the differential equation, which upon integration yields: where u > 0 and u > 0 are two different gauge coupling values and a(u)/a(u ) the ratio of the corresponding lattice spacing values. In (3+1) dimensions, the lattice gauge coupling does not explicitly depend on a and the first two beta-function coefficients for the SU(N c ) YM theory are universal and with our conventions given by As β 0 > 0, Eq. (A.48) tells us that (a(u)/a(u ) → 0) when (u → 0), and therefore the continuum limit is for SU(N c ) lattice gauge theories in (3+1) dimensions obtained by sending u = g 2 0 to zero, i.e., β g = 2 N c /g 2 0 to infinity. The same is true in (2+1) dimensions, which can be seen even without knowing the beta function: as the continuum theory coupling g 2 c is dimensionful in this case and the lattice coupling is g 2 0 = g 2 c a, the inverse lattice gauge coupling, has to diverge when (a → 0), provided g 2 c remains finite. Simulations have, of course, to be carried out at finite, i.e., non-zero a, but with a sufficiently small so that lattice artifacts are not dominant. The lattice spacing a and the linear spatial system size N s should be chosen so that for all correlation lengths, ξ c,i (in physical units), of relevant excitations one has: so that the lattice spacing is fine enough to resolve all relevant physical processes and the physical size of the lattice is sufficiently large to fit them. The lattice system is then said to be within the scaling window.
Continuum results are obtained by measuring observables on lattices of the same physical size but different values of the lattice spacing a. This means, one picks for example  Figure 18: Illustration of the scaling behavior of a correlation length ξ c in lattice units. If the lattice spacing a is sufficiently small and the physical lattice size N s sufficiently large, so that a ξ c a N s , the system is in the so-called scaling window. The physical value of ξ should then be independent of a and the value of the correlation in lattice units, ξ = ξ c /a, diverges like 1/a when (a → 0). lattices with (N s , N t ) ∈ { (12,12), (24,24), (32,32), (48,48)} and adjusts for each of them the inverse gauge coupling, β g , so that their physical system sizes, (a(β g ) N s , a(β g ) N t ), are the same (see Fig. 18). The resulting data can then be extrapolated to (a → 0), which yields the continuum result for a system contained inside a box of physical size To obtain infinite volume continuum result, this procedure has to be repeated for different physical system sizes V so that the resulting data can finally be extrapolated to infinite volume.

B Quark-anti-quark potential
In this appendix we slightly extend the discussion in the main text to cover also the Polyakov loop which give rise to a potential between quarks and anti-quarks. This is to illustrate that the D2-brane background, which lead to the rather unconventional power law behavior for the entanglement entropy S EE ∝ −4/3 , predicts also an unexpected power law behavior for the quark-anti-quark potential at small separation ∝ L −2/3 . We will see that the lattice data that we will generate seems to conform with this and thus lends more support on the usage of the gravity dual of (2+1)-dimensional super Yang-Mills theory at largish energy scales. In addition to this matching, we learned from EE studies that at finite temperature, the black hole present in the D2-background also explains the behavior of the entanglement entropy at large subregions sizes. It is therefore not completely unreasonable to ask if the large separation behavior of the qq-potential be also within grasp via D2-brane black holes in this regime. Indeed, our frugal results from the lattice seem to conform with ∝ L −10/3 behavior at large separation, extractable from the D2-brane background. If this result persists under further scrutiny, then the holographic approach predicts that the analytic continuation of the lattice results to real time cf. [39] would result in an imaginary part of the potential linear in the qq separation L.

B.1 Static quark potential from thermal D2-brane background
To derive these qq-potential results we start with the metric (3.14) with a(z) = 1 and b(z) = 1 − z z 0

5
. This is a black hole geometry with a horizon at z = z h . The potential is found by studying a string whose endpoints are at the boundary, separated by a distance L in the x 1 -direction. The action of such a string is where g 2 is the metric induced on the string worldsheet. The action can be written in terms of the profile function x 1 = x 1 (z) which determines the shape of the hanging string in the bulk as where τ is the time interval spanned by the string worldsheet. The qq-potential can be found by studying the equation of motion for the string profile x 1 (z). The Euler-Lagrange equation is The solution with the boundary conditions x 1 (z = 0) = L/2 and where F 1 is the Appell hypergeometric function. The constants L and z * are the quark separation and string turning point, respectively. They are connected by the equation The string action can be evaluated by plugging x 1 (z) from (B.4) into (B.2). One finds where is the usual UV-cutoff to regulate the otherwise divergent action. Now we can find the UV-behavior of the potential by studying (B.5) and (B.6) in the limit z * → 0 L = √ πΓ Together, these formulas imply that when z * is small the action behaves as The qq-potential is defined There are different ways to do the regularization, often the potential is computed from the difference of the string action and the action of two straight disconnected strings hanging from the boundary. This time we choose a different method, namely we find it useful to regulate so that the real part of the potential goes to zero at large separations. We will shortly show that the latter term in (B.10) is a constant. Therefore, the potential indeed is V (L) ∝ L −2/3 at small separations. In order to study the potential at large separations we need to analyze the string profile carefully. Usually the large separation limit is somewhat trivial as often the string is assumed to break at some critical separation making the potential constant at larger separations. It is possible, however, to extract non-trivial IR-behavior of the potential by analytically continuing L(z * ) and S N G (z * ) to the complex z * -plane [80]. The quark separation L is always non-negative and real but the turning point z * need not be. In Fig. 19 we show a curve along which L is real and non-negative, corresponding to complex string configurations. The curve represents the solutions to (B.5) when solving for z * and setting L to fixed nonnegative real values. At small separations z * is real as usual. The separation at which the string would break into the disconnected configuration happens in this region. Shortly after that separation the turning point becomes complex. Also the conjugate of this curve solves the equation but we choose to focus on the Im[z * ] ≥ 0 branch.
This can be used together with (B.6) to see that the action at large separations behaves as In the definition of the qq-potential, we can now see that the real part of the IR-action will cancel the 1/ -divergence present in S N G (L).

B.2 Static quark potential from the lattice
The lattice gauge theory formalism was introduced by K. G. Wilson in 1974 to show that quarks are confined in QCD [73]. The static quark potential was therefore from the beginning of primary interest in lattice Monte Carlo studies of SU(N c ) gauge theories [81] and various observables have been employed to study its properties. In particular Wilson loops, which are traces over parallel transporters along closed loops, and Polyakov loops, which are Wilson loops that wrap around a compact direction (usually the Euclidean time direction), played an important role in understanding the non-perturbative properties of QCD, as, e.g., showing that SU(2) pure gauge theory is both, confining and asymptotically free [82], and that the asymptotically free and the confining region of parameter space are separated from each other by a deconfinement transition [83].
With increasing computer performance and the development of improved simulation techniques [77] these non-perturbative studies of the properties of the static quark potential have been extended first to SU(3) [84][85][86] and then to higher N c and spacetime dimensions different from (3+1) [87]. Results at T = 0 in (2+1) dimensions can for N c = 2, for example, be found in [88,89], for N c = 3 in [90], and for N c = 5 in [91]. Recall that the static quark potential can also be computed perturbatively [37].
At T = 0 the determination of the static quark potential using Polyakov loop suffers from a bad signal to noise ratio as the magnitude of Polyakov loop correlation function drops exponentially with 1/T . At T = 0 the static quark potential is therefore usually determined from Wilson loops [85]. As we are in the present work interested not in the zero but in high temperature properties of SU(N c ) gauge theory, we will use Polyakov loop correlators to determine the static quark potential.

B.2.1 Static quarks and Polyakov loops
In this section we illustrate the relation between Polyakov loops and static quarks. Starting point is Wilson's lattice fermion action with ψ(x) andψ(x) being Euclidean Dirac spinors appropriate for d dimensions and D x,y [U ] is the corresponding Wilson-Dirac operator, Now, if we are interested in static quarks, i.e., in quarks that do not perform spatial hops but simply move along straight lines around the time direction, we can in Eq. (B.14) drop the spatial hopping terms, S x y , and keep only the temporal hopping terms, T x y . In this case, the fermion determinant in Eq. (B.16) simplifies significantly. If we stick for the moment to d = 4, it takes the form: where the subscript c of the det c -operator is meant to indicate that it acts only in color-space (i.e., that the operand is a matrix with only color indices) and the product runs over all spatial locationsx.
As can be seen from Eqn. (B.18) and (B.17), the static quarks couple to the gauge field through the Polyakov and anti-Polyakov loops. The Polyakov loop, P (x) at spatial location x, is given by the ordered product of all temporal link-variables (cf. Eq. (A.10)) overx, to form a parallel transporter that connects the site x = (x, t = 0) with itself by wrapping once on a straight temporal trajectory around the periodic time-direction: Its hermitian conjugate, P † (x), does the same thing but in opposite direction, as illustrated in Fig. 20.
As can be seen from Eqn. (B.20) by counting the powers of e µ resp. e −µ : having a static quark or anti-quark winding around the time direction over a spatial sitex is accompanied by a Polyakov loop resp. anti-Polyakov loop located on that site.

B.2.2 Static quark potential from Polyakov loops
We can now define the static quark potential, V (L, T, β g ), in terms of the spatial correlation function between a traced Polyakov loop, P , and a traced anti-Polyakov loop, P † (cf. Fig. 20), using the relation [92], with L = |ȳ −x| and T = 1/N t (in lattice units). Note that Eq. (B.21) is particularly useful at finite temperatures, i.e., when N t is not too large. As the formula indicates, at low temperatures, the magnitude of the correlation function decays quickly with N t = 1/T and develops a bad signal to noise ratio. If one is interested in the static quark potential at low or zero temperatures, one therefore typically uses a different approach which is based on so-called Wilson loops [93]. The potential V (L, T ) in Eq. (B.21) has an explicit temperature-dependency, which is particularly prominent for T /T c > 1, where the potential approaches for L T −1 a temperature-dependent plateau value. We will subtract this plateau value and consider the subtracted potential, the binding energy (in units of the critical temperature scale),

B.2.3 Setting the scale
We use the critical temperature, T c , as reference energy scale. To be able to do so, we need to determine the ratio T /T c as function of the inverse lattice coupling β g and the spatial and temporal lattice sizes N s , N t . We therefore want a function C R (β g , N t , N s ) = T /T c , which can be obtained along the lines described in [94,. Using the expression given in [94] for the pseudo-critical lattice coupling as function of N t and N s , as well as the function where [94] The subscript c in Eq. (B.27) indicates that T /T c ≈ C R,c (β g , N t , N s ) is only valid if β g is sufficiently close to its pseudo-critical value β g,c (N t , N s ).
To obtain an expression for C R (β g , N t , N s ) that is valid for arbitrary values of β g , N s and N t , We note that (a T c ) depends only on β g and N s , and the lattice spacing a itself is controlled merely by β g . At constant β g and N s (which implies also constant a(β g )), we therefore have [94]: which allows us to define the function C R (β g , N t , N s ) = T /T c for arbitrary values of β g , N s and N t by: i.e., N t (β g , N s ) is the value of N t ∈ N for which (β g − β g,c (N t , N s )) assumes its smallest, positive value for given β g and N s .

B.2.4 Data and interpolation to arbitrary temperatures
In Tab. 2 we summarize the simulation parameters for which we have produced data. As can be seen, although the data for different lattice spacings covers a similar temperature range, the simulated temperatures are in general different. In order to be able to directly compare the static quark potentials obtained at different lattice spacings, we will have to interpolate between the simulated temperature values.  We perform this interpolation as follows: we pick one of the simulated β g and a value of L, and fit for this choice the form   with the same values of α (L,βg) and γ (L,βg) . However, as our spatial lattices are relatively large, δ (L,βg) is typically very small (cf. Fig. 23, right) and whether one uses Eq. (B.33) or sticks with Eq. (B.32), does not give rise to significantly (within error bars) different results.

B.2.5 Scaling of deconfined potential at short and large distances
As discussed in Sec. B.1 holographic considerations suggest that the static quark potential in the deconfined phase of a strongly coupled gauge theory in (2 + 1) dimensions should follow different power laws at short and long distances, namely V (L T 1) (T, L) ∝ L −2/3 and V (L T 1) (T, L) ∝ L −10/3 . In Fig. 24 we fit these power law ansatzes to the short and long distance pieces of the interpolated potential (B.32) (as function of L T c ) at fixed T /T c = 2.646 (left) and fixed T /T c = 1.8 (right) for different lattice spacing values. In order to stay away from the region around L T ≈ 1, where the log-log plots in Fig. 24 have clearly visible curvature, we  (red). At short distances, well below L ∼ T −1 , the behavior of the potential appears compatible with an asymptotic power-law, ∼ L −2/3 (dashed lines). The displayed fits were performed on data for L < 0.3 T −1 (where available). At distances above L ∼ T −1 , the potential appears to behave like ∼ L −10/3 (dotted lines). Here the fits were obtained on data for L > 1.7 T −1 .
included for the short distance fits only data with L T < 0.3, and for the long distance fits only data with L T > 1.7. The short distance fit could therefore only be performed for the two smaller lattice spacing values, as for the largest one with a T c = 0.98, there is only one point satisfying the criterion L T < 0.3 .
At the current stage, these fits should of course not be taken too seriously, as the available data is of limited statistics and does also not sufficiently cover the asymptotic shortand long-distance regimes. In both regimes, we see at most the onset of a possible linear behavior of the subtracted potential in the log-log plots in Fig. 24. It should also be mentioned that at short distances, the displayed data contains also the points corresponding to L = 1 a, which can be expected to suffer from relevant finite-lattice spacing effects. At large distances one can infer that the data does not cover a sufficiently large range to decide whether the potential really approaches a power law and starts to follow a straight line in the log-log plots in Fig. 24, or whether its slope will continue to become even steeper with increasing L.
To this end, we attempted to fit also a screened power law, i.e., the form avoid the most-severe finite lattice spacing effects, and the points with L T c > 1. 35 where for the two smaller lattice spacings the finite lattice volume would start to become visible. As can be seen, these fits worked out quite well for the given quality of the data. In order to confirm that the use of Eq. (B.33) instead of Eq. (B.32) to interpolate our data for the subtracted static quark potential between the simulated temperature values, would not significantly affect the fit results, we show in Fig. 26 the analogous plots to those in Fig. 25, using the improved interpolating potential form Eq. (B.33). These latter fits might appear slightly cleaner than the ones in Fig. 25, in regards of the lattice spacing dependency of the fitted parameters. But, the fitted parameters in Fig. 25 and Fig. 26 are mutually compatible within error bars.

C Change of temporal boundary conditions and gauge invariance
Recall from Sec. 2.2 that when updating the boundary of region A by changing the temporal boundary conditions over spatial links, so that they change from region B to region A or vice versa, we undergo transitions between different partition functions {Z i } i=0,... , as introduced above Eq. (2.20). We are interested in the free energy difference between subsequent Z i , which we determine numerically by probing the overlap between the distributions of gauge field configurations of neighboring Z i . When changing from one Z i to a neighboring one, it turns out that it makes a huge difference whether the spatial link over which the temporal boundary conditions are changed has before and after the update at least one of its ends connect to a spatial site that is either completely in region A or completely in region B (cf. Fig. 27, left), or the spatial link touches before and/or after the update at both ends a spatial site at which spatial links from both regions meet (cf. Fig. 27, right).
In order to illustrate why this is the case, we look at a transition between some Z i and corresponding neighboring Z i+1 , which will require a spatial link to change from region B to region A. We recall from Fig. 5 that when a spatial link changes from region B to region A or vice versa, two temporal plaquettes over that spatial link, denoted in Fig. 5 by P 1 and P 2 , switch their top links, due to the change of temporal boundary conditions. Such an update is in general gauge-dependent: if it were, e.g., possible to use a local gauge transformation to make the to-be-swapped links the same before doing the change of temporal boundary conditions, then the Euclidean action would remain unchanged when changing form Z i to Z i+1 . Such a gauge-dependency is okay here as we are sampling the overlap between gauge  Figure 27: The two panels show the spatial links corresponding to regions A (red) and B (blue) at two intermediate stages during the process of interpolating between the situations depicted on the left-resp. right-hand side of Fig. 5. The spatial link that is highlighted by a red frame in the lefthand panel is at the moment blue and belongs therefore to region B. The spatial site x 1 is touched only by spatial links belonging to region B, while the spatial site x 2 is touched only by spatial links belonging to region A, except for the highlighted spatial link. If the highlighted link changes from region B to region A, the spatial site x 2 will be surrounded merrily by spatial links from region A, while the spatial site x 1 will now also be touched by one link from region A. However, before and after the change of the highlighted spatial link from region B to region A, one of its ends is connected to spatial site that is touched by spatial links of only one region. In contrast: the spatial sites at both ends of the highlighted spatial link in the right-hand panel are always touched by spatial sites from both regions, A and B, no matter whether the highlighted link is still in region B or changes to region A.
configuration distributions of different systems with different temporal boundary conditions and therefore different sites identified with each other. The latter, i.e., the fact that in the two systems different sites are identified, affects locally the amount of gauge-freedom the system has as one cannot perform independent local gauge-transformations on identified sites. The overlap between the gauge configuration distributions of neighboring Z i has therefore an entropic component coming from the change in the amount of gauge freedom due to the identification of different numbers of sites.
Consider as a concrete example the situation depicted on the left-hand side of Fig. 27: if the indicated link between x 1 and x 2 belongs to region B (blue), we are in the system described by Z 1 , and if the link changes to region A (red), we would be in the system described by Z 2 . If the link under consideration is currently part of region B, i.e., we are in the Z 1 -system, the spatial site x 1 is touched only by spatial links that belong to region B, and as explained in the left-hand panel of Fig. 28, we could therefore apply independent gauge transformations over x 1 to make the two links the same that would be swapped between the plaquettes P 1 and P 2 from Fig. 5 when the temporal boundary conditions are changed to those from region A. The same is true for the opposite move, where the marked spatial link from the left-hand panel of Fig. 27 between the spatial sites x 1 and x 2 is initially part of region A. In this case the site x 2 is touched only by links from region A and would allow for the required local gauge transformations to make the to-be-swapped top-links of plaquettes P 1 and P 2 identical. We can therefore define a Metropolis update that connects a gauge configuration of Z 1 to a gauge configuration of Z 2 : if the system currently belongs to Z 1 , the update proposes to pick a local gauge transformation over x 1 that makes the link-swap trivial, to then perform the link swap, and to finally do an inverse gauge transformation over x 2 ; if the system belongs currently to Z 2 , the Metropolis update would do the same thing in the opposite direction, i.e., with initial local gauge transformation over x 2 and inverse gauge transformation over x 1 after the swap. As the freedom in choosing the appropriate gauge transformations is for both directions the same, and the action does not change during the move, the corresponding transition probability is always 1. Although the way in which sites are identified changes during such an update due to the change of temporal boundary conditions, the number of identified points remains the same. Note also that, although we can achieve that the swapped links themselves are the same before and after such a move, the gauge filed around them gets changed, as the local gauge transformations that are applied at the beginning and end of the move are applied at different spatial locations. Figure 28: Left: the four sites highlighted by a red, orange, blue and green circle are fully inside of one of the two regions A (red,orange) or region B (blue,green) and therefore touched only by links that are all subject to the same type of temporal boundary conditions: either those of region A, or those of region B. If is therefore possible to apply independent local gauge transformations at the green and blue, as well as the red and orange sites, which will affect the correspondingly colored links. Right: if one applies a local gauge transformation to the site highlighted by a solid-lined black circle, which is touched by links from both regions, A and B, the situation is more complicated. In region B, a gauge transformation at this site changes the link l 1 . From the point of view from region A, however, the link l 1 is connected to the site highlighted by the dashed black circle, on which we therefore have to apply the same local gauge transformation as on the site with the solid black circle in order to leave the gauge action unchanged, which will affect all the links highlighted by dashed lines. The same conclusion is reached by focusing on the link l 2 instead of l 1 . The gauge transformation on the site with the black circle changes in region A directly the link l 2 . In region B, the link l 2 is connected to the site with the dashed black circle, on which we therefore have to apply the same gauge transformation as on the site with the solid-lined black circle. As a second concrete example consider now the situation depicted in the right-hand panel of Fig. 27. The spatial link emphasized there by the red frame connects the sites x 1 and x 3 . If the link belongs to region B, the system is described by Z 7 and if the link belongs to region A by Z 8 . Both spatial sites, x 1 and x 3 , are touched by links from both regions, A and B. As explained in the right-hand panel of Fig. 28, this means that it is not possible to find a local gauge transformation that makes the links, that are swapped between P 1 and P 2 in Fig. 5 when the temporal boundary conditions change, the same. Spatial sites that are touched by spatial links from both regions, A and B, are (indirectly) subject to both temporal boundary conditions simultaneously, meaning that the corresponding site at time t = 0 over such a spatial site is not just identified with either the site at the same spatial location at time t = N t (as inside region B) or time t = 2 N t (as inside region A) but with both.
Note that in (1 + 1) dimensions, the situation from the second example cannot occur and one is always in the analogous situation to the first example. In (1 + 1) dimensions, the boundary update of region A is always trivial and the derivative of the entanglement entropy with respect to the width, , of region A has therefore to be identically zero.