Adiabatic corrections to holographic entanglement in thermofield doubles and confining ground states

We study entanglement in states of holographic CFTs defined by Euclidean path integrals over geometries with slowly varying metrics. In particular, our CFT space-times have S1 fibers whose size b varies along one direction (x) of an ℝd−1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathbb{R}}^{{{}^d}^{-1}} $$\end{document} base. Such examples respect an ℝd−2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathbb{R}}^{{{}^d}^{-2}} $$\end{document} Euclidean symmetry. Treating the S1 direction as time leads to a thermofield double state on a spacetime with adiabatically varying redshift, while treating another direction as time leads to a confining ground state with slowly varying confinement scale. In both contexts the entropy of slab-shaped regions defined by |x − x0| ≤ L exhibits well-known phase transitions at length scales L = Lcrit characterizing the CFT entanglements. For the thermofield double, the numerical coefficients governing the effect of variations in b(x) on the transition are surprisingly small and exhibit an interesting change of sign: gradients reduce Lcrit for d ≤ 3 but increase Lcrit for d ≥ 4. This means that, while for general L > Lcrit they significantly increase the mutual information of opposing slabs as one would expect, for d ≥ 4 gradients cause a small decrease near the phase transition. In contrast, for the confining ground states gradients always decrease Lcrit, with the effect becoming more pronounced in higher dimensions.


Introduction
Entanglement is a fundamental property of quantum systems. Studying this entanglement can provide insights into the nature of quantum states, and in particular into the scale of their correlations. In the holographic context, entanglement of the dual CFT is of particular interest through its association with the Einstein-Rosen bridges of black holes [1] and perhaps more generally [2][3][4] with the emergence of bulk spacetime.
Our goal here is to generalize the analysis of holographic entanglement away from the commonly-considered highly symmetric systems. For d = 2 CFTs, much can be done exactly using conformal transformations. This fact lies behind the recent analysis [5] of the CFT states dual to asymptotically-AdS 3 mutli-boundary vacuum wormholes. In particular, it was understood there that such states admit a simple description at high temperatures where the state can be well-approximated by a thermofield double (TFD) over most of the CFT spacetime, perhaps with adiabatic variations from one point to another. While a full analysis comparable to [5] is difficult in higher dimensions, we show below that computations of entanglement in spatialy-varying holographic TFDs remains tractable in the adiabatic limit.

JHEP09(2016)058
We also investigate how entanglement in ground states of (d−1)-dimensional confining theories is affected by slow variations of the confinement scale. The particular class of confining theories we consider are those given by compactifying a d-dimensional holographic CFT on an S 1 as in [6]. Such CFT ground states are related to the above thermofield doubles, as both are given by cutting open Euclidean path integrals over geometries with S 1 × R d−1 topology. Roughly speaking, the thermofield double states are given by cutting open the S 1 factor, while ground states of confining theories are given by cutting open a direction of the R d−1 . The particular path integrals considered here will involve warped products of the S 1 over R d−1 in which the size b of the S 1 varies slowly. This gives in the first interpretation TFD states in spacetimes with spatially varying redshift, and in the second ground states of confining theories in which the confinement scale varies with position.
Since we are interested in holographic field theories, in all cases we will work directly with the dual gravitational description. Our CFT path integrals are then interpreted as integrals over all (d + 1)-dimensional asymptotically locally Anti-de Sitter (AlAdS) spacetimes with boundary geometries as above. Section 2 begins below by reviewing the Euclidean bulk geometries recently constructed in [7] that are expected to describe the dominant AlAdS saddle points. For simplicity, we allow b to vary only along one Cartesian direction of the R d−1 space. While such solutions can be constructed by Wick rotating the standard fluid-gravity correspondence [8][9][10] in the presence of a time-translation Killing field and an appropriate regularity condition at the bifurcation surface, it is more natural to follow [7] and use the U(1) symmetry to develop a related but different expansion based on standard Schwarzschild-like coordinates rather than the ingoing Eddington-Finkelstein black hole coordinates of [8][9][10].
We then proceed to compute holographic entanglement. Section 3 pursues the thermofield-double interpretation and computes the effect of varying b on the Ryu-Takayanagi (RT) entropies of slabs of thickness 2L that preserve R d−2 Euclidean symmetry on a surface fixed by a reflection of the S 1 . We include both the case of slabs contained in a single copy of the CFT and that of pairs of diametrically opposed slabs in each of the two CFTs. We thus also compute the effect of varying b on the mutual information in opposing slabs and on the critical value L crit of L at which the mutual information becomes non-zero. Section 4 then studies the effect on RT entropies for analogous slabs with S 1 × R d−3 symmetry on a surface fixed by reflecting one direction in the R d−1 . Here the interesting feature is the effect on the value L crit at which the entangling surface changes topology from connected (S 1 × [0, 1] × R d−3 ) to disconnected (two copies of the (d − 1)disk). Readers focused on the final results may wish to jump to sections 3.3 and 4.3 where the phase transitions are discussed. We close with some final discussion in section 5. The special case d = 2 is treated analytically in appendix A, and we discuss some estimation of the numerical uncertainty in appendix B.

