Lattice study of R\'enyi entanglement entropy in $SU(N_c)$ lattice Yang-Mills theory with $N_c = 2, 3, 4$

We consider the second R\'enyi entropy $S^{(2)}$ in pure lattice gauge theory with $SU(2)$, $SU(3)$ and $SU(4)$ gauge groups, which serves as a first approximation for the entanglement entropy and the entropic $C$-function. We compare the results for different gauge groups using scale setting via the string tension. We confirm that at small distances $l$ the entropic $C$-function $C(l)$ for the slab-shaped entangled region of width $l$ scales as $N_c^2 - 1$ in accordance with its interpretation in terms of free gluons. At larger distances $l$ $C(l)$ is found to approach zero for $N_c = 3, 4$, somewhat more rapidly for $N_c = 4$ than for $N_c = 3$. This finding supports the conjectured discontinuity of the entropic $C$-function in the large-$N$ limit, which was found in the context of AdS/CFT correspondence and which can be interpreted as transition between colorful quarks and gluons at small distances and colorless confined states at long distances. On the other hand, for $SU(2)$ gauge group the long-distance behavior of the entropic $C$-function is inconclusive so far. There exists a small region of lattice spacings yielding results consistent with $N_c=3,4$, while results from other lattice spacings deviate without clear systematics. We discuss several possible causes for discrepancies between our results and the behavior of entanglement entropy in holographic models.


Introduction
Entanglement entropy has become a heavily studied field of research in recent years. For example, it is widely used in quantum information theory, where many other entanglement measurements, as e.g. mutual entropy, exist [1], as well as in connection with the gauge/gravity duality [2]. It can also be used as a universal order parameter of quantum phase transitions, as it is done, e.g. for 2 + 1 dimensional topological field theories, where one relates entanglement entropy to the quantum dimension [3].
More generally, in the context of quantum field theory, entanglement entropy can be considered as a counter of the effective number of degrees of freedom, which is related to the central charge in two-dimensional conformal field theories (CFTs) [4]. Much like the central charge, entanglement entropy changes monotonically under renormalization group (RG) flow, and obeys a generalization of Zamolodchikov's c-theorem [5]. In higher-dimensional field theories, entanglement entropy can be related to Cardy's Afunction which generalizes Zamolodchikov's c-function [6] and which is also monotonous under RG flow [7,8].
In strongly coupled quantum field theories, the entanglement entropy turns out to be notoriously hard to compute from first principles either numerically or analytically unless the models are sufficiently symmetric or low-dimensional, see e.g. [9] for an overview. In the light of this difficulty, the AdS/CFT duality offers an attractive alternative route to calculate it in terms of the area of minimal surfaces in a dual gravitational theory [2]. While this computation is a priori justified only for large number of colors N c and strong coupling in N = 4 super Yang-Mills theory, it may still provide valuable guidance to understand the generic behavior of entanglement entropy in strongly coupled gauge theories. Indeed, the computation provided in [10] using a confining background predicted an abrupt transition at some l = l c in the dependence of entanglement entropy on the characteristic size l of the entangled region, at which the scaling of the entanglement entropy changes from O(N 2 c ) to O(1). Since the entanglement entropy counts the effective number of degrees of freedom at energy scale 1/l, the transition from O(N 2 c ) to O(1) scaling of the entropy can be interpreted as a transition between colorful quarks and gluons at small distances and colorless hadrons and glueballs at long distances. While such scaling laws at asymptotically short and large distances can be readily understood in terms of asymptotic freedom and confinement, the discontinuous nature of the transition is a non-trivial prediction for theories with a Hagedorn-type spectrum, which also include holographic models.
This non-trivial prediction was taken as a motivation to analyze the spatial dependence of the entanglement entropy in pure SU (2) [11] and SU (3) [12] Yang-Mills theory, approximating it with a Rényi entanglement entropy for the sake of numerical convenience. Somewhat surprisingly, there is a qualitative difference in the results: while [11] seems to show an abrupt transition after an initial enhancement, [12] is compatible with both an abrupt and a smooth transition within statistical errors, but leans more towards a smooth transition.
In this work, we continue this line of research and present a comparative highstatistics study of the Rényi entanglement entropy in lattice gauge theories with SU (2), SU (3) and SU (4) gauge groups. Our aims are to study the scaling of entanglement entropy with N c and to check whether the dependence of the entanglement entropy on the size of the entangling region approaches a discontinuous function at large N c . Next to better understanding SU (N c ) Yang-Mills theory, this is of obvious interest to quantify differences between holographic and real world gauge theories.
We work with pure SU (N c ) Yang-Mills theory, as already computing the second Rényi entropy to approximate the entanglement entropy turns out to be computationally expensive. Owing to this, we are also not able to exceed a lattice size of 16 3 × 32 or perform a continuum extrapolation 1 .
We are able to clearly establish the scaling of entanglement entropy with the number of gluon states N 2 c − 1 at short distances. While the transition between O(N 2 c ) and O(1) scalings of the entropy at short and large distances appears to be smooth for our data for all N c which we consider (as it should be for a finite-N c gauge theory on a finite lattice), we find indications that for N c = 4 the entanglement entropy changes faster in the transition region than for N c = 3. This supports the scenario of a discontinuous transition in the large-N c limit, in agreement with the arguments of [10] which predict a discontinuous behavior of entanglement entropy for any field theory with a Hagedorn-type spectrum. This paper is structured as follows: we start by reminding the reader of basic definitions and properties of entanglement entropy, also from a holographic perspective, in section 2. In section 3, we will recall how one can measure entanglement entropy on the lattice, approximating it with the second Rényi entanglement entropy. Section 4 contains our numerical results, their interpretation, and a comparison between different N c . Finally we will conclude and point out possible future directions in section 5.

