An Alternative Method for Extracting the von Neumann Entropy from Renyi Entropies

An alternative method is presented for extracting the von Neumann entropy $-\operatorname{Tr} (\rho \ln \rho)$ from $\operatorname{Tr} (\rho^n)$ for integer $n$ in a quantum system with density matrix $\rho$. Instead of relying on direct analytic continuation in $n$, the method uses a generating function $-\operatorname{Tr} \{ \rho \ln [(1-z \rho) / (1-z)] \}$ of an auxiliary complex variable $z$. The generating function has a Taylor series that is absolutely convergent within $|z|<1$, and may be analytically continued in $z$ to $z = -\infty$ where it gives the von Neumann entropy. As an example, we use the method to calculate analytically the CFT entanglement entropy of two intervals in the small cross ratio limit, reproducing a result that Calabrese et al. obtained by direct analytic continuation in $n$. Further examples are provided by numerical calculations of the entanglement entropy of two intervals for general cross ratios, and of one interval at finite temperature and finite interval length.


Introduction
The von Neumann entropy provides a quantitative measure of entanglement in a quantum system. In the special case of a thermal ensemble specified by a Hamiltonian, the von Neumann entropy reduces to the standard statistical mechanics entropy. In quantum field theory, the thermal entropy may be obtained directly from the partition function computed by standard functional integral methods, whereas the methods for calculating the entanglement entropy of a general subsystem are much less systematic, even in free field theories.
Quantum field theory does provide, however, a reasonably systematic method for evaluating the Rényi entropy [1] S ν (ρ) = 1 1 − ν ln Tr(ρ ν ) (1.1) when ν is an integer n greater than 1, by replicating the functional integral representation for ρ, the density matrix, n times. The von Neumann entropy S(ρ) = − Tr(ρ ln ρ) (1.2) is then obtained by taking the ν → 1 limit of the Rényi entropy: Before taking such a limit, note that even though the Rényi entropy S ν is well-defined for both integer and non-integer ν, it is directly computed by replicating the functional integrals only when ν is an integer n > 1. Therefore, one needs to determine the full function S ν from the set of its integer values {S n , n = 2, 3, · · · }, via some kind of analytic continuation. In practice, one achieves this analytic continuation by finding a locally holomorphic function S ν of ν ∈ C that takes the values S n = S n (ρ) for all integers n > 1 and has suitable asymptotic behaviors as ν → ∞ as required by Carlson's theorem 1 . If such a function S ν is found, Carlson's theorem guarantees its uniqueness and thus we obtain S ν (ρ) = S ν . This replica method has been applied extensively to the calculation of the entanglement entropy in many quantum systems, including quantum field theories and conformal field theories (CFTs). For example, the entanglement entropy of a single interval of length L in the vacuum state of a two-dimensional CFT on an infinite line is given by the universal formula S(ρ) = c 3 ln L ε , (1.4) where c is the central charge and ε is a UV cutoff with dimension of length.
In this paper, we present an alternative method for extracting the von Neumann entropy from the Rényi entropies S n (ρ) for integer n > 1. Our method does not rely on a direct analytic continuation in the variable n. Instead, the starting point is that once we know the traces of powers of the density matrix R n (ρ) ≡ Tr(ρ n ), (1.5) we are allowed to define a generating function of an auxiliary variable z for these R n (ρ): For a density matrix ρ, the series is absolutely convergent in the unit disc |z| < 1.
Choosing the branch cut of the logarithm to lie along the positive real axis, the function G(z; ρ) may be analytically continued from the unit disc to a holomorphic function in the cut plane C \ [1, ∞). The limit of this analytically continued function as z → −∞, which is well within the domain of holomorphicity, gives the von Neumann entropy: as may be verified directly from the definition (1.6). The existence of the analytically continued function in z beyond the unit disc may be seen from (1.6) as well. It may also be verified by rewriting the generating function in terms of a Möbius transformed variable w: Its Taylor series in powers of w is (1.10) The first few terms of the series written in terms of R n = Tr(ρ n ) are For a density matrix ρ, the Taylor series is absolutely convergent in the unit disc |w| < 1. Transforming back to the z variable, this provides explicitly the analytic continued function in Re(z) ≤ 1/2. An advantage of working with the w variable is that we may obtain the von Neumann entropy directly by taking the w → 1 limit: This may be written explicitly as a series: which provides an exact expression for evaluating the von Neumann entropy from R n = Tr(ρ n ) with integer n > 1.
In addition to extracting the von Neumann entropy, it is worth noting that R ν (ρ) = Tr(ρ ν ) for arbitrary powers ν ∈ R + may be evaluated by using similar generating functions as well, whenever the corresponding trace is convergent.
Our method is related to the resolvent method used in e.g. [4]. One advantage of our method is that it can be applied directly to numerical calculations as we will demonstrate using Eq. (1.13) shortly.
The remainder of this paper is organized as follows. In section 2 we follow the procedure described above to recover the von Neumann entropy in the simple case of a single interval, whose solution by analytic continuation in n is immediate. In section 3, we use our method to evaluate the von Neumann entropy of two intervals in the small cross ratio limit in 2d CFT starting from R n (ρ) = Tr(ρ n ), and reproduce the results of [5] for this system in the same limit. In section 4, we set up the basis for numerical calculations of the von Neumann entropy using the generating function. In section 5, we present numerical results for the entanglement entropy in the two-interval example but now for finite values of the cross ratio. In section 6, we present numerical results for the entanglement entropy of one interval in a two-dimensional CFT at finite temperature. In section 7, we analyze in detail the rate of convergence of the Taylor series for the generating function in the Möbius transformed variable w at w = 1. We conclude and comment on a few open questions in section 8.