JHEP09(2016)058 2 Preliminaries
We wish to describe holographic entanglement in CFT states defined by path integrals over geometries with topology S 1 × R d−1 and metrics of the form We take θ to have b-independent period 2π. The relevant states are constructed by slicing open the path integral along a co-dimension one surface that we identify as τ = 0 for some Euclidean time coordinate τ . To have a good translation to Lorentz signature, we require a Z 2 reflection symmetry τ → −τ . One natural choice is to take τ = θ, in which case we in fact slice the path integral along the pair of surfaces θ = 0, θ = π. The result is an entangled state on a pair of CFTs which gives an adiabatic generalization of the well-known thermofield double state. The exact timetranslation symmetry means that the state is in thermal equilibrium when viewed from the perspective of either CFT alone. However, after Wick rotation to Lorentz signature the x-dependent metric factor g θθ means that the state lives in a spacetime with x-dependent gravitational redshift. This equilibrium thus requires any local notion of temperature (such as that defined by the inverse Euclidean period) to be x-dependent as well. This interpretation is equally valid in the special case d = 2 in which there are no y directions.
For d ≥ 3, there is a second interpretation given by choosing τ to be some y direction (say, y 1 ), so that our CFT lives on a spacetime with a compact spatial S 1 . States of this theory are constructed by slicing the path integral along y 1 = 0. For small b one may Kaluza-Klein reduce on this S 1 . And as discussed in [6], one expects the result to exhibit confinement with a scale set by b. So when b varies, one may think of the result as a confining theory with a position-dependent confinement scale.
But with either interpretation, so long as b varies slowly reasoning analogous to that of [6] implies the bulk path gravitational integral with boundary conditions given by (2.1) to be dominated by a Euclidean solution to Einstein's equation in which the S 1 factor pinches off in the bulk; i.e., there will be a Killing field ∂ θ that generates a U(1) isometry with a fixed-point set of topology R d−1 .
When the function b(x) varies slowly, the construction of such solutions may be organized in a derivative expansion. Here we write b = b( x) for some small parameter . The details of this expansion were recently described in [7], where it was argued that for slowly-varying b(x) the solution should be well-approximated by the zero-order ansatz 2) where we take θ to have period 2π for all profiles b(x). For the case b = constant, the ansatz (2.2) gives the metric on the Euclidean planar AdS-Schwarzschild black hole (or, equivalently, on the Euclidean AdS soliton). The full metric is then taken to be of the form 3)

JHEP09(2016)058
where the corrections g (n) AB are determined by solving Einstein's equation with appropriate boundary conditions at each order in an adiabatic expansion and x A = (z, x, y i , θ) ranges over all bulk coordinates and similarly for x B . As shown in [7], the O( ) correction g (1) AB vanishes and, writing g y i y j = g yy δ ij , the O( 2 ) correction is of the form Here the notation makes explicit all dependence on b(x); there can be no further implicit dependence hidden in form of the coefficient functions g These coefficient functions were evaluated numerically in [7] with boundary conditions that ensure that the boundary metric remains (2.1) and that the spacetime remains regular at the fixed point set of ∂ θ (with the period of θ taken to be 2π independent of b(x)).
Below, we use the results of [7] to calculate O( 2 ) corrections to the holographic entanglement entropy. We also make use of two further results from [7]. The first is that, for d > 2, in the adiabatic expansion the Fefferman-Graham representation of our metrics takes the form The special case d = 2 is treated in appendix A. The second is that near the fixed point set of ∂ θ the metric takes the form functions of X alone, in terms of coordinates X, R that satisfy The key point of (2.6) is that it ensures the desired regularity at R = 0 (where ∂ θ = 0). In terms of the Fefferman -Graham coordinates this set is described by z =b wherẽ This is the black hole horizon for the adiabatic thermofield double interpretation and the IR floor for the confining one.