Entanglement and Rényi entropy
Intuitively, entanglement entropy measures the quantum correlations of two subsystems. Hence, we start with a bipartition Σ t = A ∪Ā of a fixed time slice Σ t as shown in the left picture of figure 1. Next we write the Hilbert space H as a product space of the bipartition, i.e. H = H A ⊗ HĀ. For non-gauge theories this factorization is unique, but for gauge theories there are ambiguities how to factorize the Hilbert space. In this work we use the so-called maximally gauge invariant way [13,14]. In this method, one cuts along the gauge links and then decides to which region, A orĀ, the link degree of freedom belongs to. This is visualized in the right picture of figure 1. To measure the entanglement between the two regions, we take a partial trace over region A and define the reduced density matrix to bê ρ A = trĀ (ρ) (2.1) whereρ = |ψ ψ| is the density matrix of the whole system. The entanglement entropy is defined as the von Neumann entropy of the reduced system, i.e.
This definition can be generalized to Rényi entropies which are given by [15,16] For most systems, as also for the system we will consider, this definition can be analytically continued to q ∈ R + . With this, one can show that the limit q → 1 of the Rényi entropies gives the entanglement entropy, i.e.
We emphasize that computing S EE is hard in Monte-Carlo simulations due to the logarithm of the reduced density matrix, whereas S (q) for q ∈ N, q ≥ 2 is relatively simple to evaluate. On the lattice, it seems most economical to use S (2) as an approximation for S EE , which can be motivated by looking at the q-dependence of S (q) in simplified models [17]. It should be noted that the alternative formula leads to the same result once the derivative is approximated as f (q = 1) = f (q = 2) − f (q = 1) on the lattice, i.e. one also measures precisely S (2) . We note that one may use larger values of q as an additional input to extrapolate to q = 1 instead of working with q = 2, but this is beyond the scope of this work as we do not expect this to change the qualitative behavior when the entangling region is varied.
In the UV limit the entanglement entropy diverges as [18,19] In what follows we will restrict ourselves to a bipartition of the three-dimensional time slice Σ t into a slab of width l and its complement. This geometry of entangling region A was also considered in most previous works [10][11][12]. With such a geometry, the leading contribution to the non-divergent part of the entanglement entropy is negative and scales as l −2 at small l for d = 4 dimensions [10]. To have an observable which is finite at all l, we define the so-called entropic C-function [20] While this function roughly corresponds to the A-function in the terminology of [6][7][8], here we follow the terminology of [11,12,20] and refer to C(l) as an entropic C-function, by analogy with central charge in two-dimensional CFTs. As already discussed in the introductory section, C(l) counts the number of effective degrees of freedom at energy scale 1/l. By virtue of the higher-dimensional generalization of Zamolodchikov's ctheorem [7,8], C(l) should monotonically decrease with l.
The entropic C-function can also be calculated holographically using the Ryu-Takayanagi formula [2] which states that entanglement entropy is proportional to the area of the minimal surface whose boundary coincides with the entangling region ∂A, i.e.
This formula was applied to confining backgrounds in [10] to find that there exists always a connected and a disconnected solution to the equations of motion which one has to consider to find the minimal area. The connected solution depends on the length l of the subsystem while the disconnected one is constant for all l. Since these two solutions do not intersect smoothly we can conclude that the entropic C-function has a jump at l = l c . At l < l c C(l) scales as N 2 c (to the leading order of 1/N c expansion), while at l > l c C(l) is of order of unity, as could be expected due to confinement at large distances.
In the next section we discuss how to calculate the entropic C-function on the lattice.