Analytical calculation for one interval
We now apply the method to our first example: to extract the von Neumann entropy of one interval of length L in the vacuum state of a two-dimensional CFT on an infinite line, from its Rényi entropies for integer n > 1. Here c is the central charge and ε is a UV cutoff. This is an example where direct analytic continuation in n is obvious, but it is nonetheless interesting to see how our method works here. This example also serves as a starting point for more complicated examples such as the one to be analyzed in the following section.
To apply our method, we start by rewriting Eq. (2.1) as where for convenience we have defined We use these R n (ρ) values for integer n to form the generating function Combining e −ny with z n and expanding the remaining exponential in powers of y, we obtain where we have defined the sum for nonnegative integer j. We will eventually take z → −∞, so we need to find the behavior of the sum (2.7) in this limit. To do this, we use the partial fraction decomposition 1 n(n + 1) k = 1 n − k s=1 1 (n + 1) s (2.8) and perform the sum (2.7) in terms of the polylogarithm function defined by Li s (z). (2.10) To find the behavior of F j (z) as z → −∞, we need the behavior of the polylogarithm as z → −∞, which is given by the Sommerfeld expansion From this we find as z → −∞. Substituting this into Eq. (2.6), we obtain as z → −∞ Taking the z → −∞ limit, we obtain the von Neumann entropy as expected.

Analytical calculation for two intervals
We now apply our method to calculate the entanglement entropy of two disjoint intervals in the vacuum state of a two-dimensional CFT on an infinite line. Let us call the The cross ratio is then defined as We will focus on the small x limit, where R n (ρ) = Tr(ρ n ) was calculated in [5,6] and given by where is a UV cutoff, α/2 is the lowest dimension in the operator spectrum of the CFT, N is the multiplicity of the lowest-dimensional operators, and · · · denotes higher order corrections. For examples, we have N = 2 for a free boson and N = 1 for the Ising model. For notational simplicity, let us define and write the sum in Eq. (3.2) in a different but equivalent way, leading to At the leading order O(x 0 ), the form of Tr(ρ n ) is similar to that of a single interval studied in section 2, and therefore the von Neumann entropy may be extracted by the same steps, resulting in We now work at the subleading order O(x α ). As our method deals with Tr(ρ n ) in a completely linear way, we may focus on the O(x α ) term in Eq. (3.4). To extract the von Neumann entropy at this order, we define the generating function at order x α where we have removed a multiplicative factor N x 4 α for notational simplicity 2 .
Our aim is to evaluate this function and analytically continue it to z = −∞. To this end we begin by deriving an integral representation.

Integral representation forG
We first combine e −ny with z n in Eq. (3.6) and expand the remaining exponential in powers of y as in section 2, obtaining (3.7) We now use the expansion proposed in [6] for the functions where p k (α) is a polynomial in α of degree k. The radius of convergence of this Taylor expansion in u is π, where the first singularity away from u = 0 is located. Substituting u = π /(n + 1), we obtain the following sum Using a representation for the factor (π ) −2α+1 in terms of an integral over an auxiliary variable t, which is convergent for all values Re(α) > 1 2 , and obtaining the factor (π ) 2k by applying a derivative in t of order 2k, we have Carrying out the finite geometric sum over , and substituting the result into the expression (3.9), we get the following integral representation forG(z; ρ): (3.12) One verifies that the above integral representation converges absolutely for Re(α) > 1 2 . Indeed, for fixed |z| < 1 and large t, the function H k (z, t) and all its derivatives tend to a finite limit, so that exponential convergence of the integral in Eq. (3.11) is assured as t → ∞. Furthermore, the functions H k (z, t) vanish linearly at t = 0, so that the ratio H k (z, t)/(e πt − 1) and derivatives of the ratio inside the integral in Eq. (3.11) are integrable at t = 0.
We may rewrite Eq. (3.12) conveniently using F j (z) defined by Eq. (2.7): Again F j+2k (z) may be expressed in terms of the polylogarithm according to Eq. (2.10).