JHEP09(2016)058 3 Adiabatic thermofield doubles
We begin with the adiabatic thermal field double (ATFD) states defined by slicing our CFT path integral along the surfaces θ = 0, θ = π fixed by the reflection symmetry θ → −θ. It is convenient to denote the union of these two surfaces by C CFT . A slight generalization of the Ryu-Takayangi proposal [11,12] then states that the von Neumann entropy of the CFT in some region R CFT ⊂ C CFT can be computed as follows. First, find the dominant saddle for the corresponding bulk path integral. One expects it to be invariant under a corresponding reflection, and that this reflection leaves fixed a co-dimension one surface that we may call C bulk . Now find the minimal-area surface Σ within C bulk that i) intersects the asymptotically AdS boundary on a set corresponding to the boundary ∂R CFT of R CFT and ii) is homologous to R CFT within C bulk [13,14].
Since the Lewkowycz-Maldacena argument [15] for the Ryu-Takayanagi proposal applies equally well to this generalization, we shall use it freely below. 1 We also note that the above prescription is equivalent to using the the covariant Hubeny-Rangamani-Takayanagi conjecture [17] in the Wick-rotated Lorentz-signature solution. 2 For simplicity, we consider slab-shaped regions R CFT defined by conditions of the form |x − x 0 | ≤ L, perhaps also restricted to one of the two boundaries (θ = 0 or θ = π). The symmetries then reduce the problem of finding the minimal surface to studying curves in the z, x plane, with the area being proportional to the volume of the y directions. For purposes of displaying a finite result we take the y coordinates to range over a torus of finite volume V . Since we are interested in the decompactified limit, we will always assume each cycle of the y-torus to have length much larger than both b and L. In particular, we assume that the dominant bulk saddle will continue to be given by (2.2).
A technical issue is that the area nevertheless remains infinite due to the divergence of the metric (2.2) at z = 0. As usual, we must renormalize this quantity in order to present finite results. Thus we define where A bare (z 0 ) is the area of the part of the surface with z > z 0 and where there is one counter-term contribution A ct (z 0 ) for each boundary of the minimal surface Σ. The general theory of such divergences is explained in [19], which shows that when the bulk is described by pure Einstein-Hilbert gravity (with no additional matter fields) one may use counter-terms determined by the boundary metric alone, 3 though these generally involve both the induced geometry on ∂Σ, the extrinsic curvature of ∂Σ [21,22], and even 1 See [16] for a discussion of the homology constraint in the context of the Lewkowycz-Maldacena argument. 2 We thank Veronika Hubeny for pointing out that this follows from the maximin construction of [18].
Since the RT surface is minimal on the Cauchy surface C bulk , its area can be no larger than that of the maximin surface. But the time-reversal symmetry means that the RT surface is also an extremal surface in the full spacetime. It can therefore have area no smaller than the maximin surface, as the latter agrees with the area of the smallest extremal surface. 3 Interestingly, this is not true in general; see [20].

JHEP09(2016)058
derivatives of such extrinsic curvatures [23] in high enough dimensions. See also [24] for a recent discussion of such counter-terms and their relation to [15]. To find a useful explicit form for our A ct (z 0 ) , we first write the area functional as for any parameter λ along the associated curve in the z, x plane. Near z = 0 it is useful to set λ = z and assume an adiabatic expansion of the form The behavior of x (0) near z = 0 is determined by the minimal surface equation of motion at order 0 . This may be written Equation (3.5) admits a power series solution of the form Indeed, the result takes the form (3.6) in any metric having the same non-zero coefficients in its Fefferman-Graham expansion. Since g AB = 0, at order the ansatz (2.2) continues to give the full metric. Noting that the endpoint conditions x(z = 0) = x 0 ±L are independent of then also gives So near z = 0 the area density (3.3) becomes as any factors x (0) (z) or x (1) (z) are of order z d and give corrections that vanish as z → 0.
Combining the Fefferman-Graham expansion of the second order metric correction (2.5) with the results above we find so we may choose