Measuring Rényi entropy on the lattice
In this section we will describe how we can measure Rényi entropy on the lattice. This has already been done for SU (2) and SU (3) pure Yang-Mills theory in [11] and [12] respectively. Therefore we will only briefly review the method and refer the interested reader to these publications.
From equation (2.3) we see that we have to calculate the trace of powers of the reduced density matrix, i.e. tr A (ρ q A ). A strategy to solve this problem using path integrals was first proposed in [17], which we will review in a graphical way adapted to the lattice. We start with the calculation of matrix elements of the density matrixρ which we can write as In this work, we always refer to the density matrix of the groundstate, i.e.ρ = |0 0|, which is obtained in the limit β → ∞. We have to understand the picture in this formula such that we integrate over all field configurations from Euclidean time τ E = −∞ to τ E = ∞ with |ψ1 /2 as boundary conditions. On the lattice, this integral is discretized and then approximated using Monte Carlo techniques. With this in mind, we can calculate the reduced density matrix by tracing over regionĀ, i.e. (3.2) In the last step, we deformed the lattice such that we automatically get rid of the sum by identifying corresponding states |ψ . It is important to note that the same number of degrees of freedom label a regular lattice and a deformed lattice. In particular, the diagonal link in the deformed lattice carries the same group element as the middle lower link. This feature will later allow to compute the derivative in the holographic C-function by comparing different cut lengths l obtained from the same underlying lattice configurations. Powers of the reduced density matrix can be computed by glueing copies of the above lattices together, i.e.
T TĀ A Figure 2: Visualization of the deformed geometry as it was already shown in [12]. The periodicity in region A (right, red) is 2T while inĀ (left, green) we have two copies with periodicity T . The crosses, squares and triangles indicate the identifications which have to be made.
This procedure easily generalizes to arbitrary integer powers as where Z[T ] is the partition function of a normal lattice and the partition function Z[A, q, T ] refers to a lattice which has a temporal period of q ×T in A and q subsystems with period T inĀ. Figure 2 provides an alternative graphical depiction already used in [12]. Using (2.3) and (3.4) we can write the Rényi entropies as where we used the relation F = − log Z between the free energy F and the partition function Z. As already mentioned, we will approximate the entanglement entropy by the second Rényi entropy S (2) . To get an approximation for the entropic C-function, we have to calculate the derivative of S (2) with respect to l. We do this by using the central finite difference approximation, i.e.
The term ∼ F [T ] in (3.5) vanishes since it does not depend on l.
So far, we have reduced the problem of calculating the entropic C-function to measuring differences of free energies on the deformed lattice geometry defined above.
Such differences can be calculated using the relation between the free energy and the partition function and the Fundamental Theorem of Calculus [21,22]. Namely, the difference between the free energies F 2 = − log Z 2 and F 1 = − log Z 1 which correspond to the lattice actions S 2 [U ] and S 1 [U ] can be represented as an integral of the form where Z(α) is the partition function which interpolates in some way between Z 2 for α = 1 and Z 1 for α = 0. The simplest partition function satisfying this requirement can be constructed by linearly interpolating between the actions S 2 [U ] and S 1 [U ]: refers to the path integral over all gauge links U . Then we can calculate the derivative with respect to α to obtain (3.8) In this formula · α denotes the average with weight of the interpolating partition function defined above. Coming back to our earlier remark, we note that it was crucial in the last step that we can interpret the same underlying lattice data as different cut lengths.
With this procedure, we can calculate the entropic C-function by the following steps:

Generate gauge configurations with interpolating action
(3.9) on the deformed lattice (see figure 2). The subscripts l and l + 1 refer to the length of the cut.
2. Measure S l+1 − S l on these configurations for several values of α 3. Integrate over α by interpolating with cubic splines.