Analytic continuation to z → −∞
To compute the limit ofG(z; ρ) as z → −∞, we need the behavior of F j+2k (z) as z → −∞, which is given by Eq. (2.12): (3.14) In this limit, Eq. (3.13) becomes 3 As a result, the evaluation ofG(z; ρ) in this limit reduces to where the coefficients g k (α) are given by the following integral representation The integral is absolutely convergent for all Re(α) > 1 2 . To evaluate it, we shall use analytic continuation in α to complex values Re(α) > k + 1 2 such that we may integrate by parts 2k times (with vanishing boundary terms), and obtain the following expression For fixed α and k + 1 2 > Re(α), the sign of g k (α) alternates as a function of k and its magnitude grows with k. Thus the sum (3.17) is an asymptotic expansion with alternating coefficients, which we will evaluate in the next subsection using a procedure similar to Borel resummation.

Matching with previous results
To make contact with the results of Calabrese, Cardy, and Tonni in [5], we begin by using the functional relation of the Riemann ζ-function to express ζ(2α − 2k) in terms of ζ(2k − 2α + 1), and then use the standard integral representation for ζ(2k − 2α + 1), to recast g k (α) as follows: Substituting this representation into (3.17), we find (3.23) The sum over k is easily recognized in terms of the defining relation (3.8) of p k (α): Substituting this into (3.23), we find a much simplified integral representation: The integral may be simplified upon integrating by parts, and we obtain The integral is absolutely convergent for −1 < Re(α) < 0 and may be evaluated exactly for this range of parameters: Applying the reflection formula for the Γ-function we see that the pre-factor cot(πα) cancels both the pre-factor sin(πα) in (3.26) and the factor cos(πα) from the reflection formula, and we find . (3.29) Restoring the multiplicative factor N x 4 α that we ignored near the beginning of the calculation, and combining this with the leading order result (3.5), we find the von Neumann entropy of two intervals where · · · denotes higher order corrections. Note that the first term is simply S A + S B , so the remaining terms give −I AB where I AB is the mutual information. This coincides exactly with the result of [5] obtained by direct analytic continuation in n.