JHEP09(2016)058
There are no explicit O( ) counter-terms since g (1) AB = 0. One may check that this choice of counterterms precisely implements the covariant counterterm prescription of [24] to O( 2 ). Following this prescription, the counterterms in d = 4 will include a logarithmic as well as a constant piece, and in d = 2 we only have the logarithmic piece. These terms are given by where no factor of V appears in d = 2 because there are no y-directions. For d = 3, the second counter-term in (3.11) vanishes; we nevertheless find that including it in the manner explained below improves the convergence of our numerics.
In practice, we find it convenient to renormalize in the following way.
for any z max . In particular, we can take z max to be the maximal value of z on our bulk extremal surface. The renormalized area (3.1) can then be written The integral in the second line now converges, and is more stable to compute numerically. The price we pay is having to add the constant term involving z max . For d = 3, we find that including the second (vanishing!) counter-term in (3.11) in this way improves our numerical convergence. This appears to be due to the fact that we perform these integrals by changing variables to integrate over x instead of z, and that the above renormalization removes an (integrable) singularity in the integrand that arises from the associated factor of z (x).
We are now ready to compute the entropies of our slabs |x − x 0 | ≤ L. For slabs contained in a single boundary, we know on general grounds that the minimal surface will remain close to the conformal boundary when L b while for L b it will track the horizon closely over almost all of the interval |x − x 0 | ≤ L. The transition between these behaviors is smooth. But if we take our slab to contain the regions |x − x 0 | ≤ L on both the θ = 0 and θ = π boundaries one finds a well-known phase transition [17,[25][26][27] when passing from the regime L b 0 to the regime L b. In the former case, the minimal surface consists of two copies of that found in the single-boundary case. In the latter case the minimal surface again has two connected components, but each component then stretches from θ = 0 to θ = π while remaining localized near x = x 0 ± L. This is the only context in which the minimal surface reaches or passes through the fixed point set of ∂ θ . In each case we find the general solution numerically below and compare it with analytic approximations for L b and b( b ) −1 L b. We also provide results for the case The effect on the phase transition itself is analyzed in section 3.3.

Entropy on a single boundary
We begin with connected slab-shaped regions R CFT of width 2L lying in a single boundary.
For generic values of the parameters, numerical calculations are required to find the extremal surface. But certain limiting behaviors can be studied analytically. We treat these cases first, and then compare the results with numerical studies of the general case. In the rest of this section, we set x 0 = 0 without loss of generality.

Analytically tractable limits
Our first special case will be the large L limit, as the fact that the minimal surface closely tracks the horizon in this regime makes it particularly easy to study. To leading order in L, the renormalized area is just the horizon area in the region |x| ≤ L. Using the induced metric on the horizon found in [7] gives yy z=b where the . . . represent terms that do not grow with L when b remains bounded. For L larger than or comparable to b/( b ), nothing more can be said without choosing to simplify (3.14).
ren + . . . , we find where ∼ indicates that we have found only the leading behavior for L b 0 . Here we were able to obtain an analytic expression at order 2 because the L 3 term comes only from the O( 2 ) term in (3.15) and thus can involve the metric only at order 0 as given by (2.2).
For L b 0 the minimal surfaces will be confined to z b 0 , so we can estimate their area by truncating the Fefferman -Graham expansion (2.5) for the metric to some order in z. The Fefferman -Graham expansion for d = 2 has a non-trivial contribution from the boundary stress tensor at order z 2 , so we treat this case separately in appendix A.
Consulting the expansion (2.5), we see that to zeroth order in the adiabatic expansion we have Poincaré AdS d+1 . So for d > 2 we find This leading term reproduces the standard result for slabs in Poincaré AdS d+1 as derived in [11].

JHEP09(2016)058
Since dθ = 0 on the surface of time reflection symmetry, the truncated induced metric (2.5) depends on b only at order 2 and there can be no O( ) correction to the minimal surface or its area. And the fact that the zero-order surface is minimal means that there is no correction at order O( 2 ) from the second-order displacement of the surface within the zeroth-order spacetime. Thus the only O( 2 ) contribution comes from evaluating the change in the area along the zeroth-order minimal surface that comes from including the O( 2 ) parts of (2.5). This correction can be computed from the integral representation of the hypergeometric function found in equation (15.6.3) of [28] and yields

Numerics and comparisons
We now consider general values of L 1/( b ). This allows us to again use (3.15) so that the results can depend only on the parameters b 0 , b 0 , and b 0 . For d = 2, 4 we write where the function form of A(L/b 0 ) is determined only by dimensionless combinations of b and its derivatives. For d = 2 and d = 4 it is useful to subtract the logarithmic dependence on coming from the regularization scheme (3.11) and write We may then use the adiabatic expansion to write Now, the correction A (1) (L/b 0 ) would have to be proportional to the first-order adiabatic parameter b 0 /b 0 . But the sign of this parameter changes under x → −x whereas the area must be invariant. So there can be no correction at this order. We thus consider only the second order corrections, which must be linear in the two dimensionless second-order adiabatic parameters (b 0 ) 2 and b 0 b 0 ; i.e., we have      Even at order 0 we require numerics to solve for the surface that extremizes the area (3.3). We use the Newton-Raphson method outlined in [29]. Figure 1 shows the solution for z (0) (x/L)/b 0 with 2 ≤ d ≤ 7 and various interval sizes. Results for the zeroth order area are shown in figure 2.
Computing the second order change in area (3.23) requires only knowledge of the surface to O( ). This is because the order-zero surface is minimal, so changes in the area computed with with zeroth order metric are quadratic in changes of the surface. The first order equation of motion is complicated, but is straightforward to work out and can be solved numerically by the same techniques as at order 0 . Results for z (1)     A ren , this correction takes the form where each line corresponds to one of the above three contributions described above. Numerical results are shown in figure 4. As a check on our numerics, we now compare with the analytic expressions of section (3.1.1). We first consider the case b/( b ) L b 0 . At order 0 we numerically compute b 0 A (0) /L for large L/b 0 , while at order 2 we compute b 3 0 A (2) /L 3 . Results are tabulated in table 1 which shows agreement with (3.16).
Turning now to the case L b 0 , we have verified that the coefficient of A (2) proportional to b 0 2 vanishes quadratically as L b 0 , and we may also numerically compute the JHEP09(2016)058

Entropy for pairs of diametrically opposed slabs
We now consider the entropy of a pair of corresponding slabs on opposing boundaries. Both slabs are defined by |x − x 0 | ≤ L, but one lies at θ = 0 while the the other lies at θ = π. Without loss of generality we again set x 0 = 0 in this section. As in [25,27], for L b the minimal surface will be simply two copies of the one found in section 3.1, so that the mutual information between these two slabs vanishes. But for L b the minimal surface represents a different phase, again having two disconnected pieces but now with each localized near x = ±L. Here the slabs share non-zero mutual information I. In this phase the entropy is independent of L and depends only on the local behavior of b(x) near x = ±L. Note that the contribution from each surface is just the entropy one would compute for a pair of half-spaces, both defined by x > L (or x < −L) but on opposite boundaries. For simplicity we thus focus on this 'half-TFD' entropy below. All quantities associated with the half-TFD problem will be marked with hats (ˆ) to distinguish them from the corresponding quantities of section 3.1.
As before, computing the area to order 2 requires only knowledge of the entangling surface to first order. It thus suffices to writê At zeroth order the entangling surface relevant to this half-TFD problem lies at preciselŷ x (0) (z) = ±L and extends from one boundary to the other, passing through to the horizon. The total area at this order may be computed analytically and we find for the half-TFD problem with 2 ≤ d ≤ 7, with bb evaluated at x = ±L. In the large d limit, one may show analytically that this function vanishes everywhere except at the horizon.
At first order we proceed numerically, withx (1) (z) satisfying the first order equation of motion We simplify the analysis by using the symmetry that relates our two boundaries. We thus compute the area for a surface extending from one boundary to the horizon and multiply by 2. The boundary conditions are thatx = ±L at z = 0 and that dx dR = 0 at the horizon R = 0, where R is the regular coordinate associated with (2.6). But since it is convenient to work in terms of the original z coordinate, we note that to order this is equivalent to imposing the boundary conditionx The second order area correction now has an additional contribution due to the O( 2 ) shift in the endpoint of the minimal surface. This contribution can be computed analytically and the full second order shift is given bŷ tt z=b (3.29) whereÃ (2) ren includes the area of only the part of the surface with z ≤ b. Note that (3.29) depends on L only through evaluating b (and its derivatives) at x = ±L. We compute (3.29) numerically. Results are displayed in table 3 in terms of dimensionless coefficients defined bŷ  Table 3. The coefficientsÂ (b 2 ) andÂ (b b ) for the half-TFD problem for 2 ≤ d ≤ 7. The numerical precision is estimated by comparing results for 100 and 150 lattice points, giving better than one part in 10 −10 .
Here b, b , b are evaluated at x = ±L. For d = 2, 4 we use analogous notation but with logarithmic subtractions as in (3.21).

Phase transition
We now analyze the transition between the I = 0 and I > 0 phases for a pair for |x−x 0 | ≤ L slabs on opposite boundaries. In particular, we compute the effect of inhomogeneities on the critical length L crit .
For this purpose, we should compare twice the area of the entangling surface for a slab |x| ≤ L on a single boundary with that of the sum of the surfaces for the half-TFD problems at x = ±L. The phase transition will occur when L is of order b, so at small we have L b/( b ) and we may expand b(±L) in (3.30) in a Taylor series. At order 0 , the surfaces at x = ±L have equal area, so we can determine the zeroth order value of L crit by comparing (twice) the numerical value of (3.26) for b = b 0 with (twice what is shown in) figure 2. Results are displayed in table 4.
As discussed in section 3.1, the first order correction to the area of the connected surface vanishes. For the disconnected surfaces, we have a first order correction from expanding (3.30). But this correction is proportional to x b 0 , so the corresponding contributions cancel between the surfaces at x = ±L; there can be no change in L crit at first order.
At second order, we can write L crit = L Here it is useful to note that Taylor expandingÂ ren | x=±Lc and then performing our adiabatic expansion giveŝ ren | x=0 +O( 4 ).
(3.33) Table 4 displays numerical results for 2 ≤ d ≤ 7 in terms of the coefficients defined by In addition, figure 6 shows result for the mutual information between the slabs using the notationÎ We find to second order thatÎ has an interesting dependence on dimension. First althoughÎ (b 0 2 ) is positive for most L > L crit , for d ≥ 4 it becomes slightly negative near L crit . As a result, a non-zero b 0 causes L crit to increase for d ≥ 4 and decrease for d = 2, 3. The effect of second derivatives depends on dimension as well: a positive b 0 increases L crit for 2 ≤ d ≤ 4 but decreases L crit for 5 ≤ d ≤ 7. For d = 2 the above behavior is derived analytically in appendix A; it would be interesting to develop an analytic understanding of the higher dimensional results as well. Due to the many interesting features in this data, we take extra care to understand the convergence of our numerics in appendix B.

States of confining theories
We now turn to the second interpretation in which our path integral computes the ground state of a confining gauge theories on the surface y 1 = 0. This necessarily restricts our discussion to d ≥ 3. We again consider slabs |x−x 0 | ≤ L. As in section 3.2, there are two possible phases for the minimal surface. For L b the minimal surface is connected and does not reach R = 0. But there is also another local extremum of the area given by a disconnected surface that consists of two disks, each localized near x − x 0 = ±L. At small L the disconnected surface has larger area, though increasing L leads to a phase transition at which the disconnected surface becomes minimal. Interestingly, at still larger values of L the connected extremum becomes singular and ceases to exist. The two phases are shown in figure 7 and will be studied in sections 4.1 and 4.2 below.
The general feature that the entanglement becomes independent of L at large L is to be expected in confining theories, as they have finite correlation lengths. But the sharp phase transition seen here is a feature of large N [30,31].
Below, we find it useful to write whereṼ is the volume of a (d−3) torus that we use to regularize the y 2 , . . . , y d−2 directions. To compute the entropy, we must as usual find the minimal surface to O( ). We will also need the explicit counterterms that renormalizing the area functional to second order. The computations are analogous to those in section 3, though now the minimal surface  equations lead to the asymptotic expansion where x B is the point of intersection with the boundary. Inserting (4.2) into (4.1) gives so for d > 4 we may take

In lower dimensions we have
where the counterterms again match the covariant prescription of [24], whose details we have again used to fix the z-independent terms for d = 4. We can now compute the area of the minimal surface for the regimes L b and L b and study the phase transition between connected and disconnected topologies. Additionally, without loss of generality we set x 0 = 0 for the rest of this section.

Narrow slabs
We begin with the regime L b, where the entropy will be given by the connected surface [30,31]. The computations proceed much as in section 3.1, though we are no longer able to obtain analytic results for the second order area in the large and small L limits. Indeed, this phase fails to exist at sufficiently large L, while for the small L limit the first order correction z (1) (x) must be computed numerically even in the approximate geometry (2.5). However, the expansion (2.5) does require the leading small L behavior of A (2) ren to be of order L 4−d . As a test of our numerics, we compare below the coefficient of L 4−d computed using the full metric against that computed using the truncated metric (2.5). At zeroth order we can compare against an analytic prediction, as at this order (2.5) is just Poincaré AdS d+1 and θ acts just like a y-coordinate with period 2πα d b. As a result, the area is given by (3.17) with V = 2πα d bṼ .
As in section 3.1, we consider the case L b/( b ) so to order 2 the inhomogeneities are described by b 0 , b 0 , and b 0 . We state our numerical results for the connected area in terms of the dimensionless function A c (L/b 0 ) defined for d = 4 by where the subscript c will denote quantities associated with the connected entangling surface. For d = 4 it is useful to explicitly remove the log( ) dependence introduced by our regularization scheme. We therefore write As before, we use the adiabatic expansion to write where symmetry under x → −x again requires the first order correction to vanish. Numerical results are displayed in figure 8. As a check on our numerics, we extract lim L→0 L d−2 A ren and compare in table 5 with the same coefficients as determined by approximating the metric to O(z 2 ) in the Fefferman -Graham expansion (2.5).

Wide slabs
For L b, the entangling surface is given by two disconnected disks each localized near x = ±L. As in section 3.2, the entropy depends on L only through the local behavior of b(x) near x = ±L. Furthermore, the contribution from each surface is just the entropy one would compute for the corresponding half-space x > L or x < −L. For simplicity we thus focus below on this notion of 'half space entropy' and choose R CFT to be the region x > ±L. Note that our geometry ends at z =b, with the extremal surface obeying the boundary condition of regularity (3.28).
Approx.  Table 5. Comparison of the numerically computed L b 0 scaling of A(L/b 0 ) for 3 ≤ d ≤ 7 from figure 8 (left columns) with that determined by truncating (2.5) at order z 2 (right columns, with "Pred." and "Approx." referring to analytic and numerical results respectively). The numerical precision is shown when it falls below three significant figures, estimated by comparing results for 100 and 150 lattice points and for fitting different ranges of L depending on the dimension.  Table 7. Numerical values of L crit /b 0 and the coefficients L (b 0 2 ) and L (b0b 0 ) from (3.33) for the RT phase transition for slabs |x| ≤ L in our confined ground state with 3 ≤ d ≤ 7. The numerical precision is shown when it falls below three figures, estimated by comparing results for 100 and 150 lattice points.

JHEP09(2016)058
The detailed computations are much as in section 3.2, so we simply display the results. The area of the disconnected surface can be written in terms of the dimensionless functions described in (4.6) and (4.8) after replacing A c (L/b 0 ) with A d . We compute the zeroth order coefficients analytically, but the second order coefficients require numerics. Half-space entropy results for 3 ≤ d ≤ 7 are tabulated in table 6 using our by-now standard notation.

Phase transition
Finally, we turn to the effect of adiabatic variations on the critical value L crit at which the dominant phase becomes disconnected. As in section 3.3, we do so by comparing the area of the connected surface (figure 8) with the area of the disconnected surface evaluated at x = ±L (table 6). Since the phase transition occurs at L b/( b ), we again expand b(x) in a Taylor's series to compute A d . The second-order coefficients of of L crit are again given by (3.33) with the replacements 2A ren → A c ,Â ren → A d . We determine L crit numerically to second order, and display these results in table 7 using the notation of (3.34).

Discussion
In the above work, we computed the leading (second order) effects of inhomogeneities on the holographic entropy of slab-shaped regions defined by |x−x 0 | ≤ L. We studied thermofielddouble states on spacetimes where the redshift changes slowly with position, and the ground JHEP09(2016)058 states of certain confining theories with corresponding slow changes in the confinement scale. In each case, we studied the effect on the length scale L crit associated with a Ryu-Takayanagi phase transition. Most of our results were numerical, though the special case d = 2 (AdS 3 ) was treated analytically in appendix A. In higher dimensions, some analytic results were also available in special limits and were used to check our numerics.
For the thermofield double, L crit is a measure of the non-locality of entanglements between opposite CFTs. The second-order coefficients (table 4) governing the response of L crit to inhomogeneities turn out to be numerical small. Some insight as to why is provided by the analytic d = 2 treatment of appendix A, which shows these coefficients to be proportional to (L crit /b) 3 . So the small values of L crit /b lead to even smaller coefficients The coefficients shown in table 4 display highly non-trivial structure with respect to the dimension d. For d ≤ 3, gradients decrease L crit , while they increase L crit for d ≥ 4. This remains true whether one studies the local response to b 0 or the average change over all x. The former is precisely the sign of L (b 0 2 ) in table 4. But averaging L (2) crit over x allows one to use either periodic boundary conditions or b → constant as x → ±∞ to integrate b 2 b by parts, giving a positive-definite quantity multiplied by (L (b 0 2 ) − 2L (b 0 b 0 ) ). It turns out that both change sign between d = 3 and d = 4. Interestingly, it is the large d behavior that corresponds to the naive expectation that that the response is given by averaging b(x) over a scale |x−x 0 | b, as such averaging would decrease L crit near a maximum of b(x) and thus require L (b 0 b 0 ) < 0. This is the opposite sign to that found analytically for d = 2 in appendix A.
One also notes that the coefficients L (b 0 b 0 ) are not monotonic with d, but appear to have a local minimum near d = 6. In contrast, L (b 0 2 ) appears to be monotonic in d but is also highly non-uniform. In particular, while most cases exhibit a clear increase in value with d, the coefficients for d = 5 and d = 6 are remarkably close. The in-depth analysis of numerical convergence in appendix B appear to confirm that these features are real and are not just numerical artifacts. It would be useful to have an analytic understanding of these dimension-dependent features; the large d limit may be worth particular study.
In contrast, the response of our confining ground states is both larger and more uniform in d; table 7 shows no changes of signs. It is nevertheless interesting that gradients -either local or averaged -always decrease L crit . This is naturally understood as a corresponding decrease in the length scale characterizing confinement. But comparing our results with [7] challenges this interpretation. For d ≤ 5, [7] found that the gradients decrease the tension of flux tubes aligned in their direction, while the increase of tension one would expect from a decrease in the confinement length scale occurred only for d ≥ 6. Furthermore, for d > 3 it found that gradients always raised the negative energy of the confining ground state -a result naturally associated with a larger confinement length scale. The main conclusion appears to be that confinement is not generally characterized by a single scale, but that changes in different confinement-related phenomenon under small perturbations are often uncorrelated. It would be interesting to develop more analytic understanding of such effects, and also to determine to what extent our results apply to other systems with spatially-varying confinement scale such as those that might be constructed in a condensed matter laboratory.

JHEP09(2016)058
coming from integrating the zeroth order surface over the region z ∈ [b, z H ] precisely cancels the second order contribution associated with the first-order shift of extremal surface within the zeroth order background. As these were the only possible contributions to this order, in agreement with our numerics we find that the full second order contribution vanishes exactly.
We may also analytically compute the entropy of a strip (analogous to our slabs in higher dimensions). We take the strip to be thin compared to the adiabatic scale (L b/ b ). Solving the equations of motion gives The numerically derived surfaces agree with the above to one part in 10 14 to zeroth order and one part in 10 7 to first order. Computing the entanglement entropy gives Comparing this result to our d = 2 numerics shows discrepancies only at the level of one part in 10 4 level for the coefficient of b 0 2 and one part in 10 15 for the coefficient of b 0 .
With these expressions for the area, we can compute the location of the phase transition between vanishing and non-vanishing mutual information to second order. To zeroth order, for the half space entangling surface we haveÂ (0) = 0 so from (A.8)Â

JHEP09(2016)058
As a final check on our d = 2 results we can solve for the diffeomorphism taking g µν + 2 g (2) µν . Working near x = 0, we find that the correct diffeomorphism beomes which indeed takes the entangling surfaces of global AdS 3 to (A.7) as desired. One may also check that (A.12) maps the phase transition for b(x) = constant (given by (A.9)) to the value specified by (A.11).

B Estimation of numerical uncertainty
We have used two distinct methods to estimate the numerical uncertainty of our results. First, for the majority of the tables we merely make a rough estimate by computing a particular coefficient using a variety numerical parameters. We then take the approximate error to be given by the standard deviation of this set. For example for the L b 0 scaling of table 1, we compare values calculated using 100 and 150 lattice points and for fitting intervals L/b 0 ∈ [40, 50] and L/b 0 ∈ [50, 60]. The estimated error is the standard deviation of this four point data set. The value displayed in the table is the mean.
However, as noted in the main text, the values tabulated in table 4 are rather less uniform than one might expect. As a result, we now take extra care to analyze the numerical results reported there. After investigating the possible sources of error by varying the precision of different parts of the computation, we find the dominant error (by far) to come from using a finite number N of lattice points in the interval [−L, L]. We now study how our results change with N .
We first compute L crit using N = [50, 300] lattice points in steps of 10. Next, we approximate the function dL crit dN by choosing an appropriate p so that the data  Table 8. We display the estimated values ofD for each of the coefficients b 0 2 and b 0 b 0 and p for 2 ≤ d ≤ 8.