Calculate the entropic C-function by
(3.10)

Numerical results
We have shown in the previous section how to calculate Rényi entropies on the lattice. In this section we will first show the details of our lattice calculation and afterwards we are going to present results for the entropic C-function in SU (2), SU (3) and SU (4) pure Yang-Mills theory. To implement the method described in the last section, we use a standard Wilson gauge action given by We are using a pseudo-heatbath algorithm to update the gauge configurations. Doing this we have to calculate the closing plaquettes for the link which we want to update. It is crucial in this step that we carefully choose the nearest neighbors of the lattice sites due to the non-trivial lattice topology.
To generate random SU (N c ) matrices we use the Cabbibo-Marinari algorithm [23]. To avoid autocorrelation we wait 100 sweeps between successive measurements. For all sets of parameters we use a lattice of size N 3 s × (q × N t ) = 16 3 × (2 × 16). Since the resolution of the entropic C-function in l would be very poor if we only considered one single β value, we have to set the scale for each lattice and compare results for several different lattice spacings. To do so we use the string tension √ σ as determined in [24] for a range of β values and then measure everything in units of √ σ. The parameters of the gauge field configurations which we have used in our analysis are summarized in table 1. The numbers in these tables are the numbers of configurations we gathered for each α value. As mentioned in the last section (cf. equation (3.10)), we have to integrate the interpolating action over α. Therefore, we have to study carefully how many interpolating points we have to consider. We did this by increasing the number of points in α and then checking the integral value for convergence. As in [12] we find that we need 11 interpolating points for SU (3), but we need 21 for SU (4). For SU (2) we also use 11 interpolating points instead of the 6 which were used in [11]. The increased number of α values can be explained by the larger lattice volume and the, therefore, larger integrand. Besides the increased computational cost due to larger color matrices, this gives an additional factor of ∼ 2 for N c > 2 in computation time.
The entropic C-functions which we obtain from our simulations are shown in figure 3. In these plots, the different colors and markers correspond to different lattice spacings.
For sufficiently small values of l, the C-function appears to be approximately constant for all values of N c . Since the entropic C-function (2.6) is defined in such a way that it yields a constant for non-interacting conformally invariant fields, this behavior could be expected in the short-distance asymptotic freedom regime where entropy is saturated by weakly interacting massless gluons. To extract the corresponding constant values of the C-function, we fitted a constant function to the first data points where l √ σ is less than 1.1, 1.3 and 0.7 for SU (2), SU (3) and SU (4) respectively. We also varied the fitting range slightly and found a mild dependence. This fit is shown by the dashed black line with the blue error band. The constant values and χ 2 /d.o.f. for these fits are given above and below this line, respectively.
At larger cut lengths l and for N c = 3, 4, we observe the expected monotonic (within statistical errors) decrease of C(l) with l down to a value which is compatible with zero within our statistical errors. This fall-off does not look like the discontinuous transition predicted by a holographic calculation [10], but this is expectable for a finite lattice and finite N c .
On the other hand, for SU (2) the data looks qualitatively different and does not exhibit any clear trend for l √ σ 1. It is striking however that for a √ σ = 0.220, the data points show the expected transition around l √ σ = 1.0. One possible explanation of the absence of a clear trend in the data for other lattice spacings is that a √ σ ≈ 0.220 gives only a small window of lattice spacings where both finite volume and discretization errors are small enough. For finite volume effects, this is supported by the temperature on our lattice being already half the critical temperature for a √ σ ≈ 0.18 [25], see also [26]. We emphasize that this discussion is only qualitative so far and needs to be revisited using larger lattices, which is however outside of the scope of the current paper. In the concluding section 5 we also discuss another possible explanation related to the different orders of finite-temperature phase transitions in SU (2) and SU (N c > 2) gauge theories.  SU (4) Figure 4: The entropic C-function for SU (2), SU (3) and SU (4) rescaled by N 2 c − 1, the number of gluons. The fits to these data are according to a regularized step function (4.2). The first three red points (SU (3)) were excluded from the fit due to large discretization errors, but show a trend towards the fit in the continuum limit.
To compare the data points for different gauge groups, we jointly plot in figure 4 the data for the entropic C-function for SU (2), SU (3) and SU (4). For a better comparison across different gauge groups, we rescale the data by a factor N 2 c − 1, which is the expected scaling in the small-cut regime. In order to keep the plot readable, we no longer distinguish data points with different values of lattice spacing. For small cut lengths, this plot again confirms the scaling of entanglement entropy with N 2 c − 1.
In order to check whether the transition between small-cut and large-cut regimes becomes steeper for larger N c , as predicted by the holographic calculation [10], we fit the data to a regularized step function given by The resulting fit with error bands is visualized in figure 4  While the statistics is not large enough to reach a definite conclusion, we observe a trend towards a steeper slope at N c = 4, which is consistent with the scenario of discontinuous transition in the large N c limit.  Table 1: SU (2), SU (3) and SU (4) configurations for the used β values and cut lengths l. The numbers are the configurations simulated for each value of α. To get the total number of configurations one has to multiply the SU (2) and SU (3) configurations by 11 and the SU (4) ones by 21.

