Entanglement entropy in cubic gravitational theories

We derive the holographic entanglement entropy functional for a generic gravitational theory whose action contains terms up to cubic order in the Riemann tensor, and in any dimension. This is the simplest case for which the so-called splitting problem manifests itself, and we explicitly show that the two common splittings present in the literature - minimal and non-minimal - produce different functionals. We apply our results to the particular examples of a boundary disk and a boundary strip in a state dual to 4-dimensional Poincar\'e AdS in Einsteinian Cubic Gravity, obtaining the bulk entanglement surface for both functionals and finding that causal wedge inclusion is respected for both splittings and a wide range of values of the cubic coupling.


Introduction
Classical stringy corrections lead to an effective higher-derivative theory of gravity. In such a theory, if the higher-derivative operators are suppressed by powers of P , we are guaranteed that the theory is well behaved. When the higher-derivative operators are unsuppressed we have to analyze each theory individually; general statements regarding the well-posedness of the theory are, alas, hard to come by. Even a fundamental property as causality has to be re-examined. The generic existence of superluminal modes implies that causality and hyperbolicity of the equations of motion are not guaranteed, and the analysis has to be carried out for each specific theory -or class of theories [1,2]. Despite, or maybe because of, all these characteristics, higher-derivative theories are interesting in more than one way. In holography they provide a testing ground where to understand more deeply how holography works in theories whose duals are more generic CFTs (different central charges, non-supersymmetric, etc). Crucial holographic constructs like the holographic entanglement entropy are modified in the case of higher derivative theories [3][4][5], and certain subtleties arise [6]. We will discuss this in detail in the next sections. The generic presence of superluminal modes in higher-derivative theories implies that previous causal constructs based on null rays have to be reexamined [7]. To add to the multitude of ways in which higher-derivative theories differ from Einstein gravity in the holographic context, there is also the question if a given higher-derivative theory is UV complete and, thus, expected to have a sensible field theory dual, or not. Holographically, relationships like causal wedge inclusion [8] and entanglement wedge nesting [9] are necessary conditions for a background to have a field theory dual and, thus, can be used to rule out certain higher-derivative theories from having unitary relativistic QFT duals. In [7] the authors showed how causal wedge inclusion can be used to arrive to the same conclusion as [10].
In the vast landscape of higher-derivative gravities, Lovelock theories are among the most studied; they have the advantage that they yield second order equations of motion. Among them, the quadratic theory, Gauss-Bonnet, has served as a prototype for many phenomena not present in Einstein gravity. From the violation of the η/s bound [11,12], to recent work related to the information paradox [13], Gauss-Bonnet theory has taught us important lessons for holography. Cubic theories of gravity have been studied in the general relativity community [14,15]. Some of their holographic properties have been explored [16,17], but the explicit form of the entanglement entropy functional, a fundamental quantity in holography, was not known. 1 In this paper we advance the understanding of cubic theories in a holographic context by deriving the holographic entanglement functional for a generic cubic gravity theory. This functional can be applied to cubic Lovelock, quasi-topological gravity and Einsteinian cubic gravity theories. Our result is applicable to general cubic gravity theories in any dimension, with the only restriction that the action does not involve derivatives of the 1 There has been previous work on some particular cubic theories, such as quasi-topological gravity [18,19].
In any case, the general cubic functional was not known, and the "splitting problem" -which we discuss in section 2 -was overlooked. curvature tensors. Obtaining the functional presents subtleties absent in quadratic theories. In [6,[20][21][22], the authors showed that, in general, there is an ambiguity in the calculation of the entanglement entropy functional. This ambiguity, known as the "splitting problem", is related to the regularization of the action near the conical singularity that appears in the Lewkowycz-Maldacena prescription [23]. We investigate this issue in detail and present the functional using two different splittings that we refer to as "minimal" and "non-minimal". The "non-minimal" prescription is known to be correct at the perturbative level. At finite coupling, determining the correct splitting is an open question. However, we illustrate our result calculating the HEE surface in a theory that was built to avoid causality problems at the perturbative level, Einsteinian Cubic Gravity (ECG).
The structure of the paper is as follows. In section 2 we review the general framework for calculating the holographic entanglement entropy in higher-derivative theories. We pay particular attention to the spliting problem and to the two different proposals that exist in the literature to solve it. Section 3 contains our main result: we derive the entanglement entropy functional for a general cubic gravity theory. We explore the result obtained using the two different splitting prescriptions. We point out that quadratic theories are insensitive to the differences between them, as also are Lovelock theories, for which the Jacobson-Myers functional is valid [24]. However, for a generic cubic theory, the minimal and non-minimal prescriptions lead to different answers. As an example of our results, in section 4 we work out in detail the entanglement functional for a particular cubic gravity theory, Einstenian Cubic Gravity (ECG) [15], and present some numerical results regarding the minimal surface they produce. Finally, in section 5 we summarize our results, its implications, and point out open directions.

Holographic entanglement entropy in higher-derivative gravity
In a holographic CFT dual to Einstein gravity, the entanglement entropy of a boundary region A is given by the area of an associated codimension-2 surface [25,26]: The surface γ A , also known as the Ryu-Takayangi or RT surface, is defined as the bulk codimension-2 surface which has the minimal area among all those homologous to the region A in the boundary (it has to end in ∂A if this is not empty). If we assume holography holds, the work of Lewkowycz and Maldacena [23] constitutes a proof that (2.1) indeed gives the entanglement entropy. We will briefly review their argument here in order to set the stage for our future discussion concerning field theories dual to higher-derivative gravities. The computation starts by considering the usual replica trick in the boundary field theory. The entanglement entropy S can be computed as the limit n → 1 of the Rényi entropies: where ρ A is the reduced density matrix of the subsystem associated with region A, and Z n is the partition function of the field theory in the n-fold cover. This is a manifold consisting of n copies of the original one, glued cyclically at the spatial region A. Z 1 is thus just the original partition function. We assume always that an analytic continuation to Euclidean signature has been performed. Notice also that Rényi entropies are defined for n ∈ N , therefore an analytic continuation in n is also assumed before taking the limit n → 1. So far, all this discussion has been restricted to the field theory, but if this is holographically dual to a gravitational one, it should be possible to find a bulk solution B n dual to the n-fold cover. Then, log Z n = −I[B n ], where I[B n ] is the on-shell gravitational Euclidean action of this dual geometry. Naturally, log Z 1 = −I[B 1 ], B 1 being just the original bulk dual. Now, the n-fold cover boundary manifold has a Z n symmetry, due to the fact that we can do permutations on the n copies of the original manifold. If we assume that this replica symmetry is respected in the bulk, we can consider the manifoldB n = B n /Z n , which has to be regular everywhere except in the codimension-2 submanifold C n consisting of fixed points of the Z n . Notice that, in the boundary, ∂A are precisely the fixed points of Z n . Also, since B n is a regular bulk solution, the orbifoldB n has a conical defect of opening angle 2π/n at C n . Due to the replica symmetry, we can write: where in I[B n ] we exclude contributions coming from the conical singularity (this is because in the left-hand side of the previous equation there are no such contributions, the geometry is regular). 2 After doing a suitable analytic continuation of thisB n to non-integer n, we can finally write: Once this expression is obtained, the computation of the entanglement entropy of the region A has been reduced to a problem in classical gravity, which can be solved in two steps: 1. The geometry with n = 1,B 1 , is a regular solution of the equations of motion. In (2.4) we seem to be doing a first order variation away from this solution, so we could naively expect that expression to vanish. This is not so because, when varying, we are changing the opening angle at C 1 = lim n→1 C n , which as mentioned should be excluded from the action integral and that procedure introduces a boundary where conditions are changing if we vary n. This localizes the computation of the entanglement entropy in C 1 , and in fact in Einstein gravity it is possible to prove from the form of the action (see [23]) that S is computed as shown in (2.1). γ A should be interpreted at this point as C 1 , where we have not proven its minimal property yet.
2. The remaining question is how we determine C 1 . Formally, it is defined by looking for C n in the analytically continued spacetimeB n , and then taking the limit n → 1. Adopting adapted coordinates at the conical singularity (see appendix of [23]), it is possible to show that the equations of motion derived from Einstein gravity forB n impose, in the limit n → 1, the minimal area condition: where K a are the traces of the extrinsic curvatures along the transverse directions to C 1 , and a is an index which runs in these two directions (this notation will be clarified in the following section). This shows that C 1 is a minimal area surface, which can then be calculated by minimization of the entanglement entropy functional (2.1). With this condition, C 1 can be characterized as the previously defined surface γ A which is homologous to the boundary region A.

The entropy functional and the splitting problem
The previous program can be carried out, with an increasing level of technical difficulty, when the gravitational theory contains higher-derivative corrections to the Einstein-Hilbert action.
Following the two steps we have just described, [4] and [5] first obtained the expression for the functional computing entanglement entropy in the presence of higher-derivative terms. Using (2.4), and considering a Lagrangian containing arbitrary contractions of the Riemann tensor (but not its derivatives), one obtains: where L E = L E (R µνρσ ) is the Euclidean version of the Lagrangian. Let us explain the notation in this expression, which employs the conventions of [4]: • The full manifold has dimension D, and we denote generic coordinates in it by x µ . Indices for this D-dimensional manifold are µ, ν, ρ, σ, . . . The surface on which the previous functional is evaluated, C 1 , is (D − 2)-dimensional. Coordinates in it are y i , and indices will be labelled i, j, k, l, . . . We assume an embedding x µ = x µ (y i ), so that we can define tangent vectors (m i ) µ ≡ ∂ i x µ , and then take two extra orthonormal vectors n a to complete the basis, G µν (n a ) µ (n b ) ν = δ ab . Indices a, b, c, d, . . . will be used for these two directions.
• The functional (2.6) is defined using a particular set of adapted coordinates for C 1 (see [4]), where tangent coordinates are x i (y) = y i , and we introduce two extra complex coordinates z,z such that the metric factorizes: with the remaining components vanishing at C 1 .
• K a ij is the extrinsic curvature of C 1 along the direction n a : There remains to explain the sum over A in the last term of (2.6), usually called the anomaly term. Its origin is subtle, coming from potentially logarithmically divergent terms in the action at C 1 when taking the limit n → 1. For those terms, a careful analysis of the limit has to be performed, and it becomes essential to understand the analytic continuation of the geometry and the regularization of the conical singularity at C 1 . 3 This regularization has to be done guaranteeing that the equations of motion of the theory are satisfied. This was at first overlooked in [4], where a minimal regularization was employed. It is nevertheless useful to quote the result obtained, since it will allow us to see more clearly the differences introduced when other regularizations are used. For the computation of the entropy functional, all the study of the behaviour of the action around C 1 boils down to the following algorithmic procedure. In a general theory, the second derivative of the Lagrangian in the last term of (2.6) will be a polynomial in curvature tensors. In this polynomial we expand every curvature tensor using the following expressions: (2.10) where indices α, β, γ, . . . denote values z orz, r ikjl is the lower-dimensional Riemann tensor, andR αβij ,R αiβj and Q αβij are defined in [4] but their particular form will not be needed here. Once the expansion is done, we label each of the individual terms with A. Considering any of them, we associate a value q A to it equal to the number of factors of Q zzij and Qzz ij plus one half the number of factors of K αij , R αβγi and R αijk . Once this is done, we divide by q A + 1, as indicated in (2.6). This expansion is what we denote by the sum over A, and when it is completed one can rewrite, if desired, everything again in terms of the original curvature tensors using (2.10). As mentioned, this prescription (which we will call minimal prescription) does not take into account the fact that the regularized metric at the conical singularity has to satisfy the equations of motion when computing (2.4). The existence of several ways to regularize the metric (with relevant consequences for the entropy functional) has been called the splitting problem in the literature, and discussions around this issue can be found in [20][21][22].
For a general higher-derivative theory, it can be quite complicated to impose the onshell condition for the regularization, but one can make a first approach to the problem by studying the constraint imposed by Einstein gravity. This produces a functional which is at least perturbatively correct, since the leading order area term is independent of the splitting chosen. This is done in great detail in [22] (and in [6], although the final result is expressed assuming the surface has K a = 0, so that it has extremal area), and it is shown that the final effect for the entropy functional is a different prescription for the A expansion (we will call this non-minimal prescription). After expanding all the curvature tensors according to (2.10), we have to rewrite R zzzz and Q zzij in terms of two new objects: and associate q A = 1/2 to K αij , R αβγi , and R αijk ; and q A = 1 to Q zzij and Qzz ij . Then we must divide by q A + 1 as indicated in the general expression (2.6), and finally, if desired, undo the expansions, writing everything in terms of curvature tensors again.
We have presented two ways to define the functional computing entanglement entropy in the presence of higher-derivative corrections to the gravitational action. Let us emphasize that for quadratic theories the functionals obtained using the minimal and non-minimal prescriptions are the same. As shown in [4], this is also the case for Lovelock theories, where the Jacobson-Myers functional is known to be correct. This is a special property of Lovelock theories, for which the entanglement entropy functional depends only on the intrinsic geometry of the surface, r ijkl . However, we will show that for general cubic gravities the functionals obtained are different.
Until now, we have only completed the first step of the general strategy outlined in the previous section, i.e., we have shown how to obtain the entanglement functional starting from the action of the theory. We still have to find the position of the surface C 1 , in other words, where to evaluate the functional. In principle, the equations of motion of the theory should determine the location of the surface, in much the same way they determine it to be a minimal area surface in General Relativity. The idea is to explicitly evaluate the equations of motion in the conically singular metric for generic n to linear order in n − 1. One does not expect divergences in the stress-energy tensor, but these generically appear in the metric part of the equations of motion, so cancelling them imposes conditions which determine the position of the surface in the limit n → 1 (in particular, GR equations of motion impose the well-known condition K a = 0). Further details can be found in [21,23]. In practice, this procedure has the drawback of requiring to deal with the equations of motion of higher-derivative theories, which can be extremely complicated. Fortunately, in [27] it is shown that the same procedure one employs in Einstein gravity is also valid in general: after computing the correct functional, one can minimize it to obtain the surface C 1 in which it is going to be evaluated. 4,5

Entanglement entropy functional in cubic gravity
This section contains our main result, we will derive the entanglement functional for a generic cubic gravity theory that does not involve explicit derivatives of the curvature tensor. We will use the method diuscussed in section 2 to compute the different contributions to the holographic entanglement entropy functional in cubic gravities.
As a warm-up exercise, let us revisit the known result for quadratic theories. As explained in section 2, there are different prescriptions to calculate the entanglement entropy functional in higher-derivative theories. In quadratic gravity theories, any prescription leads to the same result for the holographic entropy functional [29]. The reason for this is that for quadratic theories the expansion in A is trivial: after two derivatives of the Lagrangian we obtain something which does not contain curvature tensors, so essentially q A = 0 always. This guarantees that, in the case of quadratic gravities, the result obtained using either the minimal or non-minimal splitting is the same, where K a ≡ K aij g ij . This is not the case for cubic gravities. Consider the following generic cubic Lagrangian: 6 Note that the second derivative of the Lagrangian (3.3) is linear in the curvature tensor. Therefore, unlike in quadratic gravity, the two different splittings discussed in section 2 lead to different entanglement functionals. The details of these highly technical calculations are presented in Appendix A. The final result for the entropy functional obtained following the minimal splitting prescription is: Using the non-minimal prescription we obtain: where the first two terms are given by (3.5) and (3.6), and the last one is: Equations (3.4) and (3.8) are two of the main results of this paper. Let us pause and take stock of these results. First, it is natural to inquire whether we can identify where the difference in the two functionals comes from. As the separation of the functional in three different parts makes manifest, we observe that terms which are proportional to the square of background curvature tensors are equal in both prescriptions and given by (3.5). This is because they come from the first term in the general functional (2.6) (the Wald term), which is independent of the splitting. The same thing happens for terms linear in background curvature tensors (3.6): although these come from the A expansion, they are such that after the rewriting and counting procedures the result is the same for both prescriptions. At the end, differences arise only in K 4 terms. Second, as previously mentioned, determining the correct splitting for a generic cubic theory with finite coupling is an open question. It could be a splitting as yet unknown, different from the previous ones. However, after having obtained the entanglement entropy functional using the minimal and non-minimal prescriptions one can turn the question around and ask: if we use these functionals for a cubic theory with a finite coupling, will either of them produce unwanted, unphysical, behaviour? In the next section we set out to do just that, we investigate a fundamental property known as causal wedge inclusion. This property states that, in a spacetime that has a CFT dual, the causal wedge is completely contained in the entanglement wedge, C(A) ⊆ E(A). Casual wedge inclusion can be used as a criterion to constrain the space of theories with CFT duals. Studying the causal structure of a generic higher-derivative theory is usually a thorny issue because in these theories gravity can travel slower or faster than light. The causal structure is determined by characteristic hypersurfaces that are generically non-null. A thorough study of this issue was carried out for Gauss-Bonnet theory in [1,2], but the understanding of causality and hyperbolicity properties in cubic theories is an open problem. Therefore, when investigating causal wedge inclusion we will restrict ourselves to a case where we are guaranteed that the causal structure is given by null rays.

A concrete example: Einsteinian Cubic Gravity (ECG)
In the previous section we derived the entanglement entropy functional for a generic cubic gravity theory in D dimensions with the only restriction that the action does not involve derivatives of the curvature tensors. Many such theories exist in the literature [14,15]. In this section we will focus on one of these theories to carry out a concrete calculation in full detail. We will compute the entanglement entropy functional of both a disk and a strip in the state dual to vacuum AdS in 4-dimensional Einsteinian Cubic Gravity (ECG), and explore properties of the entanglement surface.

Preliminaries
Before proceeding with the calculation of the entanglement functional in ECG, let us review some aspects of this theory. We refer the reader to the original works [17,31,32] for further details.

Einsteinian Cubic Gravity
Einsteinian Cubic Gravity (ECG) [15] is the unique theory containing corrections up to cubic order in the Riemann tensor such that: 1. It has the same spectrum than Einstein gravity when linearized on a maximally symmetric background, that is, it propagates a single transverse and massless graviton. 2. The coefficients appearing in the Lagrangian are dimension-independent. 3. The cubic correction is neither trivial nor topological in four dimensions.
The action of ECG is where A similar construction exists for higher orders in the derivative expansion [33]. Black holes solutions in this theory exist in the literature and some of their properties have been studied in [31,32]. If we demand stable AdS 4 vacua, the coupling µ is constrained to a range of values determined by the equation where f ∞ relates the curvature scale of the AdS background L and the action length scale L, To have AdS 4 vacua, the roots of (4.3) should be f ∞ > 0, and thus µ ≤ 4 27 . The root f ∞ which produces a stable AdS vacuum is For µ < 0, there is a single stable (i.e., positive effective Newton constant) vacuum, but it is not possible to have black hole solutions. For 0 ≤ µ < 4 27 , there is one stable vacuum connected to Einstein gravity in the limit µ → 0, and it is possible to have black hole solutions. This root satisfies f 2 ∞ < 1 3µ . 7 Finally, holographic studies show that positivity of energy fluxes at null infinity in the CFT imposes a more stringent constraint in the coupling: −0.00322 ≤ µ ≤ 0.00312. In [17], the authors argued that the holographic dual of ECG is a non-supersymmetric CFT in three dimensions, 8 and obtained some entries of the holographic dictionary.  [8,34]. Subregion duality and AdS/CFT imply a constraint between these two holographic constructs: the causal wedge is completely contained in the entanglement wedge, C(A) ⊆ E(A). This relation is known as causal wedge inclusion. Backgrounds that do not satisfy causal wedge inclusion are not viable as holographic duals of a boundary field theory. 7 A second, positive root of the characteristic equation with f 2 ∞ > 1 3µ gives rise to an unstable vacuum, in which the effective Newton constant becomes negative. The case µ = 4 27 , f 2 ∞ = 1 3µ = 9 4 is a critical limit of the theory, in which the roots merge and the effective Newton constant blows up. We will not consider it here. 8 One of the parameters characterizing the three-point function of the stress tensor, t4, is shown to be different from zero, contrary to what happens in a supersymmetric CFT.

Entanglement entropy of a disk in 4-dimensional ECG
Consider Euclidean AdS 4 written in boundary polar coordinates: Take the boundary region to be defined as r ≤ R, φ ∈ [0, 2π) at z = 0 and some fixed τ . This is a disk which will produce an entanglement surface penetrating into the bulk which we choose to parametrize as r = ξ, z = Z(ξ) (plus the angular, symmetric direction). In Appendix B, the different geometric quantities relevant for the computation of the entanglement entropy functional associated with this surface can be found. We will present in a moment the explicit form of the functionals following both the minimal and the non-minimal prescriptions, but notice a relevant fact that can already be anticipated at this point. For both funcionals (3.4) and (3.8) there is a term (S R 2 ) which involves quadratic curvature contractions (proportional to 1/L 4 when evaluated on pure AdS, independently of the entanglement surface) and an extra piece (S K 2 R + S K 4 ) which depends on the prescription and which is at least quadratic in extrinsic curvatures. As shown at the end of Appendix B, the RT or minimal area surface with correct boundary conditions is just an spherical shell Z(ξ) = R 2 − ξ 2 , which happens to have vanishing extrinsic curvature. This surface will therefore also minimize the complete functional (whether it is (3.4), (3.8), or in fact any other prescription), because the S R 2 piece in pure AdS will be just a constant times the area of the surface, while the extra parts will produce terms in the Euler-Lagrange equations for extremization proportional to the extrinsic curvature, which therefore vanish for the spherical shell. Just as a reassuring check, we can explicitly write the functionals obtained with each of the prescriptions. For the minimal one: For the non-minimal prescription: While not immediate to discern from these complicated expressions for the functionals, one can compute their Euler-Lagrange equations and verify that Z(ξ) = R 2 − ξ 2 is indeed a solution, as the general argument presented before shows.
Since the entanglement surface coincides with the minimal area one (i.e., the one we would have obtained without higher curvature corrections), and taking into account that the causal structure around pure AdS for ECG is also the same than for Einstein gravity, we conclude that causal wedge inclusion is respected in this case just like it is respected in pure Einstein gravity. Incidentally, this means that it is marginally respected, which is to say that both the entanglement and the causal wedge coincide, as can be easily checked by verifyng that the projection of the causal wedge into the fixed τ slice coincides with the entanglement surface, z 2 + r 2 = R 2 . The simple disk geometry in the boundary we have considered is therefore not able to constrain ECG in any way, since it respects the conditions imposed by subregion duality in AdS/CFT. Let us consider a different geometry to see whether we can extract any useful information.

Entanglement entropy example: a strip in 4-dimensional ECG
Similarly to what we did in the previous section, consider Euclidean AdS 4 as the bulk metric, but written now in boundary Cartesian coordinates: and take the boundary strip to be x ∈ (− /2, /2), y ∈ (−∞, ∞) at z = 0 and some fixed τ . We parametrize the entanglement surface with (x, y) as z = Z(x), since y is a symmetry direction. We refer to Appendix C for the computation of the relevant geometric quantities of this entanglement surface. Taking those results into account, we can write the functional for a strip in ECG (4.1) following both of the splitting prescriptions. With the minimal one we obtain, starting from (3.4): (4.12) We can start from (3.8) instead to obtain the functional following the non-minimal prescription: where now S non−min ≡ − 6 + 24(Z ) 8 (4.14) A couple of comments are in order. First of all, recall that in order to have AdS with curvature radius L 2 = L 2 /f ∞ as a background in ECG, f ∞ must satisfy: and we are taking the solution of this equation which has positive effective Newton constant and connects with the GR solution f ∞ = 1 in the limit µ → 0, given by (4.4). Second, both functionals have an obvious IR divergence, since the y integral is infinite. An IRregulator must be therefore included, cutting the strip at some fixed length. Notice also another suprising feature of the previous functionals. As already mentioned, the entanglement entropy of the strip should be given by the functional obtained employing the correct splitting prescription for ECG (which might be none of the previous ones), and evaluated at the surface determined by the Z(x) which extremizes that functional. If the correct splitting happens to be one of the previous ones, then we have to obtain the Euler-Lagrange equation for Z(x) starting from the corresponding functional. In both cases, we have up to second derivatives of Z, so the corresponding differential equation for Z will be fourth order. This contrasts with the second-order character of the equations of motion for the perturbations around a maximally symmetric background in ECG, which is one of its defining properties. In the following section we disregard this fact and assume each of the functionals to be correct, as a test to see what would happen. We numerically solve for the surface profile Z(x) which extremizes the corresponding functional, 9 and we plot the result. This will serve as a probe to understand how the higher-derivative corrections to the entanglement entropy change the bulk surface in which the functional is to be evaluated.

Numerical results
In this section we present some of the curves for the numerical solutions obtained by minimizing the functionals (4.11) and (4.13). The boundary region is a strip of width = 3 in the x direction and infinite in the y direction. We solved the fourth-order equations (D.11) and (D.12) to obtain the corresponding entanglement wedge. All the details regarding the numerical procedure, along with more cases of µ values for which the solution was obtained, are presented in Appendix D.2. The resulting plots, in which we also include the causal wedge, can be found in Figures 1 and 2.   These plots show that the causal wedge is safely included in the entanglement wedge for both prescriptions, and for all values of the coupling within the range tested (in Appendix D. 2 we present examples which prove this to be the case for, at least, −10 4 ≤ µ ≤ +0.01). Thus, just like for the boundary disk, causal wedge inclusion applied to the boundary strip in Poincaré AdS does not constrain the allowed values of the coupling in ECG. Note that causal wedge inclusion in a Gauss-Bonnet black hole background does constrain the allowed values of the coupling [7]. A crucial ingredient to arrive at the results of [7] was the modified causal structure of the higher curvature theory due to the presence of superluminal modes. Here we are working in AdS and the causal structure is still governed by null rays. In addition, the disk has the same entanglement surface in both Einstein gravity and ECG, and therefore although causal wedge inclusion is marginally respected, it is so for all values of the coupling. For the strip geometry, in the case of Einstein gravity (µ = 0), causal wedge inclusion is safely respected, as the previous plots show. Thus, to violate causal wedge inclusion, it would be necessary a big modification of the surface due to turning on the coupling µ. This does not happen.
Clearly, to explore the effects of causal wedge inclusion in ECG it would be better to study a black hole background, just like the authors did in [7] in Gauss-Bonnet gravity. However, in order to construct the causal wedge in a higher curvature theory, a complete investigation of the superluminal modes that determine the causal structure is required. A different possible avenue to harness the power of causal wedge inclusion is to investigate perturbative modifications of the disk region and solve the Euler-Lagrange equations derived from (4.6) or (4.8) with the perturbed boundary condition, but still in pure AdS, similarly to what [36] did for Einstein gravity. Since causal wedge inclusion is marginally respected there, it could happen that a perturbative violation occurs for some values of the coupling, which would therefore have to be excluded. The complicated form of the Euler-Lagrange equations makes this an extremely involved problem at the technical level, which we therefore leave as an open question.

Conclusions and future directions
In the present paper we have obtained the entanglement entropy functional for a generic cubic gravitational theory not involving derivatives of the curvature tensor in the action. This constitutes our main result. We have done this using two different prescriptions: the minimal one, introduced in [4], and the non-minimal one, presented in [6,22] and which is known to be perturbatively valid. The existence of these two alternative functionals is an explicit manifestation of the splitting problem one is forced to face when deriving the entanglement entropy functional for cubic or higher order gravitational theories. Despite knowing that for a particular theory there must be a single, correct functional, we have also performed some consistency checks for both of them in a particular cubic theory. In particular, we investigated whether the functionals obtained for Einsteinian Cubic Gravity produce via minimization entanglement surfaces which satisfy the causal wedge inclusion property for both a boundary disk and a boundary strip in Poincaré AdS. We find that, for all the values of the couplings studied, causal wedge inclusion is respected.
Our work makes it possible to investigate several questions that will advance our understanding of the role cubic theories play in a holographic context.

Bit threads
In [37] the authors put forward an alternative formulation of the holographic entanglement entropy that does not rely on minimizing an area functional but invokes a divergenceless vector field, dubbed bit threads. Many aspects of bit threads have been studied [38][39][40][41][42][43], but a formulation of bit threads for higher-derivative theories was missing. Recently this gap was closed in [44], where the authors derived a bit thread formulation for a general higher-derivative theory. Now that we have the entanglement functional for cubic theories in the minimal and non-minimal splitting, it would be interesting to investigate if using bit threads we can understand better the splitting problem.

Dynamics
Holographic entanglement entropy in dynamical backgrounds has been widely studied in Einstein gravity, with and without charge. It led to interesting insights regarding the thermalization time [45][46][47]. Similar work has been carried out for Gauss-Bonnet backgrounds [48][49][50]. If a Vaidya type of solution can be written down for black holes in cubic gravity, then dynamical studies in these theories along the lines suggested would be quite interesting.

More general boundary regions
We have learned valuable lessons studying HEE for different boundary regions: as an example, the divergence structure of the holographic entanglement entropy on regions with corners uncovered cut-off independent coefficients. These coefficients were shown to be universal and to encode important field theory data in Einstein and Gauss-Bonnet theories [51,52]. It would be interesting to study regions with corners in cubic theories.
Derive the correct splitting for finite coupling Clearly, an important open question is to formally find the correct splitting for a general cubic theory with finite coupling. This highly technical problem can be approached following [6,22].
Causal structure of cubic theories In section 4, as an example of the functional we derived, we calculated the entanglement surfaces of a disk and a strip regions in an AdS background solution of ECG. More interesting phenomena regarding the causal and entanglement wedges can be expected if we were to consider black holes in any cubic gravity theory. However, black hole backgrounds in higher-derivative theories will generically have superluminal modes, and their causal structure is no longer determined by null rays. A complete study of the causal structure and hyperbolicity of equations of motion of cubic theories similar to what was done for Gauss Bonnet [1,2,53] is of interest not only for the GR communities but also in a holographic context.

Perturbative calculations in field theory
The non-minimal functional is known to be the correct one for any perturbative cubic theory which does not include covariant derivatives of the curvature tensors. Furthermore, perturbatively, we know that the entanglement entropy functional can be evaluated in the Ryu-Takayanagi surface obtained in General Relativity. This is because the area term of the functional is minimal for the Ryu-Takayanagi surface, and therefore even if the surface changes at first order, it does not produce a change of that part of the functional. The first order variation of the Ryu-Takayanagi surface can also be neglected for the remaining part of the functional coming from the higher-derivative terms. This approach will be explored elsewhere [54].

Acknowledgments
We are deeply indebted to Joan Camps for insightful explanations and correspondence. It is a pleasure to thank also Pablo Bueno, José Edelstein and Julio Oliva for useful discussions. EC and RC are supported by the National Science Foundation (NSF) Grant No. PHY1820712. RC is also supported by NSF Grant No. PHY-1620610. The work of AVL is supported by the Spanish MECD fellowship FPU16/06675, and by MINECO FPA2017-84436-P, Xunta de Galicia ED431C 2017/07, Xunta de Galicia (Centro singular de investigación de Galicia accreditation 2019-2022) and the European Union (European Regional Development Fund -ERDF), "María de Maeztu" Units of Excellence MDM-2016-0692, and the Spanish Research State Agency. AVL is also pleased to thank the University of Texas at Austin, where part of this work was done, for their warm hospitality.

A.1 A first, detailed example
In order to understand better how the entanglement entropy functional is obtained for a generic cubic theory, let us present an example in detail. Consider the following (Euclidean) Lagrangian: where λ is a constant. We will compute separately each of the terms in the general form of the functional (2.6), and for the second one we will do it following the two prescriptions presented in section 2. First of all, for the Wald term, we need the following derivative: The simplest way to do this is by first writing the expression as a contraction in the complex coordinates z andz, using the fact that the only non-vanishing components of the metric are G zz = Gz z = 1/2: where the antisymmetry of the Riemann tensor in the two pairs of indices has been taken into account. Now we use the fact that a contraction in an α-index can be substituted for a contraction in orthonormal directions, since tangent directions have no z orz components: where T is any tensor, possibly containing extra indices. The Wald term is then: This is written in a form which can be evaluated in any set of coordinates just constructing the two orthonormal vectors n a to the surface. Consider now the anomaly term. The first step is to obtain the second derivative of the Lagrangian: Now we have to expand R following the two prescriptions. For this we have to write it in terms of Riemann tensor components and expand them following (2.10): Now, in the minimal prescription the last two terms have q A = 1, while all the rest have q A = 0 (notice that G αβ Q αβij involves only Q zzij ). Therefore, when doing the A expansion we have to multiply them by 1/2. Doing that and then rewriting back everything in terms of the Ricci scalar, the minimal prescription gives: For the non-minimal prescription we need to rewrite R zzzz and Q zzij in terms of the primed versions (2.11). This is done as follows: The Ricci scalar is now written in terms of the basic objects for this prescription as: so that there are no terms with non-vanishing q A . We can then immediately rewrite everything in terms of the Ricci scalar, obtaining: This completes the A expansion for the non-minimal prescription. Notice that, in (A.7), the metric tensors are not affected by the expansion, and we can contract them with the extrinsic curvatures appearing in the general formula (2.6) as follows: All contractions in α indices can now be traded for contractions in a indices as explained when discussing the Wald term. We conclude this little example by collecting all contributions, which would produce the following entanglement entropy functionals for the Lagrangian (A.1): The following sections contain the contributions to both the Wald and anomaly terms of the functional using both prescriptions. We include all numerical factors except for the global 2π appearing in (2.6).

A.2 Minimal prescription
A.3 Non-minimal prescription r ≤ R, so the surface anchored at its boundary and going into the bulk is parametrized as: with ψ p coordinates in the (D − 3) unit sphere, ξ ∈ (0, R), and Z(R) → 0. Basis vectors tangent to the surface are then: This induces a metric on the surface of the form: We then choose our two normalized vectors orthogonal to the surface to be: This produces the following extrinsic curvature components: where γ pq is the metric on the unit round sphere and the rest of the components of the extrinsic curvature vanish. Transforming the last two indices to coordinate ones: The trace of the second extrinsic curvature is then (clearly K 1 = 0): Notice that the RT surface satisfying the prescribed boundary conditions (which is the one solving K 2 = 0) is given by Z(ξ) = R 2 − ξ 2 . In particular, this implies ZZ = −ξ and ZZ = −(1 + Z 2 ), showing that not only K 2 = 0, but the whole tensor satisfies K 2 µν = 0. This fact is relevant when discussing minimal surfaces for boundary disks in the main body of the paper.

C Geometry of the entanglement surface for a boundary strip
The aim of this short appendix is to collect the results needed to evaluate the tensors appearing in the entanglement entropy functional for a boundary strip in Poincaré AdS. Consider then the bulk metric: where the length scale L is determined by imposing this to be a solution of the equations of motion for the theory at hand. We have separated the spatial coordinates in the boundary into x and y p (with p = 1, 2 . . . , D − 3) because we will consider in this spacetime a surface anchored to a boundary strip finite in extent in the x-direction, parametrized as: Basis vectors tangent to the surface are then: This induces a metric on the surface of the form: We then choose our two normalized vectors orthogonal to the surface to be: This produces the following extrinsic curvature components: with the rest of them vanishing. Transforming the last two indices to coordinate ones: while K 2 pq = K 2 ψ p ψ q . The trace of the second extrinsic curvature is then (clearly K 1 = 0): As a final comment, we are employing a generic parametrization, but the results in the main text are presented with z = Z(x). For that case, we can just set X = 1. For numerical computations we also employed x = X(z), in which case we set Z = z and Z = 1.

D Details of the HEE in ECG calculation
This appendix contains the calculational details of the example presented in section 4. First we will obtain the equations to solve, and then present details of the numerical integration.

D.1 Minimizing the entropy functional for a strip in ECG
Consider the boundary strip in a state dual to the 4-dimensional vacuum AdS in ECG, with bulk metric given by (4.10). For this setup, the functionals (3.4) and (3.8) have been already obtained in the minimal and non-minimal prescriptions for the z = Z(x) parametrization, eqs. (4.11)-(4.14). Similar expressions are needed for the x = X(z) parametrization, since it will be used in the first two parts of the numerical procedure. With the minimal splitting regularization we obtain Minimizing this functional we obtain a fourth order differential equation for X(z) to solve, If we use the non-minimal splitting instead, the functional obtained is Minimizing this functional, the equation to solve is 8 2X + 2(X ) 3 − zX 1 + (X ) 2 5 + 3µf 2 ∞ −184(X ) 9 − 72(X ) 11 + 6z(X ) 10 X − 12(X ) 13 + 4(X ) 7 −64 + 5z 3 X X (3) + zX 8 + 13z 2 (X ) 2 + 36z 3 X X (3) + It will prove convenient for the numerical calculation to work also with the entanglement functional in terms of Z(x) instead of X(z). Recall that, with the minimal prescription, the functional reads: (D.8) And following the non-minimal prescription: where now S non−min ≡ − 6 + 12(Z ) 8 + 4(Z ) 10 (D.10) Then, the equations to solve are if we use the minimal splitting, and when using the non-minimal splitting.

D.2 Numerics
We consider an interval of width /2 = 1.5. The Z(x) parametrization is problematic if we want to start integrating from the boundary keeping the endpoints of the interval fixed. It is more convenient to work -at least initially-in terms of X(z). Our strategy will be to start with a series expansion for X(z) close to the boundary, at X(0) = /2, then numerically integrate equation (D.3) (or (D.6) depending on what splitting are we considering) and, finally, when the X(z) parametrization becomes problematic, switch to Z(x) and numerically integrate (D.11) (or (D.12)) until x = 0, where we reach the deepest point of the curve, Z(0) = z * .

D.2.1 Series expansion close to the boundary
To solve the relevant equation ((D.3) or (D.6)) starting from the boundary, we perform a series expansion of X(z) and the corresponding differential equation to order 23 for values close to z = 0. The zeroth-order term is determined by the boundary condition X(z = 0) = /2, whereas all the other coefficients depend on the value of the third-order coefficient in the expansion, related to X (0) as well as the value of µ. This series expansion is a good solution up to some small z = . The value of is chosen such that the numerical error is less than 10 −5 at any point 0 ≤ z ≤ , see Figure 3a for a representative case.

D.2.2 Numerical integration
For ≤ z ≤ z I , we integrated numerically the differential equation for X(z). The value z I was determined for each case as the largest value of z for which the errors remained below 10 −5 at any point in the interval of integration -see figure 3b. At the point (x I , z I ), we changed parametrization to z = Z(x) and we integrated numerically until we reached the z axis, thus completing the solution.
Recall that the value of the third derivative in the series expansion at z = 0 is so far an initial free parameter, so for a single value of µ we obtain a family of solutions characterized by different values of the third derivative. In figure 4, each curve corresponds to a different choice for the initial parameter. Note that this added complication arises because we are dealing with fourth order differential equations, unlike the second order ones we encounter in Einstein and Gauss-Bonnet gravity. All these solutions are good solutions of the differential equation. However, in a holographic context we expect that the curve will be smooth at x = 0. Therefore, among all the possible values of the third derivative we have to choose the one that produces a curve with Z (x = 0) = 0, and this curve will be the RT surface. In Table 1 we list the deepest point of the RT surface, z * , obtained for different values of the ECG coupling µ. The causal wedge in AdS is known to be a semicircle, thus, for a boundary region with /2 = 1.5 the deepest point of penetration of the causal wedge is z c = 1.5. If, for any value of µ considered, we find that z * < z c , this would be a clear violation of the causal wedge µ z * Minimal prescription z * Non-minimal prescription −10 4 2 inclusion. For all the µ values for which the solution was found there is no indication of such violation, see Table 1. We also show the plots for various cases of µ values for which the solution was obtained in Figures 5 and 6. Note that as we make µ more negative, the value of z * decreases for both prescriptions; even for µ = −10 4 , which is the most negative value for which the solution was obtained, we verified that z * > z c . : Causal surface (orange, dashed) and entanglement surfaces corresponding to ECG in the case of minimal prescription (blue), ECG in the case of non-minimal prescription (red), and Einstein gravity (green, dashed) [35] for the boundary strip length /2 = 1.5 and different values of µ outside the interval −0.00322 ≤ µ ≤ 0.00312. In both cases, we verify that zc < z * is satisfied for both prescriptions.