Inhomogeneous Jacobi equation and Holographic subregion complexity
Avirup Ghosh Rohit Mishra
hep-th gr-qc













    #1#1#1#1 g
    #1#1 ⟶#1#1#1#1#1
    ∂∇ tr Tr#1⟨#1⟩ L x y d1/2
    cosec sec
    avirup.ghosh@iitgn.ac.inIndian Institute of Technology, Gandhinagar, 382355, Gujarat , India.
    
    
    mrohit996@gmail.comSchool of Physical Sciences, Indian Association for the Cultivation of Science,
    2A & 2B Raja S.C. Mullick Road, Kolkata 700032, India
    We derive a general expression for obtaining Holographic subregion complexity for asymptotically AdS spacetimes, pertubatively around pure AdS using a variational technique. An essential step in finding subregion complexity is to identify the bulk minimal surface of the entangling subregion. Our method therefore heavily relies on solutions of an inhomogeneous version of Jacobi equation, used to study deformations of the entangling surface for perturbations of the bulk metric. Using this method we have obtained the change in complexity for a strip and a circular disk like subsystem for boosted black brane like perturbations over pure AdS_4. As a corollary, we find that for spherical subsytems in 3+1 dimensional bulk, the linear change of subregion complexity for boosted black brane like perturbations over pure AdS_4 , vanishes.
     Inhomogeneous Jacobi equation and Holographic subregion complexity
        Rohit Mishra,
        December 30, 2023
    ===================================================================
    
    
    
    § INTRODUCTION
    
    In classical computation the notion of complexity is a characterization of how hard a computational problem is? A sorting problem for example roughly requires N(N-1) steps to complete, where N is the no. of elements to be sorted. In the the asymptotic limit this is 𝒪(N^2) and the problem is therefore said to be of polynomial complexity. On the other hand a prime divisor problem is believed to require steps, which grow exponentially with the no. of digits. The problem is then termed to be of exponential complexity. Another, similar characterization of complexity is circuit complexity, which is a measure of the number of elementary gates required to go from one `classical' input to a `classical' output. An important corollary of such a finding is the fact that more complex a problem, the better a computer modeled along that problem, will be. 
    
    In the quest for a quantum computer one therefore needs a notion of complexity for quantum computations. In a quantum computer a classical bit is replaced by a qubit which is a two dimensional Hilbert space. A computation is then the action of a unitary operator on the states of this Hilbert space. A gate in quantum computation is just a unitary operator acting on a state. The notion of circuit complexity can then be suitably modified as the number of elementary unitary operations required to go from a given input state to another output state. An operational definition for this was given in <cit.>, following which the circuit complexity in quantum field theories was calculated in <cit.>.
    
    Considering the fact that calculating complexity in QFTs might be a difficult job, one considers the question of whether a simple gravitational dual of such a quantity exists for states which have a gravitational dual in the Holographic context. The proposal for such a quantity was first given in an attempt to establish the ER=EPR paradigm in <cit.>. A global AdS Schwarzschild black hole is described by two copies of CFT living on the two asymptotic ends of the geometry. While each copy of the CFT is in a thermal state, the state can equally be described by a thermofield double, which is a pure but entangled state in the product of these two Hilbert spaces. One can then assume that a measurement is done in one of the copies of the two Hilbert spaces. The complexity of this measurement affecting the state in the other Hilbert space was conjectured to be the volume of an extremal slice connecting two constant time slices, each in one of the asymptotic boundaries. This later came to be known as the Complexity=volume conjecture. To avoid certain shortcomings in the above conjecture the Complexity=action was put forward in <cit.>. The conjecture states that complexity is equal to the action evaluated in the domain of dependence (in the bulk) of the two asymptotic constant time slices. An essential improvement of this over the previous one was that it fixes an arbitrary multiplicative constant in the Complexity=volume conjecture. Recently by identifying the symplectic two form on the phase space of bulk gravity gravity with that of the symplectic structure on the projective space of states of the boundary field theory <cit.> it has been shown that a slight modification of Nielsen's prescription is indeed dual to the Complexity=volume prescription <cit.>.
    
    But this prescription holds for the state of the boundary subregion as a whole. A quantum computation may proceed by acting unitary operators on subregions of the full system. The state of the full system is then built by such operations performed one at a time. Hence it is sensible to define an analogous quantity for subregions. A working definition for this was given in <cit.>. The idea can be elaborated as follows. Suppose one considers a subregion A (defined by a spatial subregion of the boundary) of the boundary. One then needs to find a minimal co-dimension two surface 𝒮 in the bulk that is homologous to ∂ A. Recall that the area of this surface is nothing but the entanglement entropy of A. The subregion complexity is then the volume of a maximal co-dimension one surface bounded by 𝒮 and A.
    
    In this paper we will be mainly interested in this Complexity=volume definition for subregion. It is important to note that there are very few exact solutions for the minimal surface equation of codimension one/two surfaces in the asymptotically AdS spacetimes. For static spacetimes foliated by t say, it turns out that both 𝒮 and Σ lie on a t=constant co-dimension one slice and the problem is, fairly well, tractable. But once one introduces rotation, the calculation becomes difficult. This is also the case for time dependent geometries. The reason for the complication is that 𝒮 and Σ no more lie on the t=constant (here t is the foliating time for the boundary field theory) slices. One then resorts to perturbative schemes in which complexity/entanglement entropy is calculated by studying perturbations over AdS, for which the relevant surfaces are known. There are two approaches to do this. The first approach is a perturbative one, which requires solving the minimal surface equation perturbatively in some small approximation <cit.>. One then uses this approximate solution to evaluate the area or volume order by order<cit.>. The second approach is a variational one. It starts by taking the variation of the area and volume functional that incorporates both changes in embeddings and metric perturbation <cit.>. The changes in the embedding is obtained by solving an inhomogeneous Jacobi equation, the inhomogeneity of which is sourced by metric perturbations. It can be shown that in some known cases, the first and second variation of the area functional, obtained using the variational approach, exactly matches with the first and second order results, obtained  using the perturbative approach <cit.>. In this paper we extend this variational approach from area functional to the volume. It turns out that the differential equations that one  
    needs to solve, for calculating subregion complexity are all of second degree and second order. The first for finding out the surface 𝒮 that determines the entanglement entropy and the other for finding out the maximal slice Σ (say). 
     
    In this paper we provide a general systematic scheme for performing the variational expansion. Since the definition of subregion complexity requires finding the bulk entangling surface of the subregion, the scheme heavily relies on finding how the entangling surface deforms under the perturbation of the bulk geometry. A definite procedure for finding such deformations was given in <cit.>. They are given by a vector field which has support on the entangling surface 𝒮 and arise as a solution of a linear elliptic differential equation, which is nothing but an inhomogeneous version of the well known Jacobi equation for minimal surfaces. We find that the first order change in the holographic subregion complexity is nothing but the integral of a component of this vector field over 𝒮. Using the solutions obtained in <cit.> we compute the first order changes of complexity for BTZ like perturbations over AdS_2+1 and boosted black brane like perturbations over AdS_3+1.  
    
    We also derive general expressions for finding the second variation of complexity. As, such a calculation requires a separate discussion, we do not explicitly solve any specific examples for this case. However we do point out the advantages of finding such general expressions. In particular, how these can be put to use in the context of both holographic complexity as well as entanglement entropy, is discussed. In fact the general expressions obtained would allow one to calculate the change in complexity as well as entanglement entropy for variations of the subregion as well, i.e in cases where the boundary state does not change, rather the boundary subregion is deformed. 
     
    
    
    § NOTATIONS AND CONVENTIONS
    
    The theory of perturbation to be used here will essentially be the one in the covariant framework introduced in <cit.>. The perturbation of d+1 dimensional space time (ℳ,g) to another d+1 dimensional space time (ℳ',g') is characterized by a differentiable map Φ:ℳ→ℳ' which is however not isometric.  (ℳ',g') is then called a perturbation over (ℳ,g) if (1)P=Φ_*g'-g is a small perturbation over g. 
    
    
    Consider a surface 𝒮 isometrically embedded in ℳ and given by the function f:S→ℳ. Let the coordinates on 𝒮 be denoted by τ^a. If 𝒮 has a boundary ∂𝒮 then let its coordinates be denoted by χ ^A. This will be the general strategy followed in the text. While discussing an embedded surface, its coordinates will be denoted by τ^a and the coordinates on its boundary will be denoted by χ^a. The ambient spacetimes coordinates will always be denoted by x^μ. The induced metric on 𝒮 is given by h=f_* g, which, in the local coordinates can be written as h_ab=g(∂_a,∂_b)=∂ x^μ/∂τ^a∂ x^ν/∂τ^bg(∂_μ,∂_ν). The quantity ∂ x^μ/∂τ^a∂_μ is the push forward of the purely tangential vector field ∂_a to ℳ. `h_ab' is the first fundamental form on 𝒮.  One can decompose the tangent space at the point x ∈ 𝒮 into the tangent space of 𝒮 and the space of normal vectors as T_xℳ=T_x𝒮⊕ T_x^⊥𝒮.  Let the covariant derivative on 𝒮 be denoted by D: T𝒮⊗ T𝒮→ T𝒮. Let X,Y ∈ T𝒮. Then the Gauss decomposition allows us to write,
    
        ∇_XY=D_XY+K(X,Y),
    
    where D_XY is purely tangential and K(X,Y) is a vector in the normal bundle and is the extrinsic curvature or the second fundamental form. A connection ∇^⊥_XN^⊥ in the normal bundle can be defined as ∇^⊥:T𝒮⊗ T^⊥𝒮→ T^⊥𝒮, where X∈ T𝒮 and N^⊥∈ T^⊥𝒮. The shape operator W_N^⊥(X) is then defined as,
    
        ∇_XN^⊥=∇_X^⊥ N^⊥-W_N^⊥(X).
     
    The Weingarten equation then relates the shape operator and the extrinsic curvature,
    
        g(W_N^⊥(X),Y)=g(N^⊥,K(X,Y)),
    
    where X,Y∈ T𝒮 and N^⊥∈ T^⊥𝒮. The Riemann tensor is defined as,
    
        R(W,U)V:=[∇_W,∇_U]V-∇_[W,U]V
    
    Similarly one can define the Riemann tensor on 𝒮 as,
    
        ℛ(X,Y)Z:=[D_X,D_Y]Z-D_[X,Y]Z
    
    The Gauss and Codazzi equations relates the intrinsic Riemann tensor and the spacetime Riemann tensor. Let X,Y,Z,W∈ T𝒮 and N^⊥∈ T^⊥𝒮. Then the Gauss equation is given as,
    
        g(R(X,Y)Z,W)=g(ℛ(X,Y)Z,W)-g(K(X,Z),K(Y,W))+g(K(X,W),K(Y,Z)),
    
    and the Codazzi equation as,
    
        g(R(X,Y)N^⊥,Z)=g((∇_YK)(X,Z),N^⊥)-g((∇_XK)(Y,Z),N^⊥)
    
    Due to the presence of perturbations a variation will have two contributions, one which is a flow along a vector N∈ Tℳ, and another variation δ_g which is purely due to metric perturbations. The metric perturbation will be given by,
    
        (δ_g g)(∂_μ,∂_ν):=[Φ_* g'-g](∂_μ,∂_ν)=(1)P(∂_μ,∂_ν),   (δ_g^2g)(∂_μ,∂_ν)=(δ_g (1)P)(∂_μ,∂_ν):=(2)P(∂_μ,∂_ν)
    
    where (1)P is a symmetric bilinear form on ℳ. Note that δ_g only acts on the metric and does not change the vector fields ∂_μ. Now suppose there is a covariant derivative ∇' in ℳ' compatible with g', then for X,Y∈ Tℳ,
    
        (1)C(X,Y):=δ_g(∇_XY)=∇̃_XY-∇_XY,
    
    where ∇̃=ϕ^*∇' is the pullback connection on M. It is important to note that C(X,Y) is a vector field on ℳ.
    
        δ_g (1)C(X,Y)=(2)C(X,Y)-2 (1)P((1)C(X,Y)),
    
    where (1)P(X) is a vector defined as g((1)P(X),Y)=(1)P(X,Y).
    
     
    
    § VARIATIONS ON GENERAL MINIMAL SURFACES
    
    
     §.§ Variation of area
    
     In this section we will calculate the variation of area. These expressions are in general true for any minimal submanifold. Note that these have been obtained before, but certain boundary terms have been ignored <cit.>. We will retain the boundary terms, obtained in the process, and later see how these boundary terms are useful. Only the first order expressions will be used to find the perturbative change of holographic subregion complexity for boosted black brane like perturbations over pure AdS  in section (<ref>). The second order expressions will however not be evaluated for the examples taken but the usefulness of these expressions in obtaining properties of subregion complexity and holographic entanglement entropy will be discussed. We will also discuss the scheme which one must follow in order to evaluate the second order variation, given a perturbation. 
     
    Let us consider a one parameter flow (generated by the vector field N) of the area and an n-dimensional minimal surface 𝒮 embedded in a d+1 dimensional spacetime. First variation of area of the submanifold is given by,
     
        δ_N A=-∫_𝒮g(N^⊥,K(∂_a,∂_b))h^ab) √(h)d^nτ+1/2∫_𝒮(1)P(∂_a,∂_b)h^ab√(h) d^nτ-∫_∂𝒮 g(r̂, N^T)√(γ) d^n-1χ,
    
     where `τ^a' are coordinates on 𝒮 and `χ^A' on ∂𝒮. h_ab and γ_AB are the induced metrics on 𝒮 and ∂𝒮. The unit co-normal to ∂𝒮 w.r.t 𝒮 is denoted by r̂.
     K(∂_a,∂_b) is the extrinsic curvature of 𝒮 in ℳ. ⊥ and T denotes component perpendicular and parallel to 𝒮, respectively. Further N has been taken to commute with ∂_a's. The first term in eq. (<ref>) is however zero for a minimal surface. While deriving the second variation we will differentiate the above expression rather than the one obtained after imposing the minimal surface condition. Thus we have,
     
     
        δ^2_NA=∫_𝒮 d^nτ √(h)(h^abh^cdP(∂_b,∂_d)g(N^⊥,K(∂_a,∂_c))-h^abg(C(∂_a,∂_b),N^⊥))
        
         +∫_𝒮 d^nτ √(h)[h^ab 2(2)P(∂_a,∂_b)-1 2h^ach^bd(1)P(∂_a,∂_b)(1)P(∂_c,∂_d)+1 4h^abh^cd(1)P(∂_c,∂_d)(1)P(∂_a,∂_b)]
        
         -∫_∂𝒮 d^n-1χ √(γ)r̂_a[h^abh^cd(1)P(∂_c,∂_d)g(N,∂_b)]-2∫_∂𝒮 d^n-1χ √(γ)r̂_a[h^ab(1)P(N^⊥,∂_b)]
        
         -∫_∂𝒮 d^n-1χ √(γ)r̂_a[h^abg(N^⊥,∇^⊥_∂_bN^⊥)]
         -∫_∂𝒮 d^n-1χ √(γ)g((∇_NN)^T,r̂),
    
     The details of this calculation can be found in appendix <ref>. In section <ref> we will discuss how to evalute individual terms in the above expression. In other words we will discuss what decides the individual terms arising in the above expression. Before delving into such a discussion we will first derive the second order perturbation equation. This will help us to discuss the origin of terms like (∇_NN)^T in the above expression.
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
    
     §.§ Second order perturbation equations
    
    
      §.§.§ First order perturbation revisited
    
     Assume that under a perturbation of the ambient metric the minimal surface 𝒮 deforms to a minimal surface 𝒮'. Let the deviation be denoted by the vector field N.
     We will recall the derivation of the first order inomogeneous Jacobi equation, for 𝒮, to facilitate further discussion. Finding the first order Jacobi equation amounts to equating the total variation of the mean curvature vector H to zero, i.e,
     
        ℒ_N H+δ_g H=0
    
     However note that ℒ_N H =∇_N H-∇_H N. Since H=0 on the back-ground one only needs to calculate ∇_N H. We will see that this will also be the case for the second order perturbation equations. Let us start by calculating ∇_N∇_∂_a∂_b. N will be taken to commute with all the ∂_a's. However note that neither N^T nor N^⊥ individually commute with ∂_a's. In fact [N^T,∂_a]=-[N^⊥,∂_a] and [N^T,∂_a]^⊥=-[N^⊥,∂_a]^⊥=0. We will rewrite K(∂_a,∂_b) as ∇_∂_a∂_b-(∇_∂_a∂_b)^T.
     
        h^ab∇_N∇_∂_a∂_b=R(N,∂_a)∂_b+∇_∂_a∇_∂_bN
    
     The derivative of the tangent part gives,
     
        h^ab∇_N(g(∇_∂_a∂_b,∂_c)h^cd∂_d)=(R(N,∂_a)∂_b+∇_∂_a∇_∂_bN)^T+g(H,(∇_∂_cN)^⊥ )h^cd∂_d+(∇_(∇_∂_a∂_b)^TN)^⊥
    
     The ∇_N h^ab K(∂_a,∂_b) term gives,
     
        ∇_N h^ab K(∂_a,∂_b)=-h^cdK(∂_a,∂_c)[g(∇_∂_bN^⊥,∂_d)+g(∂_b,∇_∂_dN^⊥)]-2h^aeh^bf(g(∂_e,D_∂_fN^T))K(∂_a,∂_b)
    
     Meanwhile δ_gH= (1)C(∂_a,∂_b))^⊥-h^ach^bdK(∂_a,∂_b)(1)P(∂_c,∂_d). Using the identities in section(<ref>) and using the above expressions along with the results of appendix (<ref>), one can write δ_N H as, 
     
        h^ab(B_ab-D_ab+C(∂_a,∂_b)+S_ab-h^cdK(∂_a,∂_c)P(∂_d,∂_b))^⊥-g(H,(∇_∂_cN)^⊥ )h^cd∂_d+∇_N^TH,
    
     where the quantities B_ab,D_ab and S_ab are given by
     
        B_ab:=R(N^⊥,∂_a)∂_b+∇_∂_a∇_∂_bN^⊥,     D_ab:=∇_(∇_∂_a∂_b)^TN^⊥,
           S_ab:=-h^cdK(∂_a,∂_c)[g(∇_∂_bN^⊥,∂_d)+g(∂_b,∇_∂_dN^⊥)]
    
     For convenience we will define one more quantity L_ab, which will be useful later and is given as 
     
        L_ab=(B_ab-D_ab+C(∂_a,∂_b)+S_ab-h^cdK(∂_a,∂_c)P(∂_d,∂_b))
    
     The appearance of the ∇_N^TH might have have looked obvious but is not quite so. Due the fact that neither N^T nor N^⊥ commute with the ∂_a's there could have been terms, depending on the commutators, surviving. It however turns out that the commutator terms do not contribute (appendix <ref>). Since H=0 everywhere on 𝒮 it follows that ∇_N^TH =0. The resulting form can then be shown to reproduce the correct form of the Jacobi equation once the minimal surface condition H=0 is imposed. 
     
     
        ℒN^⊥=-C^⊥+H̃.
    
     where
     
        ℒ(N^⊥)=Δ^⊥N^⊥+Ric(N^⊥)+A(N^⊥),
    
     where we have defined  Δ^⊥N^⊥ to be the Laplacian on the normal bundle, given by h^ab(∇^⊥_∂_a∇^⊥_∂_bN^⊥-∇^⊥_(∇_∂_a∂_b)^TN^⊥), g(R(N^⊥,∂_a)∂_b,N^⊥) has been denoted by Ric(N). A(N)=h^abK(∂_a,W_N(∂_b)) is the Simon's operator. Whereas  C^⊥ is defined as C^⊥=h^ab(1)C(∂_a,∂_b)^⊥ and H̃=(1)P^abK(∂_a,∂_b).
     Thus identifying the Jacobi/Stability operator (ℒ) for minimal surfaces. This equation was derived in <cit.> in the context of the co-dimension two entangling surface, but will be applicable our case as well where the complexity surface is also minimal. In fact this equation is applicable to any embedded minimal submanifold ( 𝒮 say ) and the solutions of the above equation N^⊥ maps a minimal submanifold (𝒮) to another submanifold (𝒮' ) which is minimal in a perturbed geometry. Hence, the superscript ⊥ here implies components normal to 𝒮 and T implies components tangential to 𝒮.
     
     
     The form eq. (<ref>) will be taken to perform the second variation. The reason for taking this and not the known form is the following. Inside a derivative one cannot assume N^⊥ to be normal, hence the term S_ab cannot be written in terms of the extrinsic curvature inside a derivative. If one is dealing with a family of co-dimension two surfaces then taking N^⊥ inside a derivative amounts to making an assumption that N^⊥ is normal throughout the foliation. The simple analogue of this statement is the following. g(N^⊥,∂_a) is equal to zero, but its normal derivative ∇_N^⊥g(N^⊥,∂_a) may not be.
     
     
    
      §.§.§ Second order perturbation equations
    
     The second order equation is obtained by equation the second total variation of H to zero, i.e,
     
        (ℒ_N H+δ_g)(ℒ_N H+δ_g) H=0
    
     One can show that this is nothing but,
     
        (∇_N H+δ g)(∇_N H+δ g) H -R(N,H)N-∇_[N,H]-∇_δ_g HN-∇_Hδ_g N=0,
    
     where the condition ℒ_N H+δ_g N=0 condition has been used for the back-ground. It is evident from the the above expression that on further this condition along with the condition H=0 on the back-ground yields the equation,
     
        (∇_N H+δ g)(∇_N H+δ g) H=0
    
     Note that the equation H=0 holds everywhere on 𝒮. Therefore ∇_N^T∇_N^T H=0. Further, since the inhomogeneous Jacobi equation holds everywhere on 𝒮, it's tangent derivative is also zero. Therefore while finding the second variation one only needs to consider the normal component N^⊥ of the deviation vector and not the tangent part. The final expression that needs to be calculated is therefore (∇_N^⊥+δ_g)(L_ab^⊥). Further note that δ_g N has been taken to be zero. This is because had this been taken to be non-zero one would have got an equation for ∇_N^⊥N^⊥+δ_g N^⊥ rather than ∇_N^⊥N^⊥. Every expression including the area variations would then have contained  ∇_N^⊥N^⊥+δ_g N^⊥ rather than ∇_N^⊥N^⊥. Hence δ_g N is just a gauge redundancy and can be set to zero. Everything put together given the following equation,
     
     
     
     
    
        h^ab(∇_∂_a∇_∂_bJ^⊥-∇_(∇_∂_a∂_b)^TJ^⊥+R(J^⊥,∂_a)∂_b+2h^cdK(∂_c,∂_a)g(J^⊥,K(∂_b,∂_d)))^⊥
        
        +h^abh^cdL_ac^⊥[4g(N,K(∂_b,∂_d))-2P(∂_b,∂_d)]-2h^ab∇_(L_ab)^TN
        
        -h^abh^cdK(∂_c,∂_d)[2g(R(N,∂_b)N,∂_d)+2g(∇_b N,∇_d N)]
        
        h^ab[4R(N,∂_a)∇_∂_bN+(∇_∂_aR)(N,∂_b)N+(∇_N R)(N,∂_a)∂_b]
        
        h^ab[2(∇_N C)(∂_a,∂_b)+4C(∇_∂_aN,∂_b)]-4h^abh^cdK(∂_a,∂_c)[g(C(N,∂_b),∂_d)+P(∇_∂_b N,∂_d)]
        
        h^ab(2)C(∂_a,∂_b)^⊥-2P(C(∂_a,∂_b)^⊥)-h^abh^cd(2)P(∂_b,∂_d)K(∂_a,∂_c)=0.
    
    In order to find out the new minimal surface in ℳ' one needs to find out what N on 𝒮 is, which is provided by the solution to eq. (<ref>) and all derivatives of N along N itself, on 𝒮. It is precisely this that these perturbation equations provide. The second order perturbation equation eq. (<ref>) for example is providing what (∇_NN)^⊥ is. Further higher order equations can in general be obtained, but is not useful to our discussion.
     
     
     
    
     §.§ The perturbation scheme and solution procedure
    
     The scheme for calculating the expressions eq. (<ref>, <ref>) is the following. The first term in the first order variation eq. (<ref>) vanishes by equations of motion. Hence, one has to calculate the second term and the boundary term. The second term can easily be calculated as the embedding of 𝒮 and the perturbation (1)P(∂_a,∂_b) are known. Calculating the boundary term however requires more information. Mathematically one requires N|_∂𝒮. This  information is provided by how the boundary ∂𝒮 deforms. Let the vector field η characterise this deviation. Therefore one has N|_𝒮=η. The Jacobi equation on 𝒮 is an expression for N^⊥. The boundary conditions for solving this equation is of the Dirichlet type i.e N^⊥|_𝒮 is given. It is precisely this boundary condition that is provided by η^⊥, while η^T can be used for calculating the change in volume using eq. (<ref>). At second order one has to evaluate the terms in eq. (<ref>). All the terms except the last one in this expression can be obtained using the expressions for N^⊥, which are provided by solution of the first order IJE for 𝒮. The last term requires a discussion of the second order IJE. The second order IJE provides an expression for (∇_N N)^⊥ the boundary condition for which is fed from ∇_ηη i.e (∇_N N)^⊥|_∂𝒮=(∇_ηη)^⊥ while (∇_ηη)^T  is used in eq. (<ref>).
     
     Let us now discuss the solution procedure for the special case we will be dealing with. Note η is completely uspecified. Now, consider the special case where ∂𝒮 is also a minimal surface, but ofcourse, of one dimesion less than 𝒮. Assume that under the perturbation 𝒮 deform to a minimal surface 𝒮' such that ∂𝒮 also deforms to a minimal surface ∂𝒮'. Then it is evident that η is a solution of an IJE but now for ∂𝒮 rather than of 𝒮 i.e,
     
     
        ℒN^⊥=-C^⊥+H̃.
    
     where
     
        ℒ(N^⊥)=Δ^⊥N^⊥+Ric(N^⊥)+A(N^⊥),
    
     where we have defined  Δ^⊥N^⊥ to be the Laplacian on the normal bundle of ∂𝒮, given by γ^AB(∇^⊥_∂_A∇^⊥_∂_Bη^⊥-∇^⊥_(∇_∂_A∂_B)^Tη^⊥), g(R(η^⊥,∂_A)∂_B,η^⊥) has been denoted by Ric(η). A(η)=γ^ABK(∂_A,W_η(∂_B)) is the Simon's operator. Whereas  C^⊥ is defined as C^⊥=γ^AB(1)C(∂_A,∂_B)^⊥ and H̃=(1)P^ABK(∂_A,∂_B). All these quantities now refer to ∂𝒮 as opposed to 𝒮 and ⊥ now refers to perpendicular to ∂𝒮 rather than 𝒮. Similarly ∇_ηη is obtained from the second order IJE eq. (<ref>) but now written for ∂𝒮 rather than for 𝒮 by making the same changes as done for eq. (<ref>).
     
     
     
     
     
     
     
     
      
    
    
    § PERTURBATIVE CHANGE OF SUBREGION COMPLEXITY
    
    Given an asymptotically AdS spacetime and assuming that it is dual to a CFT, the holographic entanglement entropy of a subregion A on a fixed time slice in the boundary CFT is given according to the Ryu-Takayanagi proposal <cit.>, by the area of an extremal (minimal surface in pseudo-Riemannian geometry [From now on we will use extremal and minimal interchangeably.]) co-dimension two surface (𝒮) in the bulk, which is homologous to the boundary ∂ A of the subregion A. 
    The holographic complexity is given by the volume (V) of an extremal co-dimension one minimal surface (Σ) with 𝒮∪ A as its boundary <cit.>. 
    
        C_V =V/8π R G
    
    
    For a static geometry the t=constant slices (T say) are maximal and hence both 𝒮 and Σ lie on T and therefore all the calculations can be done on the T's. For non-static or non stationary bulk spacetimes this may not be the case. In such situations one uses the covariant Hubeny-Rangamani-Takayanagi (HRT)<cit.> proposal to calculate the holographic entanglement entropy. It is known from the variational approach that the contribution from the changes in embedding does not appear at the first order of metric perturbation <cit.>. Hence one can work with the same fixed time slice in non static spacetimes for calculating HEE upto first order. In this section we develop a similar variational scheme to handle this situation for the volume of subregions. Unlike area the first variation of the volume functional incorporates the changes in the embeddings. However in the next section we will see that the information of the deviation from the t=constant time slice does not appear in the first variation. Hence one can still work with the fixed time slice for the volume change in non static backgrounds.
    
    Let us denote by 𝒮' and Σ ' the entangling surface and the complexity surface in a perturbed geometry. We will be coordinatizing both the primed and the unprimed surfaces by the same coordinate system. Σ will now play the role of 𝒮 and S∪ A the role of ∂𝒮 of section <ref>. The area of 𝒮 is nothing but the volume (V) of the spacelike surface Σ. Σ ' can then be characterised by a deviation vector N between Σ and Σ'. Let the restriction of N to 𝒮∪ A be denoted by η. Since, the boundary subsytem is not deformed it is evident that η|_A=0. Further η must be such that it maps the minimal surface 𝒮 to the minimal surface S' in the perturbed metric. The change in volume of the complexity surface is then given by,
    
        δ_N V=-∫_Σ g(N,K^(Σ↪ℳ)(∂_a,∂_b))h^ab√(h) d^dτ+1/2∫_Σ(1)P(∂_a,∂_b)h^ab√(h) d^dτ      -∫_𝒮 ∪  A g(r̂, N)√(γ) d^d-1χ,
    
    where τ^a are the coordinates intrinsic to Σ, χ^A are the coordinates intrinsic to 𝒮, K^(Σ↪ℳ)(∂_a,∂_b) is the extrinsic curvature of Σ in ℳ and r̂ is the inward directed co-normal of 𝒮 w.r.t Σ. Note that the first term does not contribute since Σ is extremal. Moreover N, in the last term, can be replaced by its restriction η to 𝒮. The component of N perpendicular to Σ, in the above expression, is obtained by solving the inhomogeneous Jacobi equation (IJE) eq. (<ref>).
    
     To make things clear we will put another subscript to ⊥ and T (⊥_Σ, T_Σ for example) to indicate the surface to which the the components are perpendicular and tangential. As has been emphasized in section <ref>, this eq. (<ref>) does not determine the tangential component N^T_Σ. N^T_Σ which is known only at the boundary 𝒮∪ A and is obtained by how 𝒮 and the subregion A deforms. It can also be shown that the volume change does not depend on how N^T_Σ is extended to Σ. It is clear from the situation at hand that A does not deform and the deformation is non-zero only on 𝒮. However note that both 𝒮 and 𝒮' are co-dimension two minimal surfaces in ℳ ans ℳ'. Hence η:=N|_𝒮 takes a co-dimension two minimal surface in ℳ to one in ℳ'. Therefore, η^⊥_𝒮 must satisfy an inhomogeneous Jacobi equation (<ref>) for a co-dimension two minimal surface 𝒮. Note that η can be taken to be completely perpendicular to 𝒮. This is because the boundary of 𝒮 which is nothing but ∂ A, does not deform under the perturbation as the boundary subregion is kept intact. Since the normal bundle of 𝒮 is two dimensional, there are essentially two independent components of η. One of the components will be perpendicular to both Σ and 𝒮 and the other would be perpendicular to 𝒮 but tangent to Σ. It is precisely the second component that is along the co-normal r̂. Hence, (η)^⊥_Σ=(η^⊥_𝒮)^⊥_Σ:=N^⊥_Σ|_𝒮 determines the boundary condition for eq. (<ref>) applied to Σ whereas η^T_Σ=(η^⊥_𝒮)^T_Σ determines the change in volume in eq. (<ref>). The solutions η for 𝒮 was obtained for certain cases in <cit.>. We will essentially be using these solutions to determine the change in volume from eq. (<ref>). 
    
    In the next few subsections we will be evaluating the expression eq. (<ref>) for several examples, for which the deviation vector has been calculated in <cit.>. First we will consider the 2+1 dimensional case in which we will be considering rotating BTZ as a perturbation over AdS_3 as an example. Next, for the 3+1 dimensional case we will be considering boosted black brane like perturbations over AdS_4. 
    
    
    
     §.§ 2+1 dimensional case
    
    Consider the AdS_3 metric in the Poincaré patch as the unperturbed background and consider the planar rotating BTZ as a perturbation over it. The AdS_3, metric is given by
    
    ds^2=-dt^2+dz^2+dx^2/z^2𝒮 in this case is given by a spacelike geodesic in AdS_3. The equation for such a geodesic, which is homologous to an interval on a spacelike slice of the boundary is given by,
    
    z(τ)=z^(0)_*(τ),   x(τ)=z^(0)_*tanh(τ),
    .
    z_*^(0) is the AdS turning point. The parameter τ ϵ  (-∞,∞) covers the full geodesic.  Since the geodesic is affinely parametrised the metric on 𝒮 is just identity. The metric on Σ
    can be conveniently written in a co-ordinate system that includes τ and another co-ordinate c (say). To understand how c is defined note that z_*^(0) is a constant on a particular geodesic. However this can be lifted to the coordinate c (say) so that c, τ cover the region bounded by 𝒮 and A on a t=constant slice (Σ) [This follows from the nesting property of entanglement entropy <cit.>. In other words as one reduces the subregion size the bulk geodesics do not intersect each and therefore the turning point behaves monotonically and qualifies for a good coordinate.]. In other words we are foliating this region in such a way that, each leaf is a spacelike geodesic with a particular turning point, given by the value of c on that leaf. The metric h on Σ is then given by,
    
    
        h_ab dτ^adτ^b=dc^2+(τ)^2 dτ^2/c^2^2τ.
    
    
    The deviation vector η can be calculated by solving an inhomogenous Jacobi equation for the spacelike geodesic in AdS_3. Given the tetrads that are parallely propagated along the geodesic,
    
    e^μ_0=(z,0,0),  e^μ_1=(0,± z√(1-(zz^(0)_*)^2),z^2z^(0)_*),  e^μ_2=(0,z^2z^(0)_*,∓ z√(1-(zz^(0)_*)^2)),
    η can be written as η^I e_I^μ, with,
    η^0   =   1/3b z^(0)_*^2 tanh(τ)(τ)
    η^2   =   -1/6 a z^(0)_*^2cosh(2τ)(τ)^3
    
    Note that the co-normal r̂ in eq. (<ref>) is nothing but -e^μ_2. Hence the last integral in eq. (<ref>) reduces to integrating η^2 over the spacelike geodesic. The integral evaluates to,
    
    ∫ g(r̂, η)√(h) d^d-1τ=π/4 a z^(0)_*^2
    
    To evaluate the first term, we will use the τ, c coordinatisation.  Then calculating the first term reduces to calculation the integral.
    
    a∫_-∞^∞∫_0^z^(0)_*c^3(τ) dc dτ=π/4 a z^(0)_*^2
    
    The total change in volume can then be calculated as,
    
        Δ V=π/8a z^(0)_*^2
     
    
    
    The only question remaining is whether the inhomogenous Jacobi equation for the co-dimension one minimal surface admit a solution with the boundary conditions given by solutions of the co-dimension two inhomogeneous Jacobi equation. To see it, let us write down the equation first. The normalized normal one form of the t=constant slices are given by dt/z, and the normal vector is given by z∂_t. Let us denote it by n̂. Since this is the only normal and is also normalized, it follows that the ∇_a^⊥n̂=0. Now, assume that the solution is of the form αn̂. The codimension one inhomogeneous Jacobi equation, for rotating BTZ like perturabtions reduces to,
    
    
    
    
        z^2(∂_x^2 α+∂_z^2 α)-2α=0.
    
    Note that the inhomogeneous term is in fact zero for this case. This can be checked from the explicit form of (1)C^μ _νρ given in appendix (<ref>).
    To facilitate imposing boundary conditions on the entangling surface we will work in coordinates c,  τ as explained before.
    
    
    
     
        ^2τ(y^2∂_y^2α+y∂_yα)+∂_τ^2α+tanhτ∂_τα-2α=0
    
    
    
    
     Boundary conditions,
     
        α(y_0,τ)=b/3y_0^2tanhττ,   α(y,τ=∞)=0,   α(y,τ=-∞)=0
    
     Unfortunately this equation does not admit a solution by separation of variables. We therefore try to solve eq. (<ref>) on the disc x^2+z^2=1 with the first boundary condition only (The actual solution is however required on the half disc with the full set of boundary condition). If the solution satisfies the second and the third boundary condition then we can conclude that a solution on the half disc with full set of boundary conditions does exist. The plots of the solutions are given in fig. (<ref>). From the figure it is clear that this indeed is the case.
     
     
     
    
    
     §.§ 3+1 dimensional case
    
    For the 3+1 dimensional case we will consider the Boosted AdS black brane like perturbations over AdS_4. The metric and the asymptotic expansion of the metric is given in appendix <ref>. The choice of the back ground metric is therefore the AdS_4 metric in the Poincaré patch. The AdS length is set to unity and the metric is given by,
     
        ds^2= -dt^2+dx^2+dy^2+dz^2 z^2
    
     We will now consider a circular disk like subregion and calculate the change in holographic complexity from pure the AdS value.
     
     
    
      §.§.§ Circular disk subregion
    
    For the case of a circular boundary subregion of radius ℛ, the minimal surface in the AdS_d+1 is a d-1 dimensional hypersphere. Such a surface can be described by the following embedding functions, which are just a parametric form  of the equations obtained in <cit.>,
     
        x=ℛsinθcosϕ+X,  y=ℛsinθsinϕ+Y,  z=ℛcosθ,  t=constant.
    
    The coordinates θ, ϕ have ranges, 0≤θ≤π/2 and 0≤ϕ<2π. It is also clear from eq.(<ref>) that the half sphere intersects the AdS_4 boundary at θ=π/2. The intrinsic metric on 𝒮 is then given by,
    
        γ_AB dχ^Adχ^B=dθ^2 +sin^2θ  dϕ^2cos^2θ,
    
    while the metric on Σ in ℛ, θ, ϕ coordinates is given by,
    
        h_ab dτ^adτ^b=dℛ^2+ℛ^2(dθ^2+sin^2θ dϕ^2)/ℛ^2cos^2θ
    
    
    
    
    
    We will set up a local orthonormal basis on this surface to facilitate our calculations. Since the surface is space-like, the vectors
    
        e_2=cosθ∂_θ,  e_3=θ∂_ϕ.
    
    provide us with two space like bases. The other two, spanning the normal bundle will provide us with the other two basis vectors. As a matter of convention we mark the time like normal as e_0 and the space like normal as e_1.
    
        e_0=z∂_t,  e_1=z(x-X)ℛ∂_x+z(y-Y)ℛ∂_y+z^2ℛ∂_z
    
    
    The solution to the co-dimesion two inhomogeneous Jacobi equation has been obtained in <cit.> and can be represented as the components of the deviation vector in the directions normal to the co-dimension two surface eq. (<ref>) [The solution given in <cit.> has an error. The authors missed a factor of 1/2 in the first term of the inhomogeneous part of eq. (42) of <cit.>.],
    
        η^1=ℛ ^3/96(3 β ^2 γ ^2+2) cos ^2θ(cos 2 θ-5)-ℛ ^3/64cos 2 ϕ(β ^2 γ ^2 sin ^22 θ)
    
    
    
    
    
    
        η^0=-1/4βγ ^2 ℛ^3 sinθcos ^2θcosϕ
    
    The third term in eq. (<ref>) is then just the integral of -N^1 over 𝒮 and yields,
    ∫ g(r̂, η)√(h) d^d-1τ =π/9z_0^3ℛ^3 (3 β ^2 γ ^2+2) 
    
    The second term in eq. (<ref>), can be calculated by lifting ℛ to a coordinate r. The expression then reduces to,
    3 β ^2 γ ^2+2/z_0^3∫_0^ℛ∫_0^π/2∫_0^2πr^2sin(θ)dr dθ dϕ=2π/9z_0^3ℛ^3(3 β ^2 γ ^2+2)
    
    The total change in complexity is then given by,
    Δ C=0
    
    Thus the first order correction to the holographic complexity vanishes for a spherical entangling surface. This is consistent with the result for AdS black brane (β→ 0) results obtained in <cit.>. However, the finite β case has not been dealt with before. Hence, we can conclusively say that even for a boosted black brane like perturbations over AdS the first order change of subregion comlexity is zero fro a spherical subregion.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
     
    
      §.§.§ Thin Strip subregion
    
     Let us now consider a two dimensional strip like subregion on the AdS_4 boundary, given by the region [-L,L]×[-l/2,l/2] of the x-y plane, where L≫ l. The minimal entangling surface in this case, given in <cit.>, can be characterized by the following embedding functions,
     
        x=λ,  y(θ)=-z_* E(π-2θ 4.|2),  z(θ)=z_*√(sinθ),
    
    where z_* is the turning point of the minimal surface in AdS_4 and E(α,β) is the incomplete elliptic integral of the second kind. The range of the coordinates  are 0≤θ≤π and -L≤λ≤ L. The intrinsic  metric of 𝒮 therefore takes the form
    
        γ_AB dχ^Adχ^B=z_*^2dθ^2+4sinθdλ^2 4 z_*^2sinθ^2.
    
     Further, the turning point z_* can be written in terms of the width l of the subregion as z_*=Γ(1/4)l/2√(π)Γ(3/4).The metric on Σ can be written in a co-ordinate c, λ, θ, where c is obtained by elevating the constant z_* to a coordinate. We again choose a tetrad adapted to the 𝒮 as,
    
        e_2=2sinθ∂_θ,  e_3=z_*√(sinθ)∂_λ,   e_1=z(sinθ∂_z-cosθ∂_y),  e_0=z∂_t
    
    
      §.§.§  Case A: Strip along `x' boost along `x'
    
    The component of the solution that will be important for our case is η^1, which is given in <cit.>,
    
        η^1=C_1cos(θ)/√(sin(θ))+C_2sinθ  _2F_1(1/4, 1; 1/2, cos^2θ)
        
        -cos(θ)/√(sin(θ))∫_^θ[1/4(3D+3C/2)z_*^3(sinθ)^1/2-7/8Dz_*^3(sinθ)^5/2](sinθ')^3/2  _2F_1(1/4, 1; 1/2, cos^2θ')dθ'
        
        +sinθ  _2F_1(1/4, 1; 1/2, cos^2θ)∫_^θ[1/4(3D+3C/2)z_*^3(sinθ)^1/2-7/8Dz_*^3(sinθ)^5/2]cos(θ')dθ'
    
    the constants of integration being given by C_1=π z_*^3/16(2 C+D) and C_2=-Γ(1/4)^2 z_*^3(2C+D)/16 √(2 π) and constants C,D are given by C=(1 3+β^2γ^2)1 z_0^3, D=1 31 z_0^3 (appendix <ref>). The integral over Σ can be evaluated in c, λ, θ coordinate. The integral to be evaluated then reduces to,
    
    
        (2/3+β^2γ^2)/z_0^3∫_-L^L∫_0^z_*∫_0^π c (1/2√(sinθ)√(1-2 sin ^2(1/4 (π -2 θ)))+cosθ E(.1/4 (π -2 θ)|2)/2 √(sinθ))dc dθ
        
        =2Lz_*^2/z_0^3(2/3+β^2γ^2)
    
     
    The surface integral is just the integral of η^1 over the co-dimension two minimal surface (𝒮). Let us evaluate the integrals one at a time. The first term in the homogeneous part can be calculated easily and gives,
    
        ∫_-L^L ∫_0^πC_1cos(θ)/√(sin(θ))1/2z_* sin^3/2(θ)dθ dλ=-LC_1/z_*[1/sinθ]^θ→π_θ→ 0
    
    The second term in the homogeneous part has to be evaluated by an integration by parts. An integration by parts yields the following expression,
    
        ∫_-L^L ∫_0^π C_2sinθ  _2F_1(1/4, 1; 1/2, cos^2θ) 1/2z_* sin^3/2(θ) dθ=LC_2/z_*[-sin^3/2θ  _2F_1(1/4, 1; 1/2, cos^2θ)cosθ/sinθ]^θ→π_θ→ 0
        
        +LC_2/z_*∫_0^πd/dθ(sin^3/2θ  _2F_1(1/4, 1; 1/2, cos^2θ))cosθ/sinθdθ
    
    There are divergences in the integrals of both the homogeneous terms. We will therefore not evaluate it right away, but show how the divergences cancel with the terms arising from the particular solution. To evaluate the pieces in the particular solution, let us define the function,
    
        f(θ)=[1/4(3D+3C/2)z_*^3(sinθ)^1/2-7/8Dz_*^3(sinθ)^5/2](sinθ)^3/2  _2F_1(1/4, 1; 1/2, cos^2θ)
    
    The first term in the particular solutions can then be written as, 
    
        -∫_0^π2L/2z_*sinθ^3/2cos(θ)/√(sin(θ))∫_0^θ f(θ')dθ'=L/z_*[∫_0^θ f(θ')dθ'/sinθ]_θ→ 0^θ→π-L/z_*∫_0^πf(θ)/sinθdθ
    
    By virtue of the boundary conditions and consequently the values C_1 and C_2, the terms diverging as 1/sinθ cancel and we are left with,
    
        LC_2/z_*∫_0^πd/dθ[sin^3/2θ  _2F_1(1/4, 1; 1/2, cos^2θ)]cosθ/sinθdθ-L/z_*∫_0^πf(θ)/sinθdθ
    
    
    The first of the above two integrals then evaluates to -Γ(1/4)^3 Lz_*^2(2C+D)/32 Γ(3/4)√(2)[Consider the expression sin^3/2θ  _2F_1(1/4, 1; 1/2, cos^2θ) and let cos^2θ=z. Then the expression becomes (1-z)^3/4 _2F_1(1/4, 1; 1/2, z). Using Euler's identity, _2F_1(a, b; c, z)=(1-z)^c-a-b _2F_1(c-a,c-b;c,z) we have _2F_1(1/4, -1/2; 1/2, z)=(1-z)^3/4 _2F_1(1/4, 1; 1/2, z)
    Therefore
    d/dθ(sin^3/2θ _2F_1(1/4, 1; 1/2, cos^2θ))= -tanθ{  _2F_1(-1/2,1/4;1/2;cos ^2θ)-(1-cos ^2θ)^-1/4}.]
    
    
    
    and the last f(θ) integral can be calculated easily and is given as -Lz_*^2/10 (15 C+16 D).
    
    
    
    The integral of the second particular solution of 𝒮 gives 2Lz_*^2/10  (5 C+8 D).
    
    
    
    The total change in complexity therefore turns out to be,
    
        Δ C
        
        =1 8π GLz_*^2/2z_0^3(1+β ^2 γ ^2-πΓ(5/4)^2 (1+2 β ^2 γ ^2)/Γ(3/4)^2)
    
     This result is consistent with the 4- dimensional result for the `strip perpendicular to boost' case of <cit.>. This is because of the different notions of parallel and perpendicular used there.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
      §.§.§  Case B: Strip along `x' boost along `y'
    
    To handle this case we will just change the direction of the boost rather than changing the orientation of the strip. Therefore, the embeddings, eq. (<ref>) remains intach while A gets replaced by Ã and so on. Hence the solution for η^1, in this case, is same as that for eq. (<ref>) with C replaced by C̃ and D by D̃ (appendix (<ref>)). The final expression for the change in complexity is obtained just by making this replacement and is given by.
    
    
    
    
        Δ C=1 8π GLz_*^2/2z_0^3(1+2β ^2 γ ^2-πΓ(5/4)^2 (1+β ^2 γ ^2)/ Γ(3/4)^2)
    
    
    This result is consistent with the 4-dimensional result obtained in <cit.>, but now for the `strip parallel to boost' case.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    § OTHER APPLICATIONS
    
    
     §.§ Holographic entanglement entropy
    
    The general expressions found in section <ref> have many applications when it comes to finding identities for the change in `entanglement entropy'. Consider the case of the change in entanglement entropy when the boundary subregion is deformed. It may be a spatial deformation of the subregion or an evolution in time.Some noted works in this direction can be found in <cit.>. The only contribution to the change then comes from the integration over the boundary of the subsytem itself [Note that 𝒮 plays the role of 𝒮 and ∂ A plays the role of ∂𝒮 here. ].
    
    
    
    
    
        δ S=-1/4G∫_∂ Ag(ξ,r̂') √(ζ) d^d-3η,
     
    where η are coordinates on the boundary of the subsytem ∂ A and r̂' is the co-normal of ∂ A in 𝒮. Note the integration over A was absent in the previous expressions because the subsytem was kept fixed. Further, the solution of the homogeneous version of eq.(<ref>) is necessarily not a trivial solution because of non-trivial boundary conditions. We will demonstrate this for the AdS_3 case. Consider the change of enatnglement entropy for a particular deformation of the boundary subregion. Recall, that the solution of the homogeneous Jacobi equation in this case is given by <cit.>,
    
        N^0=C_1e^τ+C_2e^-τ,  N^1=B_1e^τ+B_2e^-τ,
    
    where τ is the parameter along the space-like geodesic. The tangent component N^2 in this case will be non-zero. However it is completely characterised by it's boundary value. While the normal components N^0 and N^1 are to be obtained by solving the homogeneous Jacobi equation with boundary conditions determined by nature of deformation of subregion A. The boundary ∂ A has two disjoint points. Let us denote them as ∂ A_1 and ∂ A_2. It follows that the deformation vector X (say) is given by (assuming there is no deformation into the bulk X^z=0). The deformation vector field (X=A(x,t)∂_t+B(x,t)∂_x) is then related to the solutions of the homogeneous Jacobi equation through the following relations. For the ingoing branch,
    
        g(X,e_0)=A/z_*τ  lim_τ→-∞z_*(τ)N^0=A(-l/2,t)=2z_*c_2,  lim_τ→∞(τ)N^0=A(l/2,t)=2z_*c_1
        g(X,e_2)=tanhτ/z_*τ B  lim_τ→-∞z_*τ/tanhτN^2=B(-l/2,t)=-2z_*b_2,  lim_τ→∞z_*τ/tanhτN^2=B(l/2,t)=2z_*b_1
        g(X,e_1)=B/z_*  N^1=B/z_*
    
     For the 2+1 dimensionl case the expression for the change of entanglement entropy becomes,
    
        δ S   =   1/4Gg(T,N)|_∂ Aδ^2S   =   1/4G(g(∇_NN,T)+1/2∇_Tg(N^⊥,N^⊥))|_∂ A
    
    
    The constants B,C are therefore fixed in terms of the components of the vector deforming the boundary subregion. For a general asymptotically AdS consider it's conformal completion g̃=Ω^2  g. The asymptotic boundary is the given by Ω=0. The normal is then n_a=∇_aΩ. The deformation vector field can then be considered to purely tangential to the Ω=constant hypersurfaces. On such a hypersurface one has the following decomposition,
    
        ∇_NN=𝒟_NN+𝒦(N,N)+1/2Ω^2g(N,N)∇Ω^2,
    
    where 𝒟 and 𝒦 are the intrinsic covariant derivative and extrinsic curvature of a constant Ω hypersurface. The second variation the reduces to,
    
    
        δ^2S   =   1/4G(g(𝒟_NN+𝒦(N,N),T)+1/2Ω^2∇_T(Ω^2g(N^⊥,N^⊥))+g(N^T,N^T)/2Ω^2∇_TΩ^2)|_∂ A
    
    The last two terms has to be understood as limits. For the Poincare patch of pure AdS, Ω=z. The overall variation turns out to be,
    
        Δ S=1/4G(2(b_1-b_2)^2-4(c_1-c_2)^2-2(B^2/z_*^2)|_∂ A)
    
    For a time-like deformation, more specifcally a constant constant time translation of the subregion could be initiated by X=A∂_t. In such a case b_1=b_2=0 and c_1=c_2, which leads to the result Δ S=0. For a increase in the subsytem size by δ l/2 on each side, one has, b_1=δ l/4 z_*, b_2=-δ l/4 z_*. The change in entropy is then given by,
    
        Δ S=δ l^2/4G
    
    While the approach still holds for higher dimensions, unfortunately, it will give divergent results. This because in higher dimensions the cut off dependent term in entanglement entropy itself depends on the size of the subsytem. Hence any infinitesimall change in the size will give divergent changes in entanglement entropy. This can however be cured if the integrals in the above approach is done on a bulk cut-off surface z=ϵ (say).
    
    
    
    § CONCLUSION
    
    From the information theoretic or computational perspective complexity is an important quantity. It needs to be determined what are the resources needed to perform a particular computational task. In a quest for quantum computation, it has emerged as highly relevant, when it comes to determining the states that allow for fast computation. Holographic principles can be effectively used to calculate the complexity of states that have a dual. This has reduced the daunting task of calculating complexity in field theories to the calculation of certain geometric quantities. In spite of this simplification the determination of subregion complexity geometrically amounts to solving a second order elliptic differential equation. Though, this is possible for a few bulk geometries which are highly symmetric, for more general space-times solving such an equation is difficult. The general consensus is then to do it perturbatively around the solutions known for pure AdS.
    
    One of the approaches is to identify a parameter of the bulk space-time and perform an asymptotic expansion of the volume integral in terms of the parameter. The asymptotic limit of the parameter gives the AdS solution. But this is inherently problem specific and  the expansion parameter of-course depends on the bulk. Our approach on the other hand is a variational one that computes terms in order of the strength of a perturbation, which is conveniently taken to be in the Fefferman-Graham gauge. Our approach is in no way advantageous over the previous one in terms of algebraic simplification. However, the advantage is that it provides universal expressions applicable to any space-time that is a perturbation over AdS or for that matter any back-ground. What it does is to cast the problem in a form that might give further insights.
    
    The expression eq. (<ref>) for example refers to a general perturbation P(∂_a,∂_b) and the restriction of the deviation vector N to the boundary ∂𝒮. The restriction of the deviation vector N to the boundary ∂𝒮, η satisfies a linear elliptic equation. There is considerable amount of literature in the theory of differential equations of the elliptic type. We hope that various bounds on the solutions of equations of this type may yield interesting inequalities in the context of complexity and entanglement entropy.
    
    It is important to note what the ingredients required to calculate the first order change are and how it compares to the known approaches. It is clear that one requires the pull back of the perturbation P on Σ and the deviation of 𝒮 along Σ or the initial constant t slice of pure AdS. It might therefore seem that even by restricting oneself to a Σ such calculations can be done. In other words deviations away from the initial constant t slice are redundant and may just consider a co-dimension one Jacbi equation on Σ . But this may not be completely correct. For example it is not known whether the component of the deviation vector along r̂ satisfies a co-dimension one Jacobi equation on Σ itself (and not on the full spacetime ℳ). If the answer is in the negative then the full co-dimension two Jacobi equation on 𝒮 must be solved, as has been done in our case, in order to find the correct component along r̂. Hence completely restricting oneself to Σ may not be a good approximation after all.  
    
    Finally, we believe that the approach taken by us might have important consequences in the context of bulk reconstruction. Specifically, in the absence of a perturbation, the solutions of the Jacobi equation tell us how the minimal surface deforms with the deformation of the boundary subregion. Hence the Jacobi equation tells us how the bulk is traced as one deforms the boundary subregion. This might be useful to find the flow of the metric induced on a given minimal surface, thus reconstructing a region of space-time.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    § ACKNOWLEDGEMENTS
    
    Most of the integrals and solutions of the differential equations were done using Maple <cit.> and Mathematica <cit.>. Tensor calculations were done with the GRTensor package <cit.> on a Maple platform. AG would like to thank Sudipta Sarkar for helpful discussions. RM would like to thank Koushik Ray for helpful discussions. AG is supported by SERB, government of India through the NPDF grant (PDF/2017/000533). 
    
    
    
    
    
    
    
    
    
    § FEFFERMAN GRAHAM EXPANSION
    
    
     §.§ Asymptotic expansion of BTZ
    
    The rotating BTZ metric is given as,
    
        ds^2=-(r^2-r_+^2)(r^2-r_-^2)/r^2 dt^2
        +r^2/(r^2-r_+^2)(r^2-r_-^2) dr^2+r^2(dx-r_+r_-/r^2 dt)^2,
    
    where r_+,r_- are the radii of the outer and inner horizon respectively. An asymptotic expansion of the above metric can be realized by defining coordinate ρ such that dρ/ρ=dr/r√(f(r)). In terms of ρ this metric becomes,
    
    
        ds^2=dρ^2/ρ^2+ρ^2 [(-dt^2+dx^2)+1/ρ^2((r_+^2 + r_-^2)/2dt^2 - 2 r_+ r_- dt dx+(r_+^2 + r_-^2)/2dx^2)
        +1/ρ^4(-(r_-^2-r_+^2)^2/16dt^2+(r_-^2-r_+^2)^2/16dx^2)]
    
    By writing the metric in coordinates (ρ=1/z), one can identify the perturbations as,
    
        (1)h_μν=
        [  (r_+^2 + r_-^2)/2                  0           -r_+ r_-;                  0                  0                  0;           -r_+ r_-                  0 (r_+^2 + r_-^2)/2. ]
    
    The components of C^μ _νρ can be obtained with the expressions for (1)h_μν. Denoting a=(r_+^2 + r_-^2)/2 and b=-r_+ r_- one has,
    
        (1)C ^t  _z  t =-z  a, (1)C^x  _z  t =z  b, (1)C ^t  _t  z =-z  a, (1)C ^x _t  z =z  b(1)C ^t  _x  z =-z  b, (1)C^x  _x  z =z  a, C ^t  _z  x =-z  b, (1)C ^x  _z  x =z  a
    
     
    
     §.§ Asymptotic expansion of Boosted Black brane
    
     The boosted black brane metric in holographic coordinates is of the form 
    
        ds^2=R^2 z^2[-𝒜(z)dt^2+ℬ(z)dx^2+𝒞(z)dtdx+dx^2+dz^2 f(z)],
    
    where,
    
        𝒜(z)=1-γ^2(z z_0)^3,  ℬ(z)=1+β^2γ^2(z z_0)^3,
        𝒞(z)=2βγ^2(z z_0)^3,  f(z)=1-(z z_0)^3
    z_0 is the location of the horizon while 0≤β≤1 is the boost parameter. As usual γ is defined as 1√(1-β^2). The boost here is along x direction. By writing the metric in asymptotic coordinates one can view this as a perturbation over AdS_4. This is achieved by defining a coordinate ρ such that,
    
        dz z√(f(z))=dρρ
    
    Integrating this and setting the integration constant to (ρ_0^3=4z_0^3) we get,
    
        1 z^2=1ρ^2(1+(ρρ_0)^3)^4 3=1ρ^2g(ρ)^4 3
    
    Now we expand the metric coefficient upto second order in (ρρ_0)^3, Substituting this back in the metric we get
    
        ds^2
         =R^2ρ^2[dρ^2+(η_μν+ρ^3γ^(3)_μν+ρ^6γ^(6)_μν)dx^μdx^ν]
    
    Where
    
        γ^(3)_μν=[   -(1 3-γ^2)(1 z_0)^3         βγ^2(1 z_0)^3                     0;         βγ^2(1 z_0)^3 (1 3+β^2γ^2)(1 z_0)^3                     0;                     0                     0          1 3(1 z_0)^3 ]
    
    One can check that Tr( γ^(3)_μν)=0 and
    
        γ^(6)_μν=[    -(2 9+8 3γ^2)116 z_0^6         -1 6βγ^2(1 z_0)^6                         0;         -1 6βγ^2(1 z_0)^6 (2 9-8 3β^2γ^2)1 16 z_0^6                         0;                         0                         0              2 91 16z_0^6 ]
    
    The perturbation (1)P_μν  and (2)P_μν can be read off as, (1)P_μν=γ^(3)_μνz and 1/2(2)P_μν=γ^(6)_μνz^4 respectively. Also, the following notation will be used in the main text C=(1 3+β^2γ^2)1 z_0^3, D=1 31 z_0^3.
    
    
    For boost along the y axis, (1)P(∂_μ,∂_ν) and (2)P(∂_μ,∂_ν) is of the form,
    
    
        (1)P _μν=([  Ã z    0 B̃ z    0;    0 C̃ z    0    0; B̃ z    0 D̃ z    0;    0    0    0    0;  ])    1/2(2)P _μν=([  Ã'  z^4        0 B̃'  z^4        0;        0 C̃'  z^4        0        0; B̃'  z^4        0  D̃' z^4        0;        0        0        0        0;  ])
    
    
     
    where C̃=1 3(1 z_0)^3,  D̃=(1 3+β^2γ^2)(1 z_0)^3,  C̃^'=2 91 16z_0^6,  D̃^'=(2 9-8 3β^2γ^2)1 16 z_0^6,  B=B̃=βγ^2(1 z_0)^3. 
    
    
    
    § SECOND ORDER VOLUME VARIATION
    
     In this section we will briefly go through the derivation of the second variation of the area functional including metric perturbation and deviation of the surface. Starting from the area functional and computing its second variation gives the following expression,
     
        δ^2A=∫ d^nτ (δ√(h))h^ab(g(∇_∂_aN,∂_b)+1 2(1)P(∂_a,∂_b))_I+∫ d^nτ √(h)(δ h^ab)(g(∇_∂_aN,∂_b)+1 2(1)P(∂_a,∂_b))_II
        
        +∫ d^nτ √(h)h^abδ(g(∇_∂_aN,∂_b)+1 2(1)P(∂_a,∂_b))_III.
    
    The above expression now has to be restructured so as to extract the first order inhomogeneous Jacobi equations. In order to achieve this we will first separate this expression into terms with perturbations and the ones without. We will then separate the terms containing the normal component of the deviation vector from those containing the tangent component. At the very end we will isolate the boundary terms from the bulk ones. To carry out this formidable task we mark the individual terms as I, II and III. Term I in the above expression can be written as,
     
        I=∫ d^nτ √(h)h^cdh^ab[g(∇_∂_cN,∂_d)g(∇_∂_aN,∂_b)]_A_1
        +∫ d^nτ √(h)h^cdh^ab[P(∂_a,∂_b)g(∇_∂_cN,∂_d)+1 4h^abh^cdP(∂_c,∂_d)(1)P(∂_a,∂_b)]_A_2,
    
     where we have used the first variation of h_ab. On using the expression h^ach^bd P(∂_a,∂_b)g(∇_∂_cN,∂_d)=P((∇_∂_cN)^T,∂_a), term II can be written as, 
    
        II=-∫ d^nτ √(h)[h^ach^bd(g(∇_∂_cN,∂_d)+g(∇_∂_dN,∂_c))g(∇_∂_aN,∂_b)]_B_1
        -∫ d^nτ √(h)[2h^ab(1)P((∇_∂_aN)^T,∂_b)+1 2h^ach^bd(1)P(∂_a,∂_b)(1)P(∂_c,∂_d)]_B_2
    
    Similarly term-III can be written as,
    
        III=∫ d^nτ √(h)h^ab(g(∇_N∇_∂_aN,∂_b)+g(∇_∂_aN,∇_∂_bN))_C_1
        +∫ d^nτ √(h)h^ab(2(1)P(∇_∂_aN,∂_b)+2g(C(∂_a,N),∂_b)+1 2(2)P(∂_a,∂_b))_C_2
    
     From expression of eq. (<ref>), eq. (<ref>) and eq. (<ref>) we can see that we have been successful in separating terms with perturbation from those without. Now we will try to further simplify these. Let us consider the terms without perturbation first.
     
     
     
     
     
     
     
    
    
     §.§ Terms without perturbations
    
    The terms without perturbation have been grouped as A_1 ,B_1 and C_1. Their sum is therefore total contribution from terms not containing perturbation quantities.
    
        A_1+B_1+C_1= ∫ d^nτ √(h)h^cdh^abg(∇_∂_cN,∂_d)g(∇_∂_aN,∂_b)_(1)
        
         +∫ d^nτ √(h)h^abg(∇_∂_aN,∇_∂_bN)-∫ d^nτ √(h)[h^ach^bd(g(∇_∂_cN,∂_d)+g(∇_∂_dN,∂_c))g(∇_∂_aN,∂_b)]_(2)
        +∫ d^nτ √(h)h^abg(∇_N∇_∂_aN,∂_b)_(3)
    
    For our convenience we have seperated the terms and identified them as (1),(2) and (3). Now let us first consider the term (2)
        h^abg(∇_∂_aN,∇_∂_bN)-[h^ach^bd(g(∇_∂_cN,∂_d)+g(∇_∂_dN,∂_c))g(∇_∂_aN,∂_b)]
        
        =h^abg((∇_∂_aN)^⊥,(∇_∂_bN)^⊥)-h^ach^bdg(D_∂_dN^T,∂_c)g(D_∂_aN^T,∂_b)+2h^abg(D_∂_aN^T,W_N^⊥(∂_b))
        
        -h^abg(W_N^⊥(∂_a),W_N^⊥(∂_b)),
    
    where we have used the decompositions
    ∇_∂_aN^T=D_∂_aN^T+K(∂_a,N^T); ∇_∂_aN^⊥=-W_N^⊥(∂_a)+∇^⊥_∂_aN^⊥. As described before K(∂_a,∂_b) is the extrinsic curvature and W_N^⊥(∂_a) is the shape operator. They are related by the Weingarten equation g(W_N^⊥(∂_a))=g(N^⊥,K(∂_a,∂_b)). The shape operator van also be shown to be self adjoint i.e g(W_N^⊥(∂_a),∂_b)=g(W_N^⊥(∂_b),∂_a). Using this we can write the above expression as
    
        (2)=h^ab(g(K(∂_a,N^T),K(∂_b,N^T))+2g(K(∂_a,N^T),∇^⊥_∂_bN^⊥)+g(∇^⊥_∂_aN^⊥,∇^⊥_∂_bN^⊥)
        +2g(D_∂_aN^T,W_N^⊥(∂_b))-g(W_N^⊥(∂_a),W_N^⊥(∂_b))-h^ach^bdg(D_∂_dN^T,∂_c)g(D_∂_aN^T,∂_b)
    
    In the above expression we have been able to seperate tangent contributions from normal contributions. Next we consider the term (3) in the expression for δ^2_NA
    
        (3)= h^abg(∇_N∇_∂_aN,∂_b)
        =-[Ric(N^T,N^T)+2Ric(N^T,N^⊥)+Ric(N^⊥,N^⊥)]
        +h^ab[g(D_∂_a(∇_N^TN^T)^T,∂_b)-2g(D_∂_aW_N^⊥(N^T),∂_b)
        +g(D_∂_a([N^⊥,N^T])^T,∂_b)+g(D_∂_a(∇_N^⊥N^⊥)^T,∂_b)]
    
    A similar exercise for term (1) yields,
    
        ∫ d^nτ √(h)h^cdh^abg(∇_∂_cN,∂_d)g(∇_∂_aN,∂_b)
        =∫ d^nτ √(h)h^abh^cdg(D_∂_aN^T,∂_b)g(D_∂_cN^T,∂_d)
    
    We see that the expression for (1) has contributions only from tangent terms. Adding (1), (2) and (3) we get the full expression. Note that the equations of motion H=0 and the Gauss and Codazzi identities (<ref>) has been used. The full expression δ^2_NA coming from terms that do not contain perturbation contributions is,
    
    
        δ^2_NA=∫ d^nτ √(h)∇_∂_a[h^abg(N^⊥,∇^⊥_∂_bN^⊥)]-∫ d^nτ √(h)g(ℒ(N^⊥),N^⊥)
         
        +∫ d^nτ √(h)h^cdg(D_∂_d[N^Th^abg(D_∂_aN^T,∂_b)],∂_c)+∫ d^nτ √(h)h^abg(D_∂_a([N^⊥,N^T])^T,∂_b)
        +∫ d^nτ √(h)h^abg(D_∂_a(∇_N^⊥N^⊥)^T,∂_b),
    
    where ℒ(N^⊥) is the Jacobi operator defined in section <ref>.
    Thus we have achieved the decomposition of the terms without perturbation in terms of tangent and normal contributions. Next we consider the terms with perturbations.
    
    
    
     §.§ Terms containing perturbations
    
    The terms containing perturbations have been grouped together. Thus the contribution to the area variation from terms containing perturbations are,
    
        δ^2_PA
        =∫ d^nτ √(h)[2h^abg((1)C(∂_a,N),∂_b)+2h^ab(1)P((∇_∂_aN)^⊥,∂_b)+h^cdh^ab(1)P(∂_c,∂_d)g(∇_∂_aN,∂_b)]
        +∫ d^nτ √(h)[h^ab 2(2)P(∂_a,∂_b)-1 2h^ach^bd(1)P(∂_a,∂_b)(1)P(∂_c,∂_d)+1 4h^abh^cd(1)P(∂_c,∂_d)(1)P(∂_a,∂_b)]
    
    Now our objective is to decompose the above expression in terms of tangent and normal contributions.
    By using the identities in appendix <ref> we can write the above expression as
    
        δ^2_PA
         =∫ d^nτ √(h)h^ab∇_∂_a[h^cd(1)P(∂_c,∂_d)g(N,∂_b)]
        
         +∫ d^nτ √(h)[2h^abg((1)C(∂_a,N^⊥),∂_b)+2h^ab(1)P((∇_∂_aN^⊥)^⊥,∂_b)-h^abh^cd(1)P(∂_c,∂_d)g(N,∇_∂_a∂_b)]
        
         +∫ d^nτ √(h)[h^ab 2(2)P(∂_a,∂_b) -1 2h^ach^bd(1)P(∂_a,∂_b)(1)P(∂_c,∂_d)+1 4h^abh^cd(1)P(∂_c,∂_d)(1)P(∂_a,∂_b)]
    
    Now our task is to decompose the above expression further in terms of bulk and boundary contributions
    
        (∇_∂_a(1)P)(N^⊥,∂_b)=g((1)C(∂_a,N^⊥),∂_b)+g((1)C(∂_a,∂_b),N^⊥)
    
    Substituting this we get the term containing (1)C as
    
        h^ab(∇_∂_a(1)P)(N^⊥,∂_b)+h^ab(1)P((∇_∂_aN^⊥)^⊥,∂_b)=h^ab∇_∂_a[(1)P(N^⊥,∂_b)]+ h^abh^cd(1)P(∂_b,∂_d)g(N^⊥,K(∂_a,∂_c))
        -h^ab(1)P(N^⊥,(∇_∂_a∂_b)^T)
    
    Now we should note that 
    
        ∇_∂_a[√(h)h^ab(1)P(N^⊥,∂_a)]
        =√(h)h^ab∇_∂_a[(1)P(N^⊥,∂_b)]-√(h)h^ab(1)P(N^⊥,(∇_∂_a∂_b)^T)
    
    Which is same as the term (2)P. Similarly,
    
        ∇_∂_a[√(h)h^abh^cd(1)P(∂_c,∂_d)g(N,∂_b)]
         =√(h)h^ab∇_∂_a[h^cd(1)P(∂_c,∂_d)g(N,∂_b)]-√(h)h^abh^cd(1)P(∂_c,∂_d)g(N,∇_∂_a∂_b)
    
    Which is same as the term (1)P. Where we have used the fact that ([∂_a,∂_b]=0). Substituting this back in the expression for δ^2_PA we get,
    
        δ^2_PA=∫ d^nτ ∂_a[√(h)h^abh^cd(1)P(∂_c,∂_d)g(N,∂_b)]+2∫ d^nτ ∂_a[√(h)h^ab(1)P(N^⊥,∂_b)]
        
         +2∫ d^nτ √(h)(h^abh^cd(1)P(∂_b,∂_d)g(N^⊥,K(∂_a,∂_c))-h^abg((1)C(∂_a,∂_b),N^⊥))
        
         +∫ d^nτ √(h)[h^ab 2(2)P(∂_a,∂_b)-1 2h^ach^bd(1)P(∂_a,∂_b)(1)P(∂_c,∂_d)+1 4h^abh^cd(1)P(∂_c,∂_d)(1)P(∂_a,∂_b)]
    
    Thus we have been able to decompose the expression into terms containing normal and tangent contribution and further we have separated the terms contributing to the boundary and bulk. Next we will add the contribution from the terms with and without perturbation to obtain the total variation.
    
    
    
     §.§ Total variation
    
     The expression for the total variation upto second order is obtained by adding eq. (<ref>) and eq.(<ref>). It ontains terms which involve ∇_NN. We would however like to write it in terms of ∇_N^⊥N^⊥. This can be done with the following steps,
    ∇_N N=∇_N^⊥ N^⊥+∇_N^⊥ N^T+∇_N^T N^⊥+∇_N^T N^T
    
    Note that the third surface term can be written as,
    
    h^abg(D_∂_aN^T,∂_b)g(N^T,r̂)   =   γ^ABg(D_∂_AN^T,∂_B)g(N^T,r̂)+g(∇_N^TN^T,r̂)
       =   γ^ABg(N^T,K^(𝒮↪Σ)(∂_A,∂_B))g(N^T,r̂)+g(∇_N^TN^T,r̂)
    
    Denote by K^(𝒮↪ℳ)(∂_A,∂_B) as the extrinsic curvature of 𝒮 in ℳ. Similarly denote the other extrinsic curvature. Then from Gauss's decomposition it follows that,  
    
    K^(𝒮↪ℳ)(∂_A,∂_B)=K^(𝒮↪Σ)(∂_A,∂_B)+K^(Σ↪ℳ)(∂_A,∂_B)
    
    Since 𝒮 is minimal in ℳ we have,
    γ^ABK^(𝒮↪ℳ)(∂_A,∂_B)=γ^ABK^(𝒮↪Σ)(∂_A,∂_B)=γ^ABK^(Σ↪ℳ)(∂_A,∂_B)=0
    
    from which it follows that K^(Σ↪ℳ)(r̂,r̂)=0, since h^abK^(Σ↪ℳ)(∂_a,∂_b)=0. The following relations can be obtained by using this constraint.
    
    g([N^⊥,N^T],r̂)=g(∇_N^⊥N^T,r̂)+g(N^⊥,K^(Σ↪ℳ)(r̂,r̂))g(N^T,r̂)
    
    g(∇_N^⊥N^⊥,r̂)=g(∇_N N-∇_N^⊥N^T-g(∇_N^TN^T,r̂)+g(N^⊥,K^(Σ↪ℳ)(r̂,r̂))g(N^T,r̂)
    
    Finally putting everything back into the total variation one obtains the following expression. This expression is more useful because it has been expressed in terms of quantities which are directly known from boundary conditions or as solutions of the IJEs.
    
        δ^2A=∫ d^nτ √(h)(h^abh^cd(1)P(∂_b,∂_d)g(N^⊥,K(∂_a,∂_c))-h^abg((1)C(∂_a,∂_b),N^⊥))
        
         +∫ d^nτ √(h)[h^ab 2(2)P(∂_a,∂_b)-1 2h^ach^bd(1)P(∂_a,∂_b)(1)P(∂_c,∂_d)+1 4h^abh^cd(1)P(∂_c,∂_d)(1)P(∂_a,∂_b)]
        
         +∫ d^nτ ∂_a[√(h)h^abh^cd(1)P(∂_c,∂_d)g(N,∂_b)]+2∫ d^nτ ∂_a[√(h)h^ab(1)P(N^⊥,∂_b)]
        
         +∫ d^nτ √(h)h^abg(D_∂_a(∇_NN)^T,∂_b)+∫ d^nτ ∇_∂_a[√(h)h^abg(N^⊥,∇^⊥_∂_bN^⊥)]
    
    
    § SECOND ORDER JACOBI EQUATION
    
    In this appendix we will try to give a brief derivation of the second order Jacobi equation. Since the first order Jacobi equation is obtained by taking first variation of the mean curvature (H) and imposing H=0. Similarly the second order Jacobi equation is obtained by taking a second variation of the mean curvature. Now the first variation of the mean curvature is given by
    
        (∇_N+δ_P) H=h^abL_ab^⊥-g(H,(∇_∂_cN)^⊥ )h^cd∂_d+∇_N^TH
    
    Where
    
        L_ab=B_ab-D_ab+S^(N)_ab=∇_∂_a∇_∂_bN+R(N,∂_a)∂_b-∇_(∇_∂_a∂_b)^T N+2 h^cdK(∂_a,∂_c)g(N,K(∂_b,∂_d))
    
    Note that in general [N^⊥, ∂_a]≠ 0. However in this appendix for simplicity we are only considering variation where [N^⊥,∂_a]=0 and assume that N^⊥  is perpendicular to the surface 𝒮. In appendix (<ref>) we will show that the commutator terms actually cancel among each other. One can easily check from eq. (<ref>) that by using (H=0)condition we get back the first order Jacobi equation. Now we need to take a first variation of eq. (<ref>) to obtain the second order Jacobi equation. 
    
        (∇_N+δ_P)(h^abL_ab^⊥-g(H,(∇_∂_cN)^⊥ )h^cd∂_d+∇_N^TH)=(∇_N^⊥+δ_P)(h^abL_ab^⊥)
    
    Where we have used the on shell condition (H=0). Thus our first task is to take the covariant derivative of the quantities B^⊥_ab, D^⊥_ab and S^(N)_ab respectively.
    
    
    
    
      §.§.§ Second order Jacobi equation without perturbation
     We will drop the superscript from N^⊥, to avoid clutter. For convenience let  J:=∇_N N. First we will start with B_ab=∇_∂_a∇_∂_bN+R(N,∂_a)∂_b. Thus
    
        h^ab∇_N B_ab^⊥=h^ab∇_N(B_ab-B_ab^T)=h^ab∇_N(B_ab-g(B_ab,∂_c)h^cd∂_d)
        
        =h^abM_ab^⊥-h^abg((B_ab)^⊥,(∇_∂_cN)^⊥)h^cd∂_d-(∇_G N)^⊥.
    
    Since we can only work with the tangent vectors, the normal part is obtained by subtracting the tangent part from the whole. The rest of the expression above follows from straight forward computation of the covariant derivatives. In the above expression the quantities M_ab and G are given below
    
        M_ab=[4 R(N,∂_a)∇_∂_bN+(∇_∂_a R)(N,∂_b)N+(∇_N R)(N,∂_a)∂_b+R(N,∇_∂_a∂_b)N.
        .+∇_∂_a∇_∂_bJ+R(J,∂_a)∂_b], 
         G=h^abB_ab^T=h^abg(B_ab,∂_c)h^ce∂_e.
    
    In the same manner we can take covariant derivative of D_ab=∇_(∇_∂_a∂_b)^T N
        h^ab∇_N D_ab^⊥=h^ab∇_N(D_ab-D^T_ab)=h^abQ_ab^⊥+h^ab(∇_L^T_abN)^⊥-h^abg(D_ab^⊥,(∇_∂_cN)^⊥)h^cd∂_d-(∇_F N)^⊥
    
    As before this expression is obtained by subtracting the derivative of the tangent part from the derivative of the total. In order to achieve this form the following identity has been used [N,(∇_∂_a∂_b)^T]=L^T_ab. In the above expression the form of Q_ab and F is given below
    
        Q_ab=R(N,(∇_∂_a∂_b)^T)N+∇_(∇_∂_a∂_b)^TN,   
         F=h^abD^T_ab=h^abg(D_ab,∂_c)h^ce∂_e
    
    The last term is the Simon's term S^(N)_ab=-h^cdK(∂_a,∂_c)[g(∇_∂_bN,∂_d)+g(∂_b,∇_∂_dN)] this is a normal vector hence its covariant derivative is obtained directly and the following lengthy expression is obtained
    
        h^ab∇_N S^(N)_ab=2h^abh^cdg(N,K(∂_b,∂_d)) L_ac^⊥-h^abg(S^(N)_ab,(∇_∂_eN)^⊥)h^ef∂_f+h^abS^(J^⊥)_ab
           -2h^abh^cdK(∂_a,∂_c)(g(R(N,∂_b)N,∂_d)+g(∇_∂_bJ^T,∂_d)+g(∇_b N,∇_d N))
    
    Following results from eq.(<ref>), (<ref>) and eq. (<ref>) above we see that
    
        h^ab(∇_NL^⊥_ab)=h^ab(∇_N B_ab^⊥-∇_N D_ab^⊥+∇_N S^N_ab)
        =h^ab(M_ab^⊥-Q_ab^⊥+S^(J^⊥)_ab)-2h^ab(∇_L^T_abN)^⊥
           +2h^abh^cdg(N,K(∂_b,∂_d)) L_ac^⊥-2h^abh^cdK(∂_a,∂_c)(g(R(N,∂_b)N,∂_d)+g(∇_∂_bJ^T,∂_d)+g(∇_b N,∇_d N))
    
    Where we have used the first order Jacobi equation(h^abL^⊥_ab=0) in obtaining the above result. Now one can further simplify this by using eq(<ref>),(<ref>) and the onshell condition (H=0) to obtain the following expression
    
        h^ab(∇_NL^⊥_ab)=h^ab((∇_∂_a∇_∂_bJ)^⊥-(∇_(∇_∂_a∂_b)^T J)^⊥+(R(J,∂_a)∂_b)^⊥+S^J_ab)+h^ab(4 R(N,∂_a)∇_∂_bN.
        .+(∇_∂_a R)(N,∂_b)N+(∇_N R)(N,∂_a)∂_b)^⊥-2h^ab(∇_L^T_abN)^⊥+2h^abh^cdg(N,K(∂_b,∂_d)) L_ac^⊥
        -2h^abh^cdK(∂_a,∂_c)(g(R(N,∂_b)N,∂_d)+g(∇_∂_bJ^T,∂_d)+g(∇_b N,∇_d N))
    
    Finally one can check that 
    
        ∇_N(h^abL_ab^⊥)=h^ab(∇_NL^⊥_ab)+(∇_N h^ab)L^⊥_ab=h^abL^(J)_ab+h^ab(4 R(N,∂_a)∇_∂_bN+(∇_∂_a R)(N,∂_b)N.
        .+(∇_N R)(N,∂_a)∂_b)^⊥-2h^ab(∇_L^T_abN)^⊥+4h^abh^cdg(N,K(∂_b,∂_d)) L_ac^⊥
        -2h^abh^cdK(∂_a,∂_c)(g(R(N,∂_b)N,∂_d)+g(∇_∂_bJ^T,∂_d)+g(∇_b N,∇_d N))
    
    Where we have used eq. (<ref>) and L^(J)_ab is the first order Jacobi equation satisfied by J=∇_NN. This is consistent with any second order equation. One will get eq. (<ref>) as the second order Jacobi equation when there is no perturbation. Next we will include the effect of perturbations and check how eq. (<ref>) changes.
    
    
    
      §.§.§ Second order Jacobi equation with perturbation
    
     When metric perturbations are included our definition of first order equation changes to L_ab=B_ab-D_ab+S^(N)_ab+(1)C(∂_a,∂_b)+h^cd(1)P(K(∂_a,∂_b),∂_c)∂_d. Naturally, the second order deviation equation (eq.(<ref>)) will now include derivatives of perturbations. The new terms will come from derivative of C(∂_a,∂_b)
        h^ab∇_N (1)C(∂_a,∂_b)^⊥=h^abU_ab^⊥-h^abg((1)C(∂_a,∂_b)^⊥,(∇_∂_cN)^⊥)h^cd∂_d-(∇_E N)^⊥,
    
    where E:=h^ab(1)C(∂_a,∂_b)^T and U_ab:=[(∇_N (1)C)(∂_a,∂_b)+ (1)C(∇_∂_aN,∂_b)+(1)C(∇_∂_aN,∂_b)]
    As before the derivative of the normal term is obtained by subtracting the derivative of the tangent part from the derivative of the total term. The covariant derivative of the final term gives,
    
        h^ab∇_N(h^cdK(∂_a,∂_c)(1)P(∂_b,∂_d))=h^abh^cd(1)P(∂_b,∂_d)(L^⊥_ac-g(K(∂_a,∂_c),(∇_m N)^⊥)h^mn∂_n)
        
        +h^cdh^abK(∂_a,∂_c)(2g((1)C(N,∂_b),∂_d)+(1)P(∇_∂_bN,∂_d)+(1)P(∂_b,∇_∂_dN))
    
    Now we will provide the effect of action of δ_p on the terms. First we will check the action on L_ab. We will tabulate the action of δ_p on various terms contained in L_ab such as B_ab, D_ab etc.
    
        h^abδ_P B_ab=(1)C(N,∇_∂_a∂_b)+∇_N (1)C(∂_a,∂_b)=:h^abℳ_ab,  
        h^abδ_P B_ab^T
        =h^abℳ_ab^T +h^abP(B_ab^⊥,∂_c)h^cd∂_d
    
    The variation of the normal component therefore takes the form,
    
        h^abℳ_ab^⊥ -h^abP(B_ab^⊥,∂_c)h^cd∂_d.
    
    Similarly one can obtain the variation of the other terms,
    
        h^abδ_P D_ab=h^ab𝒬_ab^⊥-h^abP(D_ab^⊥,∂_c)h^cd∂_d,
    
    where h^ab𝒬_ab:=∇_E N+h^ab(1)C(N,(∇_∂_a∂_b)^T) and E=P(∇_∂_a∂_b,∂_c)h^cd∂_c+g((1)C(∂_a,∂_b),∂_c)h^cd∂_c-h^ceh^dfg(∇_∂_a∂_b,∂_c)P(∂_e,∂_f)∂_c. The variation of the S_ab term is algebraically a bit more tedious to compute and yields,
    
        h^abδ_P S_ab=-2h^abh^ceh^dfP(∂_e,∂_f)K(∂_a,∂_c)g(N,K(∂_b,∂_d))
        
        +2h^abh^cd[(1)C(∂_a,∂_c)^⊥-h^mnP(K(∂_a,∂_c),∂_m)∂_n]g(N,K(∂_b,∂_d))
        
        -h^abh^cdK(∂_a,∂_c)[2P(∇_bN,∂_d)+2g(C(∂_b,N),∂_d)]
    
    The variation of the (1)C(∂_a,∂_b) term gives
    
        h^abδ_P (1)C(∂_a,∂_b)^⊥=(2)C(∂_a,∂_b)^⊥-2 P((1)C(∂_a,∂_b))^⊥-h^abP((1)C(∂_a,∂_b)^⊥,∂_c)h^cd∂_d.
    
    Note that we have introduced a notation (2)C(∂_a,∂_b) which is nothing but the expression for C(∂_a,∂_b) with (1)P(∂_a,∂_b) replaced by (2)P(∂_a,∂_b). The final term yields,
    
        h^abδ_P(h^cdK(∂_a,∂_c)P(∂_b,∂_d))
        
        =h^ab[-h^ceh^dfP(∂_e,∂_f)K(∂_a,∂_c)P(∂_b,∂_d)+h^cd((1)C(∂_a,∂_c)^⊥-P(K(∂_a,∂_c),∂_m)h^mn∂_n)P(∂_b,∂_d)
        
        +h^cdK(∂_a,∂_c)(2)P(∂_b,∂_d)]
    
    The variation of the h^ab is trivial and gives δ_P h^abL_ab^⊥=-h^aeh^bfP(∂_e,∂_f)L_ab^⊥. Througout this calculation we have used identitites given in appendix (<ref>). 
    
    Further, we have used the fact that [N,∂_a]=0. However if we didn't assume this then several additional tangent contribution would arise. However it will be clearfrom the next appendix  that these tangent terms won't contribute to the final equation. 
    
    
    
    
    § TANGENT VARIATIONS
    
    In the following two sections we will derive certain identities for tangent variations of the extrinsic curvature and mean curvature vector H. These will be used in the main text in section <ref>.
    
    
     §.§ Tangent variation of extrinsic curvature
    
    Let us start by calculationg the tangent variation of the extrinsic curvature of 𝒮. Consider a vector X ∈𝒯_p𝒮 such that [X,∂_a]≠ 0. Note that, 𝒮 being an integral sub manifold [X,∂_a]^⊥=0. After a few algebraic steps it follows that,
    
        ∇_X K(∂_a,∂_b)=(R(X,∂_a)∂_b^⊥+(∇_∂_a∇_∂_bX)^⊥-(∇_∇_(∂_a∂_b)^TX)^⊥+g((∇_∂_a∂_b)^⊥,(∇_∂_cX)^⊥)h^cd∂_d 
        
        +K(∂_a,[∂_b,X])+K([X,∂_a],∂_b))
    
    
     §.§ Tangent variation of H
    
    From the above expression eq. (<ref>), it follows that h^ab[∇_∂_a∇_∂_bX+R(X,∂_a)∂_b-∇_(∇_∂_a∂_b)^TX ]^⊥-2h^aeh^bfg(∇_e X,∂_f)K(∂_a,∂_b)=∇_X H. However there is another elegant way to show it using Gauss's equation.
    
        h^ab[∇_∂_a∇_∂_bX+R(X,∂_a)∂_b-∇_(∇_∂_a∂_b)^TX ]^⊥-2h^aeh^bfg(∇_e X,∂_f)K(∂_a,∂_b)
        
        =h^ab[∇_∂_a^⊥ K(∂_b,X)+K(∂_a,D_∂_bX)
        +∇_X^⊥ K(∂_a,∂_b)-K(∂_a,∇_X∂_b)-K(∇_X∂_a,∂_b)
        
        -∇_∂_a^⊥ K(X,∂_b)+K(X,∇_∂_a∂_b)+K(∇_∂_aX,∂_b)-K((∇_∂_a∂_b)^T,X) ]^⊥-2h^abK(∂_a,∇_∂_b X)=∇_X H
    
    
    § FURTHER IDENTITIES
    
      *  The metric variation of the Riemann tensor can be calculated using the above identities and is given as,
    
        δ_P R(W,U)V=(1)C(W,∇_U V)+∇_W (1)C(U,V)-(1)C(U,∇_W V)-∇_U (1)C(W,V)-(1)C([W,U],V)
    
      *  The metric variation of the extrinsic curvature is given by,
    
    h^abδ_p K(∂_a,∂_b)=(1)C(∂_a,∂_b)^⊥
      *  The following identity will be useful in the calculation of second order perturbation equation.
     
        h^ab[N,(∇_∂_a∂_b)^T]=h^ab[N,g((∇_∂_a∂_b),∂_c)h^cd∂_d]=h^ab∇_N[g((∇_∂_a∂_b),∂_c)h^cd]∂_d
        
        =g(B_ab-∇_(∇_∂_a∂_b)^TN,∂_c)h^cd∂_d
    
      *  The Codazzi Identity is given by
    
        Ric(N^T,N^⊥)=h^abg((∇_N^TK)(∂_a,∂_b),N^⊥)-h^abg(D_∂_aW_N^⊥(N^T),∂_b)
        
         +h^abg(K(∂_a,N^T),∇^⊥_∂_bN^⊥)+h^abg(D_∂_aN^T,W_N^⊥(∂_b))
    
      *  The Gauss equation is given as
    
        -Ric(N^T,N^T)+h^abg(K(∂_a,N^T),K(∂_b,N^T))=-R̅i̅c̅(N^T,N^T)
    
      *  An identity for the change in the affine connection (1)C is,
    
        2h^cdh^abg((1)C(∂_a,∂_c),∂_d)g(N,∂_b)=h^cdg((1)C(∂_c,N^T),∂_d)
    
      *  To extract a surface term from a (1)P dependent term one needs the following identity,
    
        h^cdP(∂_c,∂_d)g(∇_∂_aN,∂_b)=∇_∂_a[h^cdP(∂_c,∂_d)g(N,∂_b)]
        -h^cd[2P((∇_∂_a∂_c)^⊥,∂_d)+2g((1)C(∂_a,∂_c),∂_d)]g(N,∂_b)-h^cdP(∂_c,∂_d)g(N,∇_∂_a∂_b)
    99nielsen2005geometric
      Nielsen,  Michael  A,
      arXiv preprint quant-ph/050207
    
      
    Jefferson:2017sdb 
      R. Jefferson and R. C. Myers,
      
      JHEP  1710, 107 (2017)
      doi:10.1007/JHEP10(2017)107
      [arXiv:1707.08570 [hep-th]].
      
      
    Susskind:2014rva 
      L. Susskind,
      
      [Fortsch. Phys.  64, 24 (2016)]
      Addendum: Fortsch. Phys.  64, 44 (2016)
      doi:10.1002/prop.201500093, 10.1002/prop.201500092
      [arXiv:1403.5695 [hep-th], arXiv:1402.5674 [hep-th]].
      
      
    Brown:2015bva 
      A. R. Brown, D. A. Roberts, L. Susskind, B. Swingle and Y. Zhao,
      
      Phys. Rev. Lett.  116, no. 19, 191301 (2016)
      doi:10.1103/PhysRevLett.116.191301
      [arXiv:1509.07876 [hep-th]].
      
      
    Belin:2018fxe 
      A. Belin, A. Lewkowycz and G. Sárosi,
      
      Phys. Lett. B  789, 71 (2019)
      doi:10.1016/j.physletb.2018.10.071
      [arXiv:1806.10144 [hep-th]].
      
      
    Belin:2018bpg 
      A. Belin, A. Lewkowycz and G. Sárosi,
      
      JHEP  1903, 044 (2019)
      doi:10.1007/JHEP03(2019)044
      [arXiv:1811.03097 [hep-th]].
      
      
    Alishahiha:2015rta 
      M. Alishahiha,
      
      Phys. Rev. D  92, no. 12, 126009 (2015)
      doi:10.1103/PhysRevD.92.126009
      [arXiv:1509.06614 [hep-th]].
      
      bender2013advanced
    C. M. Bender and S. A. Orszag, Advanced mathematical methods for
      scientists and engineers I: Asymptotic methods and perturbation theory.
      Springer Science & Business Media, 2013.
    
    lin1988mathematics
    C.-C. Lin and L. A. Segel, Mathematics applied to deterministic problems
      in the natural sciences, vol. 1. Siam, 1988.
    
      
    Bhattacharya:2012mi 
      J. Bhattacharya, M. Nozaki, T. Takayanagi and T. Ugajin,
      
      Phys. Rev. Lett.  110, no. 9, 091602 (2013)
      doi:10.1103/PhysRevLett.110.091602
      [arXiv:1212.1164 [hep-th]].
      
      
    Mishra:2015cpa 
      R. Mishra and H. Singh,
      
      JHEP  1510, 129 (2015)
      doi:10.1007/JHEP10(2015)129
      [arXiv:1507.03836 [hep-th]].
      
      
    Quijada:2017zif 
      E. Quijada and H. Boschi-Filho,
      
      arXiv:1711.08505 [hep-th].
      
      
    Bhattacharya:2019zkb 
    A. Bhattacharya, K. T. Grosvenor and S. Roy,
    
    arXiv:1905.02220 [hep-th].
    
    Lewkowycz:2018sgn 
      A. Lewkowycz and O. Parrikar,
      
      JHEP  1805, 147 (2018)
      doi:10.1007/JHEP05(2018)147
      [arXiv:1802.10103 [hep-th]].
      
      
    Mosk:2017vsz 
      B. Mosk,
      
      Class. Quant. Grav.  35, no. 4, 045013 (2018)
      doi:10.1088/1361-6382/aaa4e9
      [arXiv:1710.01316 [hep-th]].
      
      
    Bao:2019bib 
      N. Bao, C. Cao, S. Fischetti and C. Keeler,
      
      arXiv:1904.04834 [hep-th].
      
      
    Engelhardt:2019hmr 
      N. Engelhardt and S. Fischetti,
      
      arXiv:1904.08423 [hep-th].
      
    Speranza:2019hkr 
      A. J. Speranza,
      
      arXiv:1904.08012 [gr-qc].
      
      
    Ghosh:2017ygi 
      A. Ghosh and R. Mishra,
      
      Phys. Rev. D  97, no. 8, 086012 (2018)
      doi:10.1103/PhysRevD.97.086012
      [arXiv:1710.02088 [hep-th]].
      
      
    Ryu:2006ef 
      S. Ryu and T. Takayanagi,
      
      JHEP  0608, 045 (2006)
      doi:10.1088/1126-6708/2006/08/045
      [hep-th/0605073].
      
      
    Stanford:2014jda 
      D. Stanford and L. Susskind,
      
      Phys. Rev. D  90, no. 12, 126007 (2014)
      doi:10.1103/PhysRevD.90.126007
      [arXiv:1406.2678 [hep-th]].
      
      
    Alishahiha:2015rta 
      M. Alishahiha,
      
      Phys. Rev. D  92, no. 12, 126009 (2015)
      doi:10.1103/PhysRevD.92.126009
      [arXiv:1509.06614 [hep-th]].
      
      
      
    
    Ghosh:2016fop 
      A. Ghosh and R. Mishra,
      
      Phys. Rev. D  94, no. 12, 126005 (2016)
      doi:10.1103/PhysRevD.94.126005
      [arXiv:1607.01178 [hep-th]].
      
      
    Stewart:1974uz 
      J. M. Stewart and M. Walker,
      
      Proc. Roy. Soc. Lond. A  341, 49 (1974).
      doi:10.1098/rspa.1974.0172
      
      
    Hubeny:2007xt 
      V. E. Hubeny, M. Rangamani and T. Takayanagi,
      
      JHEP  0707, 062 (2007)
      doi:10.1088/1126-6708/2007/07/062
      [arXiv:0705.0016 [hep-th]].
      
      
    Nozaki:2013vta 
      M. Nozaki, T. Numasawa, A. Prudenziati and T. Takayanagi,
      
      Phys. Rev. D  88, no. 2, 026012 (2013)
      doi:10.1103/PhysRevD.88.026012
      [arXiv:1304.7100 [hep-th]].
      
      
    Bhattacharya:2013bna 
      J. Bhattacharya and T. Takayanagi,
      
      JHEP  1310, 219 (2013)
      doi:10.1007/JHEP10(2013)219
      [arXiv:1308.3792 [hep-th]].
      
      
    Wall:2012uf 
      A. C. Wall,
      
      Class. Quant. Grav.  31, no. 22, 225007 (2014)
      doi:10.1088/0264-9381/31/22/225007
      [arXiv:1211.3494 [hep-th]].
      
      
    Fonda:2014cca 
      P. Fonda, L. Giomi, A. Salvio and E. Tonni,
      
      JHEP  1502, 005 (2015)
      doi:10.1007/JHEP02(2015)005
      [arXiv:1411.3608 [hep-th]].
      
      
    Ben-Ami:2016qex 
      O. Ben-Ami and D. Carmi,
      
      JHEP  1611, 129 (2016)
      doi:10.1007/JHEP11(2016)129
      [arXiv:1609.02514 [hep-th]].
      
      
    Karar:2019bwy 
      S. Karar, R. Mishra and S. Gangopadhyay,
      
      arXiv:1904.13090 [hep-th].
      
    Cavini:2019wyb 
      G. Cavini, D. Seminara, J. Sisti and E. Tonni,
      
      arXiv:1907.10030 [hep-th].
      
      
    
    Carmi:2015dla 
      D. Carmi,
      
      JHEP  1512, 043 (2015)
      doi:10.1007/JHEP12(2015)043
      [arXiv:1506.07528 [hep-th]].
      
      
      
    
    Faulkner:2015csl 
      T. Faulkner, R. G. Leigh and O. Parrikar,
      
      JHEP  1604, 088 (2016)
      doi:10.1007/JHEP04(2016)088
      [arXiv:1511.05179 [hep-th]].
      
      
       
       maple
     Maple (13.0). Maplesoft, a division of Waterloo Maple Inc., Waterloo, Ontario.
      
      mathematica
      Wolfram Research, Inc., Mathematica, Version 11.0.1.0, Champaign, IL (2016)
      
      
      
      grtensor
      Peter Musgrave, Denis Pollney and Kayll Lake, 
              GRTensorII software (Version 1.79 (R4), 2001). Queen's University, Kingston, Ontario, Canada.