Discussion
In this paper, we have calculated the second Rényi entropy S (2) as function of the width l of the slab-like entangled region in pure SU (N c ) Yang-Mills theory for N c = 2, 3, 4 using lattice techniques. S (2) here serves as an approximation to the entanglement entropy, which is hard to compute directly.
We have confirmed that at short distances the UV-finite part of the entanglement entropy scales proportional to the number of gluon states, N 2 c − 1. At larger values of l the entropy drops significantly, in agreement with the absence of colorful states in the low-energy confinement regime. We found that the decrease of the Rényi entropy at larger l is steeper for N c = 4 than for N c = 3, which might be an indication that in the large-N c limit the transition between small-l and large-l regimes is indeed discontinuous, as predicted by holographic models [10].
Let us note that since for finite N c phase transitions cannot happen in a gauge theory on a finite lattice, we could not expect a discontinuous transition in our data. On the other hand, since for sufficiently large lattices the large-N limit mimics the infinitevolume limit [27,28], our findings can imply two different scenarios. One possibility is that the entanglement entropy is a discontinuous function of l for any N c , but the discontinuity can only be seen by working on sufficiently large lattices and taking the proper infinite-volume limit. Another possibility is that the discontinuity appears only for infinite N c , since the predictions of [10] based on the AdS/CFT correspondence and a Hagedorn-type spectrum are strictly speaking only valid in this limit. In this case, a careful study of even larger values of N c is required. Needless to say, one also cannot exclude that the predictions of [10] are not directly applicable to pure gauge theories, but rather require supersymmetry.
An interesting observation is that the data for SU (2) gauge theory appear to be qualitatively different from that for the SU (3) and SU (4) in that the entropic Cfunction does not decrease to zero at large cut lengths, but rather exhibits a fluctuating behavior without any clear trend. While this difference might well originate in lattice artifacts, it is also possible that a different behavior of entanglement entropy in SU (2) and SU (N c > 2) gauge theories is a physical feature related to the different orders of finite-temperature deconfinement phase transitions in these theories. It is not unreasonable to conjecture that the behavior of entanglement entropy near critical cut length might be related to the behavior of thermal entropy in the vicinity of a finitetemperature phase transition, as both transitions appear to be qualitatively similar for theories with a Hagedorn-type spectrum of states [10]. In the vicinity of a second-order thermal phase transition in SU (2) gauge theory one can expect an enhancement of entanglement entropy and its fluctuations due to the divergence of correlation length [29], which is absent for first-order finite-temperature transitions in SU (N c > 2) gauge theories.
For future work, it would be of obvious interest to increase the statistics and perform a continuum limit. Going to larger N c may reveal a steeper slope and thus a more abrupt transition, which our data already hints at. However, the computational cost to accurately evaluate the integral (3.8) increases not only due to the obvious N cdependence, but also due to the need for a finer lattice of α-values at larger N c . For N c = 3, 4, this so far resulted in an additional factor of 2 Nc . This observation makes it questionable whether accurate data at larger N c will be available in the near future.
On the other hand, improving our understanding of the behavior of the entropic C-function for SU (2) Yang-Mills theory is in reach using moderate computational resources and certainly interesting to pursue.
Also the holographic calculations of the entanglement entropy can be brought closer to the lattice setup by considering corrections due to finite 't Hooft coupling. Instead of entanglement entropy, one can also calculate the Rényi entropy within holographic models [30], which might exhibit a behavior closer to what we have found numerically.