Setup of the numerical method
While the application of our method in various examples of conformal field theory is highly nontrivial, it turns out that the numerical application of our method provides promising approximations in a variety of complicated models. In the present section, we shall establish the numerical procedure for the next two sections. Given our generating function where we have defined the goal is to find the von Neumann entropy via the limit Similarly, for two disjoint regions A and B, we are often given but are interested in the mutual information We conjecture 4 that our method applies to the mutual information as well. In other words, we expect to be able to use the same generating function (4.1) with f (k) given by Eq. (4.4) and obtain the mutual information where ρ can be thought of as the density matrix on A ∪ B. The Möbius transformation maps z = 1 to w = ∞ and z = ±∞ to w = 1, allowing us to rewrite the generating function as (4.8) Supposing that G(z; ρ) has singularities within [z 1 , Therefore we have mapped a difficult problem of evaluating G(z; ρ) well outside its radius of convergence (which is z 1 ) to an easy one of evaluating h(1).
The radius of convergence for h(w) is w(z 2 ), so in the worst scenario of z 2 = ∞ we would evaluate h(1) at the edge of the radius of convergence. If z 2 is finite, we would evaluate h(1) as an absolutely convergent series. The value of z 2 is equal to the inverse of the smallest eigenvalue ρ min of the density matrix ρ. Therefore we reorganize the terms in (4.8) as a manifest Taylor series in w: may be determined as a linear combination of f (1), f (2), · · · , f (k). Numerically it is straightforward to evaluate (4.9) by truncating the Taylor series at some large order k max . Therefore by knowing the Rényi entropies up to order k max + 1, we can numerically estimate the von Neumann entropy with good precision.
In all the field theory examples that we will consider below, the sum (4.9) exhibits stable power-law behaviors at large k. In other words,f (k)/k ≈ Ck −p with some power p > 1 at large k. Intuitively, this is because in a quantum field theory, the smallest eigenvalue ρ min → 0 in the continuum limit and we are therefore evaluating the sum (4.9) at the edge of its radius of convergence, leading to power-law convergence. We will determine the power p numerically in each example below, and give an interpretation of this power in section 7.

Numerical studies of two intervals
Many two-dimensional CFTs may be constructed in terms of the free field theory of scalar bosons using the Coulomb gas representation and bosonization. For those theories the basic building blocks are the CFTs of a free boson with central charge c = 1, which are further distinguished by the compactification radius. In this case, the von Neumann entropies for the union of two intervals on an infinite line are characterized by the cross ratio x, as well as a universal critical exponent η, which is proportional to the square of the compactification radius.
For finite x and η in a free boson CFT, the general form of Tr(ρ n ) was derived in [6]: where is a UV cutoff and c n is a non-universal, model-dependent coefficient with c 1 = 1 [7]. In our numerical calculations below, we will choose c n = 1 for simplicity and choose a reasonable value of . 5 Note F n (x, η) is defined as for generic integers n ≥ 1. The Riemann-Siegel theta function Θ(z|Γ) is defined as where Γ is an (n − 1) × (n − 1) matrix with elements Γ rs = 2i n n−1 k=1 sin πk n β k/n cos 2πk n (r − s) , and with 2 F 1 being the hypergeometric function. Note that (5.2) is manifestly invariant under η ↔ 1/η. Currently, the analytic continuation of the von Neumann entropy for general finite x and η is not analytically known. But given our method, one can calculate the von Neumann entropy with high accuracy numerically. In the following two subsections, we will present the numerical studies of two intervals in the small x or the decompactification η → ∞ limit, where analytic perturbative expansions are available for comparison with our method. We will look at the general case of finite x, η in the third subsection.

Two intervals at small cross ratio
For general η = 1, 6 F n (x) has the following small x expansion [5,6]: where α is twice the lowest operator dimension in the CFT, N denotes its multiplicity, and · · · denotes higher-order terms. For a free boson, we have α = min[η, 1/η] and N = 2. We will numerically calculate the von Neumann entropy and confirm its small x expansion where c 1 is minus the n-derivative of c n at n = 1, which is determined by matching S A∪B to S A + S B in the limit of x 21 , x 43 x 31 , x 42 . This small x expansion of the von Neumann entropy agrees with our analytical calculation in section 3.
To set up the numerics, we model the system of two disjoint intervals A = [x 1 , x 2 ] and B = [x 3 , x 4 ] on an infinite line 7 with |x 12 | = r = |x 34 |, and set the distance between the centers of A and B to be L = 14. In this case, the cross-ratio is Note also that |x 14 | = L + r = L(1 + √ x) and |x 23 | = L − r = L(1 − √ x). Thus we may rewrite the "easy factor" in front of F n (x, η) in Eq. (5.1) in terms of x, L.
We have performed the numerical calculation with many choices of numerical constants. For example, with x = 0.25, α = η = 0.295, and 2 = 0.3, the result is S A∪B ≈ 1.216 by summing up to k max = 800, within 10 −2 of the answer S A∪B ≈ 1.224 approximated by (5.7). The sum (4.9) exhibits a stable power p ≈ 1.744, with around 10 −3 relative error. We have probed more regimes in the parameter space; see Figure 1, where we have plotted the entanglement entropies calculated by our numerical procedure compared with the analytical result given by (5.7).
Furthermore, the second order contribution in x given in [5] can be numerically evaluated, even though we are not aware of an explicit analytic continuation 8 to n ≈ 1.
In this case, we use Note that the second order terms only start contributing at n ≥ 4. Now again with x = 0.25, α = 0.295, and 2 = 0.3, the result is S A∪B = 1.105 by summing up to k max = 70. In particular, the second order correction to the entanglement entropy is negative as suggested by the last minus sign in (5.7). Since we are currently unaware of an analytical calculation at second order in x for the free boson, we do not have a way to analytically confirm our numerical results in general, but our numerical study is consistent with the holographic case [9].

Two intervals in the decompactification limit
We now consider a different limit [6] than the small x case. In the decompactification limit η → ∞, we have for each fixed value of x However, the holographic case was studied in [9][10][11]. where F k/n is defined as in Eq. (5.5). To see F n (0, η) = 1 in the η → ∞ limit, we note that both the numerator and the denominator go to infinity.
We will use the symmetry η ↔ 1/η to study the result for η 1 instead: For the decompactification regime, our numerical result will be tested against the von Neumann entropy approximated by the following expansion [6] S A∪B (η 1) S W AB + where S W AB is the von Neumann entropy computed from (5.1) without F n (x, η), and ln F z (x). (5.14) We have studied various cases; see Figure 2. For example, with x = 0.25, η = 0.295, and 2 = 0.1 (where it is best approximated [6,8]), we get S A∪B ≈ 1.584 by summing up to k max = 400, which is very close to the analytical answer S A∪B ≈ 1.608 approximated by (5.13). The sum has a stable power p ≈ 1.695, with around 10 −3 relative error. Note that even with the same choices of numerical constants, the small x expansion would not in general agree with the decompactification limit, which was already indicated in [5]. This is expected as we are not in the exact limit of either x → 0 or η → 0.  Finally, we may also study the mutual information in the decompactification limit, where we consider the following generating function that according to (4.5) should give the mutual information The mutual information is approximated in [6] by where again I W AB is the mutual information computed from (5.15) without F n (x, η). See Figure 3 for our numerical results. As an example, if we take x = 0.25, η = 0.295, and 2 = 1, we find numerically I AB ≈ 0.456 by summing up to k max = 400, within 10 −5 of (5.16) as well as the answer approximated by (5.17). The sum exhibits a stable power p ≈ 2.457, within 1% relative error.

Two intervals at finite cross ratio and compactification radius
For the most general case of finite x and η, one need to evaluate directly the Riemann-Siegel theta function in (5.1) and (5.2) numerically. In this case, an analytical expression for the von Neumann entropy is not yet known.
We start with (5.2) which we reproduce here for convenience The expression has a symmetry under η ↔ 1/η. To obtain more accurate results in our numerical procedure, one would need to evaluate higher dimensional matrices within the Riemann-Siegel theta function for specific choices of x and η. This is difficult 9 to deal with numerically for k max 1. Our approach is to use the identity (5.19) to rewrite (5.18) in the following way [6]: . (5.20) We use this formula to perform the numerical calculation. For example, with x = 0.25, η = 0.295, and 2 = 1, summing up to k max = 15 we get S A∪B ≈ 0.747.
6 Numerical studies of one interval at finite temperature and length Our next nontrivial example is a single interval at finite temperature and finite length. This example was studied for a 2D free Dirac fermion on a circle using bosonization in [14]. For such a finite system we would need to consider periodic boundary conditions for both space and imaginary time, corresponding to finite size and finite temperature, respectively. Setting the spatial size to 1, the two dimensional Euclidean theory thus lives on a torus defined by z ∼ z + 1 and z ∼ z + τ , with τ = iβ for temperature β −1 . We use to denote the length of the interval. Using Tr(ρ n ) calculated in [14], we find that the generating function is given by where is the UV cutoff 10 and ν is determined by the boundary condition for the fermion. In particular, we will study the case of ν = 3 that corresponds to the Neveu-Schwarz (NS-NS) sector. Note that η(τ ) is the Dedekind eta function defined as where q = e 2πiτ . The Jacobi theta functions θ 1 and θ 3 are defined as The exact expressions for the von Neumann entropy are only known in hightemperature and low-temperature expansions. For the high-temperature expansion, the von Neumann entropy is Note that the first term reproduces the infinite length, finite temperature von Neumann entropy [7,15,16] which is universal. Numerically, we take the choices of β = 0.9, = 0.5, and = 0.1. By summing up to k max = 700 for (6.1), we get S A ≈ 0.580, which is within 10 −3 of the analytical answer S A ≈ 0.582 from (6.4), with a stable power p ≈ 1.873, within 10 −2 relative error. For the low-temperature expansion, the von Neumann entropy is Similarly, we see that the first term reproduces the finite size, zero temperature von Neumann entropy [7,15,16] S L A = 1 3 ln 1 π sin π , (6.7) which again is universal. Numerically, we take the choices of β = 10, = 0.5, and = 0.1. By summing up to k max = 700 for (6.1), we get S A ≈ 0.385, which is within 10 −3 of the analytic result S A ≈ 0.386 from (6.6), with a stable power p ≈ 1.922, within 10 −3 relative error.
Even though the analytic expressions are known only for the high and low β regimes, we can numerically interpolate between the two regimes using the generating function (see Figure 4).  Figure 4. Numerical studies of a single interval A at finite temperature and length. We fix the interval length = 0.5 and the cutoff = 0.1, and plot S A as a function of β. The blue curve is plotted by the high-temperature expansion (6.4). The green curve is plotted by the low-temperature expansion (6.6). Each point is calculated by our numerical procedure up to k max = 200. One can see that the high-temperature expansion fits well in the plot, but the deviation with our numerical method increases slowly for larger β not shown here.

Power-law convergence of the generating function
The various field theory examples that we have considered so far all exhibit power-law convergence at w = 1 of the series (4.9). In other words, we havẽ at large k with a power p that depends on each specific example. This leads to a natural question: what type of eigenvalue distributions for a density matrix ρ would exhibit such power-law behaviors?
To answer this question, we define the eigenvalue distribution P (x) so that the number of eigenvalues within a small range [x, x + dx] is P (x)dx. We may rewrite Tr(ρ n ) in terms of P (x): In particular, Tr ρ = 1 means In the case of a finite-dimensional density matrix with eigenvalues {ρ i }, we have P (x) = i δ(x − ρ i ). Using Eq. (7.2), we may write Eq. (4.10) as We would like to find the behavior of this integral for large k. We expect the integral to be dominated by small x when k is large, due to the presence of (1 − x) k . Let us therefore consider a power-law ansatz for the behavior of the eigenvalue distribution P (x) near zero eigenvalue: Convergence of the integral in Eq. (7.3) then requires γ > −2. Substituting this into Eq. (7.4), we find at large k This is indeed a power law for large k. Comparing it with our definition of the power p in Eq. (7.1), we find p = γ + 3. (7.7) As we mentioned previously, Eq. (7.3) requires γ > −2, giving p > 1, which is precisely the necessary and sufficient condition for the sum (4.9) to converge at w = 1. Eq. (7.7) describes how the power-law behavior of the sum (4.9) is related to the small eigenvalue behavior of the eigenvalue distribution P (x). Therefore, using the power p that we have determined numerically in the field theory examples studied in previous sections, we may now predict that their eigenvalue distribution P (x) must scale like x p−3 for small eigenvalue x.
As we alluded to briefly in section 4, this power-law convergence is closely related to the infinite-dimensional Hilbert space of quantum field theories. Let us use {ρ i } to denote the eigenvalues of the density matrix. The generating function (1.6) becomes It has branch cuts along z ∈ [1, 1/ρ min ] where ρ min is the smallest eigenvalue of the density matrix ρ. After the Möbius transformation to w, the branch cuts are along w ∈ [ 1 1−ρ min , ∞). A necessary condition for the series (4.9) to be power-law convergent at w = 1 is that its radius of convergence must be 1. This means 1 1−ρ min = 1, or ρ min = 0. 11 This is true for density matrices in continuous quantum field theories with infinite-dimensional Hilbert spaces.

Discussion
In this paper we have shown that the von Neumann entropy may be obtained by assembling the traces Tr(ρ n ) for all positive integers n into a generating function of an auxiliary complex parameter z, analytically continuing in z, and then taking the limit as z → −∞. Our construction demonstrates that the analytic continuation in z exists when ρ is a density matrix. We showed how the procedure may be carried out analytically for the cases of one interval and of two colinear intervals in certain limits, and we demonstrated that a simple variant of the method also leads to numerical evaluations of the von Neumann entropy in a practical, reliable way.
Many open questions and future directions for investigation remain. Most of the cases that we studied were selected because they are simple enough to exhibit the method clearly, but complicated enough so that the standard replica "analytic continuation" in n is not automatic. Firstly, it is urgent to see whether our method can analytically produce the von Neumann entropy in cases that have eluded direct analytic continuation in n. One such case is the example of two intervals at finite cross ratio for a free boson with a general compactification radius. Another interesting example involves the contributions of replica non-symmetric saddle point solutions in holography, which are difficult to analytically continue in n but have been shown in [17] to give important enhanced corrections to the Ryu-Takayanagi formula [18,19] at holographic entanglement phase transitions.
Secondly, our method appears to have a broader regime of applicability than what might be expected from our derivation -in particular, the method often applies to parts of Tr(ρ n ) with "easy factors" stripped off (as in footnote 2) and we also conjecture that it applies to the mutual information (as in section 4). It would be very interesting to understand precisely how broadly our method applies and why.
Finally, an alternative starting point for the replica method, which recently appeared in [3], is from the functions Tr(ρ 1/n ) instead of Tr(ρ n ). It will be interesting to see whether our method can be adapted to handle such cases as well.