Lifshitz holography: The whole shebang

We provide a general algorithm for constructing the holographic dictionary for any asymptotically locally Lifshitz background, with or without hyperscaling violation, and for any values of the dynamical exponents $z$ and $\theta$, as well as the vector hyperscaling violating exponent, that are compatible with the null energy condition. The analysis is carried out for a very general bottom up model of gravity coupled to a massive vector field and a dilaton with arbitrary scalar couplings. The solution of the radial Hamilton-Jacobi equation is obtained recursively in the form of a graded expansion in eigenfunctions of two commuting operators, which are the appropriate generalization of the dilatation operator for non scale invariant and Lorentz violating boundary conditions. The Fefferman-Graham expansions, the sources and 1-point functions of the dual operators, the Ward identities, as well as the local counterterms required for holographic renormalization all follow from this asymptotic solution of the radial Hamilton-Jacobi equation. We also find a family of exact backgrounds with $z>1$ and $\theta>0$ corresponding to a marginal deformation shifting the vector hyperscaling violating parameter and we present an example where the conformal anomaly contains the only $z=2$ conformal invariant in $d=2$ with four spatial derivatives.


Introduction
The use of holographic techniques in order to gain insight into the strongly coupled dynamics of condensed matter systems has attracted considerable interest in the last few years. Gravity duals to quantum critical points exhibiting Lifshitz [3][4][5] or Schrödinger [6,7] symmetry have been put forward and studied extensively. More recently, scaling geometries where translations in the radial coordinate is not an isometry but only a conformal isometry have been proposed as gravity duals to non-relativistic systems exhibiting hyperscaling violation [1,[8][9][10][11][12][13][14]. Hyperscaling violating Lifshitz (hvLf) geometries are characterized by two dynamical exponents, the Lorentz violating exponent z and the hyperscaling violating parameter θ, and take the form where d is the number of spatial dimensions, a = 1, . . . , d, and is the Lifshitz radius. This metric is invariant under time and spatial translations, as well as spatial rotations, but under the anisotropic scaling transformation x a → λx a , t → λ z t, u → λu, (1.2) it transforms homogeneously according to Hence, (1.2) is only a conformal isometry of (1.1) unless θ = 0, which corresponds to the scale invariant Lifshitz (Lif) geometry. For z = 1 the metric (1.1) coincides with the (non-compact part of the) near horizon geometry of relativistic Dp branes [15][16][17][18][19], with the hyperscaling violating exponent θ given by This special case not only provides insight into the physics described by hyperscaling violating backgrounds, but also is an important guide in developing the holographic dictionary for such backgrounds.
As for Dp branes, the holographic relation between the energy scale of the dual field theory and the radial coordinate u can be unambiguously identified through a supergravity probe calculation [20,21]. This determines that the ultraviolet (UV) of the dual theory is located at u = 0, independently of the value of θ, in agreement with the relativistic case z = 1 [15][16][17][18][19]. It follows that the proper identification of the boundary of the geometry (1.1) through a conformal compactification requires a Weyl transformation to the "dual frame" [18,22], where the metric becomes Lifshitz, thus providing an unambiguous definition of the boundary. In the conformal case, θ = 0, such a potential ambiguity does not arise since no field redefinition (including Weyl frame transformations) change the asymptotic behavior of the metric. Given that the curvature invariants scale with u as one might be tempted to conclude that e.g. for θ > 0 there is a curvature singularity as we approach the UV at u = 0. However, given that geometries of the form (1.1) with θ = 0 generically require the presence of a linear dilaton that tends to ±∞ as u → 0, such statements are not well defined since we can tune the curvature singularity at will by changing Weyl frame. In particular, in the dual frame the curvature singularity is completely absorbed in the dilaton. Since this is the proper holographic frame in the case θ = 0, there are no restrictions on θ imposed by requiring absence of curvature singularities in the UV. In the IR one can apply the criterion of [23], which again provides an unambiguous statement about curvature singularities in the presence of scalars.
Restrictions on θ and z do arise, however, from the null energy condition (NEC) where k µ is an arbitrary null vector field, i.e. k µ k µ = 0. The NEC leads to the two constraints Including the relativistic case, z = 1, the solutions of these constraints are: For θ = 0 all cases except I and II admit solutions, which leads to the condition z ≥ 1. A comparison with the relativistic case is instructive. From (1.4) follows that for p ≤ 4 we have θ ≤ 0, corresponding to case IIIa. For p = 5 (1.4) is ill defined but it can be understood as the limit θ → −∞ or θ → +∞, corresponding respectively to cases IIIa and II. Finally, p = 6 gives θ = 9 > d + z = 7 and so it belongs to case II. It is well known that there are no well defined Fefferman-Graham asymptotic expansions in the case of D6 branes [18], which reflects the fact that there is no decoupling limit [15]. A general criterion for the existence of well defined asymptotic expansions is the volume divergence of the on-shell action. For the metric (1.1) in the Einstein frame we get S ∼ˆd u u d+z+1−θ , (1.9) which diverges as u → 0 provided θ ≤ d + z. (1.10) This criterion is independent of the choice of Weyl frame. It follows that all cases except I and II admit well defined asymptotic expansions. Asymptotic expansions, therefore, exist for z > 1, but not for z < 1, and so we will mostly focus on the case z > 1 in the following.
For an extensive list of references on non-relativistic backgrounds, their hyperscaling violating versions and possible string theory embeddings we refer the reader to the following recent papers and references therein [1,[24][25][26]. The body of literature most relevant to us here, however, concerns earlier work on holographic renormalization and the holographic dictionary for asymptotically Lifshitz spacetimes [27][28][29][30][31][32][33][34]. These papers focus mainly on the Einstein-Proca theory, i.e. gravity coupled to a massive vector field, mostly without any scalars and only with conformal (Lifshitz) boundary conditions. Moreover, the emphasis is often on the physically interesting but rather special case d = z = 2. Our aim here is to extend these analyses to the case of general hvLf boundary conditions.
Besides the aforementioned studies on the first principles construction of the holographic dictionary for asymptotically Lifshitz backgrounds of the Einstein-Proca theory, there are few examples where the non-relativistic dictionary has been inferred from a related relativistic dictionary for asymptotically AdS backgrounds. In [35] a 4-dimensional model that admits z = 2 Lifshitz backgrounds was obtained by a dimensional reduction of an axion-dilaton system in 5 dimensions that can be embedded in Type IIB supergravity. In particular, the z = 2 Lifshitz backgrounds are obtained from the reduction of 5-dimensional Schrödinger solutions of the axion-dilaton theory with z = 0, which are asymptotically AdS 5 . This connection was utilized in [36] in order to deduce the holographic dictionary for the Lifshitz backgrounds from the dictionary for asymptotically locally AdS solutions of the axion-dilaton theory developed in [37]. The same model was revisited in [38,39] using the vielbein formalism and a connection between the structure of the sources and Newton-Cartan geometry on the boundary was proposed. Another way to relate the Lifshitz and AdS boundary conditions is a scaling limit where z → ∞. The resulting asymptotic geometry is AdS 2 ×R d−1 . This limit, however, is not very useful in practice because the holographic dictionary for the limiting spacetime is not fully understood -due to the non-compact R d−1 directions and the well-known subtleties associated with AdS 2 holography. Finally, one can study Lifshitz backgrounds with dynamical exponent infinitesimally close to the relativistic value, i.e. z = 1 + , where is small [40,41]. This corresponds to deforming the relativistic CFT with an irrelevant operator and so the analysis must be done with a UV cut-off.
The main goal of the present paper is a systematic derivation of the holographic dictionary for general asymptotically Lif and hvLf backgrounds, for generic values of the dynamical exponents z and θ. In particular, the aim here is not a detailed discussion of the physics of a specific model, but rather the construction of a general algorithm from which the physics can be systematically extracted for any model that admits Lif and hvLf backgrounds. Moreover, throughout this paper we adopt the point of view that the field theory exhibiting Lifshitz or hyperscaling violating Lifshitz symmetry is at the UV -not in the IR -since the physics of Lif or hvLf geometries in the IR can be simply extracted by studying the corresponding UV theory. The IR physics of a geometry which, for example, starts as AdS in the UV and runs to hvLf in the IR (or at some intermediate energy scale) can be studied using standard well known tools for asymptotically locally AdS holography. There is no need for new machinery in that case. Here we are therefore concerned exclusively with backgrounds which are asymptotically locally Lif or hvLf in the UV. For θ > d + z such backgrounds will generically require a different UV completion, but we will not be concerned with this case here.
Our algorithm for constructing the holographic dictionary hinges upon a certain asymptotic solution of the radial Hamilton-Jacobi (HJ) equation [42][43][44][45], subject to asymptotically Lif or hvLf boundary conditions. This asymptotic solution of the radial HJ equation not only provides the necessary local boundary counterterms to render the on-shell action finite, but also is required in order for the variational problem to be well defined both for asymptotically locally AdS [46] and asymptotically non AdS [45] backgrounds. Moreover, the procedure of holographic renormalization based on such an asymptotic solution of the HJ equation is completely equivalent to the traditional method based on asymptotic solutions of the equations of motion [47][48][49]. However, there are two crucial differences between our use of the radial HJ equation and the way it is used in most of the literature. Firstly, we do not need to make an ansatz for the solution of the HJ solution. Finding the correct ansatz becomes increasingly difficult in the presence of matter fields and especially when non AdS boundary conditions are imposed. Moreover, the number of equations obtained by inserting an ansatz into the HJ equation is in general greater than the number of unknown parameters of the ansatz and so the system is overdetermined. Instead, the way we solve the HJ equation is by setting up a recursion procedure based on the covariant expansion of the HJ solution in eigenfunctions of a suitable operator. For scale invariant boundary conditions this operator is usually the relativistic [44] or non-relativistic [29,33] dilatation operator. For more general boundary conditions, such as non-conformal branes or hvLf backgrounds, a generalized dilatation operator is required, such as the one discussed in [37] for relativistic non scale invariant boundary conditions. One of the main results of the present paper is the identification of a suitable set of commuting operators that lead to a recursive solution of the HJ equation with Lif or hvLf boundary conditions [2]. A second point where our approach differs from other approaches to the holographic dictionary is that at no point do we use the general second order equations of motion. In particular, the asymptotic Fefferman-Graham expansions are obtained by integrating the first order flow equations corresponding to the asymptotic solution of the HJ equation. In this way there is no need for making an ansatz for the asymptotic solutions of the equations of motion -the asymptotic form is determined algorithmically by integrating order by order the flow equations. This is particularly useful in the case of non AdS boundary conditions where the form of the asymptotic expansions is a priori unknown and may even contain multiple scales [37].
The paper is organized as follows. In Section 2 we present a general bottom up model that admits both Lif and hvLf backgrounds and we formulate its dynamics in the radial Hamiltonian formalism, which we use later in order to develop the holographic dictionary. Section 3 concerns exclusively homogeneous but anisotropic background solutions of the model presented in Section 2. Both Lif and hvLf backgrounds are discussed in detail and the holographic dictionary for the minisuperspace of homogeneous asymptotically Lif and hvLf backgrounds is obtained. This serves as a self contained warm up for the derivation of the general dictionary for asymptotically locally Lif and hvLf backgrounds that will follow, but also it provides a general description of anisotropic holographic renormalization group (RG) flows. In Section 4 we discuss the boundary conditions corresponding to asymptotically locally Lif and hvLf backgrounds and we present a general algorithm for solving the radial HJ equation iteratively for such backgrounds. This is achieved by covariantly expanding the solution of the HJ equation in simultaneous eigenfunctions of two commuting operators, which as we show are the appropriate generalization of the dilatation operator for anisotropic and non scale invariant boundary conditions. The full holographic dictionary, i.e. the Fefferman-Graham asymptotic expansions, the identification of the sources and 1-point functions of the dual operators, the holographic Ward identities and the conformal anomalies, as well as the covariant boundary counterterms that render the on-shell action finite all follow directly from general asymptotic solution of the HJ equation as is discussed in Section 5. Finally, a number of examples are worked out in Section 6, and a few technical results are presented in the appendices.

The model and radial Hamiltonian formalism
The minimal field content that supports Lifshitz solutions is a massive vector field, or a massless vector field and a scalar, coupled to Einstein-Hilbert gravity. A more general model that includes both these cases and supports in addition hyperscaling violating solutions is the action introduced in [1], namely where κ 2 = 8πG d+2 is the gravitational constant in d + 2 dimensions and S GH is the Gibbons-Hawking term The functions Z(φ), W (φ) and V (φ) are arbitrary, subject only to the condition that the equations of motion admit the desired asymptotic solutions. We will derive these conditions in detail in the subsequent analysis. Moreover, the parameter α > 0 can be removed by a rescaling of the scalar field, but we keep it to facilitate direct comparison with the existing literature, where different conventions are used. Finally, we do not include Chern-Simons terms here in order to keep the spacetime dimension arbitrary throughout most of our analysis. Such terms can be incorporated in the analysis though, once a choice of spacetime dimension has been made. We want to generalize the action (2.1) in two crucial ways, however. Firstly, in order to consistently describe this theory in a Hamiltonian language we need to maintain the U (1) gauge invariance in the presence of a mass term for the vector field. This can be done straightforwardly by introducing a Stückelberg field ω and replacing so that the U (1) gauge transformation leaves B µ invariant. As it turns out, the preservation of the U (1) gauge invariance has important implications for the holographic dictionary. Secondly, in order to be able to develop the holographic dictionary for asymptotically Lifshitz and hyperscaling violating Lifshitz backgrounds simultaneously, it is necessary to go to a generic Weyl frame by means of the Weyl transformation g → e 2ξφ g, (2.5) of the action (2.1), with ξ an arbitrary parameter. As we shall see later, ξ is related to the hyperscaling violation exponent θ in the Einstein frame. With these generalizations, the model we will study is defined by the action where and the Gibbons-Hawking term now takes the form The equations of motion following from this action are (2.9) We will not need these equations in the subsequent analysis, except for demonstrating that the first order equations we will derive for background homogeneous solutions solve these equations.

Radial Hamiltonian formalism
The starting point for the derivation of the holographic dictionary for the action (2.6) is a radial Hamiltonian description of the dynamics, where the radial coordinate is interpreted as the Hamiltonian 'time'. We start by the standard ADM decomposition of the metric [50] as 10) where N and N i are respectively the shift and lapse functions, and γ ij is the induced metric on the radial slices Σ r . In terms of these variables the action (2.6) can be written as a radial integral over the Lagrangian where the extrinsic curvature K ij is given by 12) and D i denotes the covariant derivative with respect to the induced metric γ ij . Moreover, we will use the notation K = γ ij K ij to denote the trace of the extrinsic curvature. Since no radial derivatives of N , N i or A r appear in this Lagrangian, the corresponding canonical momenta vanish identically and these fields play the role of Lagrange multipliers, imposing the usual first class constraints which we will derive shortly. The canonical momenta for the rest of the fields are These relations can be inverted to obtain the generalized velocities in terms of the canonical momentȧ 14) The Hamiltonian is then obtained as the Legendre transform of the Lagrangian, namely where the local densities H, H i and F are given by (2.16) These three quantities appear in the Hamiltonian as coefficients of the three Lagrange multipliers N , N i , and A r respectively, and so the corresponding Hamilton equations yield the three constraints (2.17) These first class constraints reflect the full diffeomorphism and U (1) gauge invariance of the action (2.6).
In particular, this would not have been the case had we not used the Stückelberg mechanism to preserve the U (1) symmetry in the presence of a mass for the vector field. This plays a critical role in our construction of the holographic dictionary. The constraints (2.17) are the basis of the radial Hamilton-Jacobi formulation of the model (2.6). The key new ingredient provided by the Hamilton-Jacobi formalism is the alternative expression for the canonical momenta as gradients of a functional S[γ, A, φ, ω] of the induced fields, namely Inserting these expressions for the momenta in the constraints (2.16) leads to a set of functional partial differential equations for S[γ, A, φ, ω], which is often known as Hamilton's principal function. A fundamental property of the Hamilton-Jacobi approach to the dynamical problem is that the Hamilton-Jacobi equations, i.e. the constraints (2.17), together with the relations (2.18) expressing the momenta as gradients of a 'potential' S[γ, A, φ, ω], provide a full description of the dynamics. In particular, there is no need to consider the second order equations of motion (2.9). By constructing suitable solutions of the Hamilton-Jacobi equations, therefore, we can provide a complete description of the classical dynamical problem, and hence of the holographic dictionary. Our main objective in the subsequent analysis will therefore be to develop a systematic algorithm for solving the Hamilton-Jacobi equations (2.17), subject to the desired boundary conditions. In fact, we only need to focus on the Hamiltonian constraint H = 0, as the other two can be satisfied by construction. In particular, the momentum constraint H i = 0 simply requires the functional S to be invariant with respect to diffeomorphisms on the radial slices Σ r , while the constraint F = 0 imposes U (1) invariance, i.e. it simply requires that S depends on A i and ω only through the gauge-invariant filed B i . Provided then we look for Diff Σr -invariant solutions S[γ, B, φ], the only equation we need to solve is the Hamiltonian constraint H = 0. Of course, the other two constraints will also play a crucial role in the construction of the holographic dictionary, giving rise to certain Ward identities.
Given a solution S[γ, B, φ] of the Hamilton-Jacobi equations, the radial trajectories of the induced fields can be obtained by integrating the first order equations (2.14), where the canonical momenta are expressed as gradients of the given solution of the Hamilton-Jacobi equations as in (2.18). With the gauge choice which we will adopt from now on, these first order equations take the forṁ (2.20) We will use these first order equations in two different but complementary ways. Firstly, making an ansatz for a class of background solutions, these first order equations become analogous to first order BPS equations, while Hamilton's principal function S plays the role of a fake superpotential [51]. We will discuss this in detail in Section 3. The second major application of these equations will be to obtain the asymptotic Fefferman-Graham expansions of the fields, and as a result the holographic dictionary, from the general asymptotic solution of the Hamilton-Jacobi equation subject to specified boundary conditions. The systematic construction of this general asymptotic solution of the Hamilton-Jacobi equation is the subject of Section 4. As we shall see, the general asymptotic solution contains a number of undetermined integration functions. In the Hamilton-Jacobi language these are the 'initial' momenta contained in a complete integral of the Hamilton-Jacobi equation, while in the holographic context they correspond to the renormalized momenta. Via the flow equations (2.20) these undetermined functions give rise to the normalizable modes in the Fefferman-Graham expansions of the fields. The non-normalizable modes, on the other hand, appear as the integration functions of the first order flow equations themselves. The Hamilton-Jacobi formalism, therefore, provides a natural qualitative division of the asymptotic data into two classes, data arising from the integration of the Hamilton-Jacobi equation and data arising from the integration of the first order flow equations. This division in most cases coincides with the separation of the asymptotic data into sources and 1-point functions in the holographic context, but there are exceptions to this rule. An obvious exception is the case of scalars or vector fields with two normalizable modes. More generally, the symplectic form on the space of asymptotic solutions, parameterized by the modes arising from the integration of the Hamilton-Jacobi equation and the first order flow equations, will not be diagonal. The way to identify the sources and 1-point functions out of these asymptotic data in such cases is to diagonalize the symplectic form [45].

Holography for homogeneous anisotropic backgrounds
As a prelude to the general analysis of asymptotically locally Lif and hvLf backgrounds, and in order to outline several of the key steps of our method, it is very instructive to start by discussing the Hamiltonian formalism and the holographic dictionary within the minisuperspace of homogeneous, yet anisotropic, background solutions of the equations of motion.
In particular, in this section we will consider solutions described by the ansatz where a, b = 1, . . . , d. Inserting this ansatz in the equations of motion (2.9) gives the set of equations These equations, except the first and the last, are the equations of motion following from the effective point particle Lagrangian The values of these conserved quantities are zero in the gravitational context, which can be derived by keeping the Lagrange multipliers N and A r in the effective point particle Lagrangian. The canonical momenta following from the Lagrangian (3.3) are and the corresponding Hamiltonian is This Hamiltonian is conserved, but invariance under radial reparameterizations -which would be manifest in (3.3) had we not gauge-fixed the einbein -requires that it is in fact zero. The Hamilton-Jacobi equation therefore is with the canonical momenta expressed as gradients of a function S ef f (f, h, a, φ, ω) of the generalized coordinates so that (3.7) becomes a partial differential equation (PDE) for the function S ef f (f, h, a, φ, ω).

Hamiltonian algorithm for the holographic dictionary
The full holographic dictionary for the backgrounds (3.1) can be constructed from suitable solutions S ef f (f, h, a, φ, ω) of the HJ equation (3.7), without ever using the second order equations (3.2). To this end it is very important to understand the relation between solutions of the HJ equation and solutions of the equations of motion. In particular, the most general solution of the equations of motion can be obtained from a complete integral of the Hamilton-Jacobi equation, i.e. a solution S ef f (f, h, a, φ, ω; π f , π h , π a , π φ , π ω ) that contains as many integration constants as generalized coordinates. These integration constants will eventually be identified with the renormalized momenta, i.e. the renormalized 1-point functions [45]. Such a complete integral is clearly not the most general solution of the HJ equation, but it is all that is needed in order to describe the general solution of the equations of motion. However, the solutions of the HJ equation generically contain branch cuts in field space, and so a given complete integral may not cover the entire solution space, but rather a subset. A discrete set of complete integrals is sufficient to cover the entire space of solutions of the second order equations of motion. There are two types of solutions of the HJ equations we will need: • Exact solutions of the HJ equation These are special but exact solutions of the HJ equations that can be understood as 'fake superpotentials' [51]. Typically they are obtained by finding suitable ansätze that render the HJ equation tractable. Moreover, any discrete branch of the HJ equation is acceptable. 1 The corresponding exact backgrounds that solve the equations of motion are obtained by integrating the flow equations (2.20). Such solutions may or may not contain any integration parameters and they are generically interpreted as RG flows of the dual theory.
• An asymptotic complete integral of the HJ equation This type of solution is the main tool in the construction of the holographic map. It is only required to be an asymptotic solution of the HJ equation, in the sense explained in Fig. 1, but must contain all integration constants required of a complete integral. In order to include these integration constants the asymptotic solution must be obtained up to and including the finite terms in S ef f (f, h, a, φ, ω; π f , π h , π a , π φ , π ω ). These finite terms are exactly the terms that are not completely determined in the asymptotic solution and so are parameterized in terms of a number of undetermined integration constants. Moreover, the condition that the solution must be valid in the asymptotic region A in configuration space requires that a particular branch of the Hamilton-Jacobi solution be chosen. In the Poincaré domain wall example this is the well known fact that only a superpotential with a quadratic term that corresponds to a deformation can be used to construct the holographic dictionary [52]. Constructing such an asymptotic complete integral and deriving the holographic map for asymptotically Lifshitz and hvLf backgrounds is the main purpose of this paper. We now describe this construction within the minisuperspace (3.1) of homogeneous backgrounds, postponing the general case for Section 4.

Asymptotic complete integral and the Fefferman-Graham expansions
Although we are focusing on homogeneous solutions for now, the asymptotic complete integral we want to construct must still correspond to the zero-derivative asymptotic solution of the HJ equation in the full theory, even when the fields have arbitrary spacetime dependent sources. Since for a renormalizable holographic dual the divergent part of the on-shell action must be local in these sources, as well as diffeomorphism and gauge invariant, it follows that the most general form of the divergent part of the HJ solution in the full theory must be of the form for some 'superpotential' U . This restriction, however, does not apply to the finite part of the asymptotic complete integral, for which there is no requirement of locality. This observation is crucial in order to obtain the full complete integral with the correct number of integration constants, which clearly cannot be obtained from the superpotential U that contains up to two integration constants. However, once the divergent part is determined, the finite part can be obtained in terms of a number of undetermined integration constants, as we will show shortly. The form (3.8) of the divergent part of the general asymptotic HJ solution implies that the divergent part of the complete integral S ef f we are interested in for the homogeneous backgrounds takes the form Defining X := φ, Y := −e −2f a 2 , and inserting this point particle HJ function in the Hamiltonian leads to the following PDE for the superpotential U (X, Y ): where the subscripts X and Y denote partial derivatives w.r.t. the corresponding variable. The superpotential equation (3.10) significantly simplifies the problem of determining the divergent part of the general asymptotic complete integral, since we have to solve a PDE in only two variables, but can also be used to obtain exact solutions. Identifying the canonical momenta (3.5) with the gradients of (3.7) and the ansatz (3.9) leads to the first order flow equationṡ Given any solution of the superpotential equation (3.10), asymptotic or exact, the flow equations(3.12) can be integrated to obtain the trajectories of X and Y . Inserting those in turn in (3.11), f , h and a can be determined as well. As we stressed earlier, solutions obtained in this way automatically satisfy the second order equations of motion (3.2).
A last point we must address is the finite part of the asymptotic complete integral, which as we explained cannot be assumed to be of the form (3.9). To this end let us consider a solution S o of the HJ equation, which without loss of generality can be taken to be of the form (3.9). We then seek to determine the possible infinitesimal deformations of this solution, which should give us the full set of integration constants that parameterize a complete integral. Inserting (3.13) in (3.7) and keeping terms up to linear order in δS gives the linear PDE Comparing this with the flow equations (3.11) and (3.12) we see that this equation can be written in the form 15) which shows that only the finite part of the solution S o can be deformed. To determine the complete set of deformations it suffices to consider this equation in the leading asymptotic limit as r → ∞ so that the radial derivative is replaced by the dilatation operator δ D [44]: The characteristic surfaces of this linear first order PDE determine the deformation parameters of the solution S o , which correspond to the full set of normalizable modes. Various solutions of the superpotential equation (3.10) will be discussed in detail in Section 3.6, including the derivation of the general asymptotic complete integral for Lif and hvLf backgrounds.

Lif solutions
In order for the equations (3.2) to admit Lifshitz solutions, the potentials in the action (2.6) must be of the form at least asymptotically, where the various constants are constrained in a way we will specify momentarily. In this section we will assume that this is the exact form of the potentials, but more general potentials will be considered later on. The Lifshitz solutions take the form where the various parameters are related as follows: Note that a possible additive constant in the scalar field has been absorbed in the Lifshitz radius , which we set to 1. These solutions are related in the Einstein frame to the hvLf solutions of [1]. We will discuss the connection of these solutions to hvLf solutions shortly. Moreover, various limits of these solutions deserve special attention.

Special limits
This case is interesting because it corresponds to a massless U (1) gauge field, and so the Action becomes the Einstein-Maxwell-Dilaton (EMD) action. The values of the parameters in this case simplify as follows: In the Einstein frame this case corresponds to hvLf solutions with which are compatible with the NEC solutions III-V provided also θ < d + z. Setting ξ = 0 in these solutions we recover the anisotropic solutions obtained in [5]. Note that necessarily µ = 0 in this case, and so a running scalar is required to support these solutions. The limiting case θ = d + z leads to Q = 0 and was discussed in [24]. However there are more solutions with Q = 0 which we discuss now.
ii) W o = 0, Q = 0: This case also corresponds to a massless U (1) gauge field but now the gauge field in not switched on in the background. The values of the parameters in this case are: As we shall see, these solutions in the Einstein frame are hvLf solutions with These solutions include the zero vector field solution with θ = d + z discussed in [24], but the fact that any θ in the range (3.23) leads to a solution with W o = 0 and Q = 0 was missed in [24] because only the case = d + z + dµξ was considered there.
iii) µ = 0: This is another important special case, where non-relativistic conformal invariance is recovered at least asymptotically. The parameters of the solution now take the simpler form The scalar can be set identically to zero in this case, so that the action (2.6) reduces to Einstein-Proca theory [4]. The scalar is not identically zero necessarily in this case, however, and so it is important to keep ξ as a parameter. Firstly, when we generalize these solutions to inhomogeneous solutions with dependence on the transverse coordinates we will see that there can be non-zero subleading terms in the scalar. Moreover, if the potentials (3.17) are suitably modified at subleading orders, then the scalar can acquire not trivial radial dependence. Both cases of constant scalar and and non-constant scalar with µ = 0 will be studied in detail in Section 6. iv) Dp branes in the dual frame: Finally, it is useful as a reference to obtain the relativistic Dp brane solutions by setting z = 1 in (3.18). The resulting family of solutions with parameters , (3.25) corresponds to Dp branes in the dual frame [15,16].

hvLf solutions
By means of the coordinate transformation and a suitable rescaling of the time and spatial coordinates, the hvLf metric (1.1) takes the form where (3.28) Note that in this coordinate system the UV is located at r → ∞ for θ < 0 and at r = 0 for θ > 0. Inserting this ansatz in the equations of motion (3.2), together with the homogeneous ansatz for the rest of the fields, we find that such solutions exist provided As for the Lifshitz solutions, the additive constant in the scalar field has been absorbed into the Lifshitz radius, which we set to 1. Note that these solutions do not exist for µ = 0, and so they always require a running dilaton. Moreover, the parameter ξ in these solutions is somewhat redundant as we can always set it to zero by a redefinition of θ. For d = 2 and ξ = 0 they reduce to the solutions discussed in Section 3.2.2 of [1]. Note in particular that the independent metric and gauge field hyperscaling violating parameters discussed in [1] are related to our parameters θ and µ respectively.

Special limits
As for the Lifshitz solutions, there are two cases with massless vector. Namely Q = 0 and Q = 0. In the former case the hvLf solutions of the EMD model satisfy the following conditions: These solutions are related to the finite charge density solutions in [53]. Note that, as for the Lifshitz solutions, there is a limiting case of this class of solutions that has Q = 0 and ν z + dν 1 − 1 + dξµ = 0. For ξ = 0 this is the corresponding Lifshitz solution we discussed above but now in the Einstein frame, and it is also the Q = 0 solution discussed in [24]. However, as in the Lifshitz case, there are more solutions with Q = 0.
ii) W o = 0, Q = 0: The class of hvLf solutions with Q = 0 corresponds to the parameter space Setting ξ = 0 in these solutions we reproduce the Einstein frame version of the Lifshitz solutions (3.22) with θ in the range (3.23).
iii) Dp branes in the Einstein frame: Finally, from the relativistic limit z = 1 of the hvLf solutions (3.30) we recover the Einstein frame version of the Dp brane solutions with parameters

Weyl transforming hvLf solutions to Lif solutions
As we have already mentioned, hvLf and Lif solutions are conformally related. This is immediately obvious from the metric (1.1), but it is useful to see how all the parameters of the solutions transform under the relevant conformal transformation, and in particular to clarify the role of the Weyl frame parameter ξ.
Starting with the hvLf (3.27) metric and introducing the new coordinates we obtain while the scalar is given by Note that the UV is located atr → ∞ for all values of θ = 0. It follows that the hvLf metric (3.27) can be written as where µ L = −θµ h /d and g L is a Lifshitz metric with radius L = |θ| h /d. We now observe that if a metric g o solves the equations of motion (2.9) with ξ = 0, then g = e −2ξφ g o solves the equations of motion with non-zero ξ. In particular, let g L = e −2ξ L φ g o be a Lifshitz metric and g h = e −2ξ h φ g o a hvLf one with hyperscaling violating parameter θ that solve the equations of motion corresponding respectively to ξ = ξ L and ξ = ξ h . The two metrics are therefore related as Comparing this with (3.37), we arrive at the following mapping of the parameters of the dual frame Lifshitz background corresponding to a given hvLf background: (3.39) In practice we are interested mostly in the case ξ h = 0, so that the hvLf metric solves the equations of motion in the Einstein frame. This relation between Lifshitz and hvLf solutions can be utilized in order to transform such hvLf backgrounds into Lifshitz backgrounds in a different Weyl frame. This is exactly analogous to the way Dp branes with p ≤ 4, were studied in [18] by going to a Weyl frame where the geometry is asymptotically locally AdS. The method we develop in the following in order to systematically construct the holographic dictionary is directly applicable to Lif backgrounds in any Weyl frame and to hvLf backgrounds in the Einstein frame with θ < 0. This restriction for hvLf in the Einstein frame is related to the fact in the coordinate system (3.27) the UV is located at r = 0 for θ > 0. However, for any θ, we can work in the dual frame where the hvLf backgrounds become Lifshitz. We will therefore work entirely in the dual frame from now on and consider Lifshitz asymptotics only. In this way we are able to develop the holographic dictionary for both Lif and hvLf with any θ simultaneously. This is the reason for allowing for a non-zero Weyl parameter ξ throughout our analysis. It is useful to keep in mind that the combination of parameters dµξ from now on can be understood as dµξ = −θ, (3.40) where θ is the hyperscaling violating parameter of the corresponding hvLf background in the Einstein frame.

Lif boundary conditions as a second class constraint
From the solutions (3.18) follows that Lifshitz boundary conditions amount to the asymptotic relationṡ Inserting these asymptotic expressions in the flow equations (3.12) and (3.11), one finds that the resulting set of linear PDEs for U (X, Y ) admit an asymptotic solution for the superpotential U (X, Y ) provided asymptotically The corresponding superpotential U (X, Y ) takes the form It is important to pause for a moment and clarify the significance of these asymptotic conditions since they play a key role in the construction of the holographic dictionary for anisotropic backgrounds and throughout the subsequent analysis. Using the definition of the variable Y we can express the time component of the vector field as This expression can be seen as a change of variables in configuration space (a special canonical transformation), trading the variable a in favor of Y − Y o , without any physical significance. The non-trivial condition, however, comes from demanding Lif asymptotics, i.e. that asymptotically Y − Y o → 0. The reason why this is particularly significant is that setting is not compatible with any integral of motion of the equations (3.2) and so amounts to a second class constraint. Another way this constraint can be deduced is the fact there is no superpotential U (X)crucially without any dependence on Y -that leads to the asymptotics (3.41) via the flow equations (3.11) and (3.12). In Appendix A we show how such a constrained system can be described in a Hamiltonian language, either by solving explicitly the constraint at the start, or by using Dirac's algorithm for constrained systems. As long as we keep at least the linear term in Y − Y o in (3.43), which corresponds to a deviation from the constraint surface (3.45), the standard Hamiltonian analysis applies, however. Demanding that a Taylor expansion in Y − Y o be compatible with the dynamics is equivalent to requiring that (3.45) be a consistent truncation of the theory. In other words, we are asking that the effective potential 2 for the fluctuation Y − Y o has no linear term and that the quadratic term (mass) is such that Y − Y o sources a relevant operator. As we will see shortly, this leads to further conditions for the potentials parameterizing the Lagrangian (3.3), besides the leading asymptotic form (3.17).

Fefferman-Graham expansions and anisotropic RG flows from a superpotential
In the previous subsection we determined that imposing Lifshitz asymptotics requires the superpotential U (X, Y ) to have the asymptotic form (3.43). In order to obtain asymptotically Lif backgrounds that correspond to deformations of the 'ground states' (3.18), such as anisotropic renormalization group (RG) flows, we need an exact solution of the superpotential equation (3.10) that satisfies the asymptotic condition (3.43). In this subsection we make use of various ansätze to simplify the superpotential equation and we present a class of exact solutions corresponding to a certain marginal deformation of the backgrounds (3.18).
We also obtain the general solution to the superpotential equation (3.10) with the asymptotic condition (3.43) in the form of a Taylor expansion in Y − Y o , which can be used to determine the general asymptotic complete integral and the Fefferman-Graham expansions. It is worth pointing out that a solution U (X, Y ) of (3.10) cannot be polynomial in Y for the physical range of the various parameters. Combined with the asymptotic condition (3.43), this implies that any superpotential can be expressed as a non-truncating Taylor series in Y − Y o , although there can be non-analytic terms starting at the normalizable order.

Superpotential I:
An important special case of the Lagrangian (3.3) occurs when the potentials are exactly -not merely asymptotically -exponentials as in (3.17), i.e.
with the various parameters satisfying the relations (3.19). Since this holds asymptotically anyway, this example captures the essential physics for general asymptotically Lif and hvLf backgrounds. The superpotential equation (3.10) in this case can be reduced to an ordinary differential equation for some function w(y) of y ≡ Y Z ξ (X). Inserting the ansatz (3.47) into (3.10) we get a first order ODE for w(y): where The asymptotic condition (3.43) determines that w(y) must satisfy Equation (3.48) can be transformed into an Abel equation of the first kind [54], which is in general non-integrable. For special ranges of the parameters it admits analytic solutions of the form w = √ a + by, which are special cases of the more general class of solutions derived from superpotential II below. For generic values of the parameters, however, we can obtain the solution to (3.48) subject to the initial conditions (3.50) in the form of a Taylor expansion around y o , including potential non-analytic terms at normalizable order. In particular, for generic values of the parameters the solution of (3.48) subject to the initial conditions (3.50) takes the form 3 where w 0 and w 1 are determined by (3.50), w 0 is an integration constant, and the scaling dimension When d + z + dµξ − ∆ + = 0 we have instead y − y o ∼ r −1 and so y − y o is the source of a marginally relevant operator in this case. As we shall see, the value of ∆ + is related to w 2 , which is determined by the quadratic equation The two roots of this equation are For ξ = µ = 0 these roots reduce to The coefficients w n with 2 < n < n c : can be obtained recursively from the linear equations These are all the terms that are needed to determine the asymptotic solutions of the fields via the flow equations, since the terms w n with n > n c , as well as the terms w n with n ≥ 1, are subleading relative to the normalizable modes. When ∆ + = d + z + dµξ, however, the mode y − y o ∼ 1/r goes to zero only logarithmically and n c → ∞, which means that all terms in the solution (3.51) must be kept in this case to obtain the correct asymptotic solution of the HJ equation. This is reminiscent of what happens in the case of Improved Holographic QCD [37,57,58] and it is important in order to correctly renormalize the often studied Einstein-Proca theory for d = z = 2 when the marginally relevant deformation y − y o is turned on as in e.g. [59]. All terms must also be determined in order to obtain an exact background solution. Backgrounds with w n = 0 can be obtained through the recursion relations (3.58) applied to any n > 2.
These results are in agreement with those of [1,27,31], which were obtained through an analysis of the linearized fluctuations of the equations of motion. Here we have derived these in a simpler way using only the superpotential equation (3.10). There is no need for studying linearized perturbations of the equations of motion (except for computing 2-point functions, of course), or indeed using the second order equations, since the full asymptotic expansions can be obtained from an asymptotic complete integral of the HJ equation.
Inserting the solution (3.51) in the flow equations (3.11) and (3.12) leads to the first order equationṡ Note that for ξ = µ = 0 we have = z and so These first order equations can be integrated to obtain the full set of asymptotic expansions, including the normalizable and non-normalizable modes. In particular, the non-normalizable modes appear as integration constants of these first order equations. Namely, the leading asymptotic form of the fields takes the form where φ o , c 3 , c 4 , c 5 and a o are integration constants, and we have kept the notation of [27] to facilitate the comparison of the modes. However, a o is fixed by the boundary condition (3.45) in terms of the other parameters as It corresponds to a source of a marginal operator with respect to Lifshitz boundary conditions, which do not want to turn on. Moreover, if ∆ + ≥ d + z + dµξ, then the mode c 3 must also be set to zero since otherwise is not vanish asymptotically and the Taylor expansion in Y − Y o breaks down. In terms of the dual theory, in that case c 3 sources a marginal or irrelevant operator relative to the Lifshitz theory. Finally, this asymptotic form of the scalar is valid assuming µ = 0. If µ = 0 then one has to look at subleading terms of the potential, and in particular at the mass term, to determine the asymptotic form of the scalar.
To determine the normalizable modes we need to consider the most general deformations of the solution (3.51) of the HJ equation, as was discussed in Section 3.1. We showed that this can be done by finding the characteristics of the linear PDE defined by the dilatation operator. The dilatation operator itself is obtained from the asymptotic form of the non-normalizable modes through Note that the Lifshitz boundary condition has changed the form of the dilatation operator, replacing the derivative with respect to a with a derivative with respect to Y − Y o . This reflects the fact that Lifshitz boundary conditions fix the mode a o and so we cannot consider variations with respect to a o without changing the variational problem. To determine the normalizable modes, therefore, we need to find the characteristics of the linear PDE Assuming µ = 0, a convenient basis for the three independent characteristics is and so the most general 5 deformation of the solution (3.51) of the HJ equation can be written in the form where q i are the normalizable modes. 6 Note that the parameter w 0 in (3.51) can be expressed in terms of these deformation parameters. The fact that there are only three independent normalizable modes, while there are apparently four sources is due to the fact that we consider homogeneous solutions. A fourth deformation of the HJ solution is the energy, but such a deformation is not allowed in a model that comes from gravity since the Hamiltonian vanishes due to diffeomorphism invariance. The source conjugate to the energy is the radial cutoff r o , which can be used to eliminate one of the sources for homogeneous solutions. We choose to eliminate φ o . From (3.67) we find that the symplectic form on the space of asymptotic solutions [45] takes the form As we shall see in Section 5, the modes q 1 , q 2 and q 3 are related respectively to the energy density, spatial stress tensor and scalar operator dual to Y − Y o [27]. Finally, from the momenta (3.5) we see that the deformations (3.67) will modify the flow equations (3.59) according to The most general deformation, of course, corresponds to adding arbitrary functions of these characteristics. However, we are only interested in a complete integral and for this it suffices to consider constant coefficients multiplying a given function of the characteristics. 6 This is a special case of Sreg in (5.20) in Section 5 for general asymptotically hvLf backgrounds.
where we have used the fact that the sources a o and φ o have been fixed. Since these terms correspond to the normalizable modes in the asymptotic expansions, the latter are only needed up to this order.

Superpotential II:
We now consider an ansatz that allows us to separate variables in the superpotential equation (3.10), and as a result, to obtain exact hvLf solutions that correspond to marginal deformations of the backgrounds (3.18). Inserting the ansatz where ε 0,1,2 = ±1 are independent signs, in the superpotential equation (3.10) leads to the three equations (3.71) The first and second equations can be integrated directly to obtain v = ± √ αˆX dX ε 2 W (X ), However, u must also satisfy the last equation in (3.71), which leads to a constraint relating V (X), Z(X) and W (X). Any solution of these equations is a solution to the original superpotential equation (3.10), but in order for this superpotential to correspond to Lif or hvLf solutions the asymptotic conditions (3.43) must also be satisfied. Expanding the ansatz (3.70) around the asymptotic curve (X, Y o (X)) we obtain Comparing these with the asymptotic conditions following from (3.43) determines Inserting the asymptotic condition for v 2 in the first equation in (3.71) leads to a constraint on the parameters of the solutions, namely Before determining the possible solutions of this constraint, it is instructive to derive it in an alternative way. Inserting the ansatz (3.70) in the flow equations (3.12) (3.11) we obtain where Y ≡ e −2ξX Y and the radial coordinate ρ is defined by Combining the first two flow equations leads to a first order equation for Y as a function of X: This is an Abel equation of the second kind [54], which is in general non-integrable but there are known integrable classes. In particular, this equation can be solved for the u and v in (3.75). The solution is where c is an integration constant. Since d − θ = −(z − 1) (otherwise u and v vanish identically), the only way this solution can be compatible with the asymptotic condition (3.42) is that the parameters of the solution satisfy α + ν(ν + ξ) = 0 and the integration constant is chosen appropriately so that Y = Y o (X) identically. It can be checked that this condition on the parameters is precisely the constraint (3.76). It is also the condition for the dimension ∆ + in (3.60) to be equal to d + z − θ and therefore, the operator dual to the deformation Y − Y o is a marginal operator. Indeed, (3.80) can be written as for an arbitrary constant c and so The boundary condition (3.42), however requires that we turn off the source for this operator and so we must set c = 1. With the source for Y − Y o set to zero the corresponding background solutions are identical to the backgrounds (3.18), but for the specific set of parameters that satisfy (3.76). However, turning on a source for Y − Y o in this case leads to a marginal deformation of the dual theory, which can be seen as a shift in the exponent . The parameter space allowed by the marginality condition (3.76) turns out to be rather restricted, but non-empty. One can show that there is no solution with µ = 0 and finite ξ, or with ξ = 0. Solving the constraint for αµ 2 in terms of µξ, d and z we get Recall that −dµξ is the hyperscaling violating exponent θ in the Einstein frame, while αµ 2 ≥ 0 is related to the independent vector hyperscaling violating parameter discussed in [1].
i) z > 1: For z > 1 we must choose the plus sign in (3.82). The quantity inside the bracket then is positive provided either of the following two conditions holds: The first condition requires which is compatible with the NEC provided In terms of θ these solutions can be summarized as follows: ii) z < 1: For z < 1 the minus sign in (3.82) must be chosen. The RHS of (3.82) is then positive provided which violates the NEC except for the limiting case θ = d as above, but now with z ≤ 0.
Superpotential III: As a final example, we consider the Taylor expansion of the general superpotential U (X, Y ), without any simplifying assumptions for the potentials of the Lagrangian except for the asymptotic conditions (3.17). However, as we already anticipated, additional consistency conditions will arise by requiring that a Taylor expansion in Y − Y o be consistent with the asymptotic expansion, as required by the Lifshitz boundary conditions. The analysis here is a straightforward generalization of the analysis for superpotential I above.
We start by expanding the superpotential In order to simplify the subsequent formulas we reparameterize the coefficients U m (X) as In fact, there are three distinct requirements this superpotential must fulfill in general: The asymptotic form (3.43) of the superpotential determines the asymptotic behavior of the coefficients u 0 (φ) and u 1 (φ) to be More generally, where w n are the coefficients of the Taylor expansion (3.51).
ii) Hamilton-Jacobi equation Inserting the formal Taylor expansion in the superpotential equation (3.10) leads to a set of equations for the coefficients u m (φ). The first three orders in Y − Y o give respectively Note that these equations alone do not completely determine the functions u n (φ) in the Taylor expansion of the superpotential.

iii) Consistency of the Taylor expansion
A final condition on the functions u n (φ) is imposed by requiring that the Taylor expansion is consistent with the asymptotic expansion. To derive this consistency condition we need to write the flow equations (3.11) and (3.12) in terms of the functions u n (φ), namelẏ The consistency condition comes from the inhomogeneous term in the flow equation for Y − Y o , which must vanish identically in order for the Taylor expansion to be well defined. Note that if the inhomogeneous term is not zero then does not vanish asymptotically. This condition holds automatically for the asymptotic form (3.91) of u 0 and u 1 and the leading form of Z in (3.17), but it imposes a non-trivial condition on the subleading terms of u 0 and u 1 (or of Z if one views this as an equation for Z.) These three conditions on the superpotential completely determine the coefficients u n (φ) in the Taylor expansion. Notice that the inhomogeneous term in the Y −Y o flow equation is identical to the coefficient of u 2 and u 3 in (3.94) and (3.95) respectively. Since this term must vanish, u 2 is eliminated from (3.94) and u 3 from (3.95). Equations (3.93) and (3.94) then become two equations for u 0 (φ) and u 1 (φ), while (3.95) becomes a Riccati equation for u 2 (φ). Higher order terms are determined by first order linear equations that are derived from higher orders in Y − Y o of the HJ equation. Since u 0 (φ) and u 1 (φ) must also satisfy the constraint coming from the consistency of the Taylor expansion, there are three equations for these two functions, and hence there is an implicit constraint on the three potentials V , W and Z. The three equations are 7 However, in a bottom up approach the potentials V , W and Z are a priori unspecified and so we can in fact define the potentials in terms of the two functions u 0 (φ) and u 1 (φ) of the superpotential, which are only subject to the asymptotic conditions (3.91). Given these functions, the Riccati equation (3.95) can be solved for u 2 and the higher order coefficients u n are determined by solving the linear equations coming from the higher order terms in the Taylor expansion of the HJ equation. The leading asymptotic form of these will be identical to the one obtained from the superpotential I above, but they can potentially differ at subleading orders due to the choice of subleading terms in u 0 (φ) and u 1 (φ). Finally, the Fefferman-Graham asymptotic expansions are obtained by integrating the flow equations (3.96) and (3.97). Note that since the leading asymptotic form of these expansions is the same as for the superpotential I above, the non-normalizable modes remain the same as in that case. Moreover, since the form of dilatation operator is determined by the non-normalizable modes, it follows that the analysis of the finite part of the asymptotic complete integral, and hence the normalizable modes, are again the same as in the superpotential I case.
The only exception occurs in the case µ = 0, where the subleading terms in u 0 (φ) and u 1 (φ) determine the asymptotic form of the scalar. But the corresponding normalizable and non-normalizable modes can be determined by the same procedure in that case too.

Recursive solution of the HJ equation for asymptotically locally Lif backgrounds
In the previous section we considered exclusively homogeneous backgrounds, for which we obtained the general asymptotic solution of the Hamilton-Jacobi equation, the Fefferman-Graham expansions, as well as the non-normalizable and normalizable modes corresponding respectively to the sources and 1-point functions of the dual operators. We now extend this analysis to incorporate sources with arbitrary spatial and time dependence. Note that the solution of the HJ equation we obtained in Section 3 is still relevant in the presence of arbitrary spacetime-dependent sources, since it appears as the leading zero derivative solution of the HJ equation. What we will be mainly concerned with in this section, therefore, is the systematic construction of the subleading terms in the HJ solution that contain transverse derivatives.

Locally Lif boundary conditions
Before we address the derivative terms in the solution of the HJ equation, however, we need to identify the most general spacetime-dependent sources allowed by Lifshitz boundary conditions. To this end we consider again the most general diffeomorphism and gauge invariant solution of the general HJ equation (2.17), containing no transverse derivatives. As we have argued in the previous section this takes the form should appear in the superpotential and not A i A i , and so S (0) in fact contains transverse derivatives, but in a rather trivial way. The relation between the superpotential U (X, Y ) and the asymptotic form of the fields is provided by the flow equations (2.20), which now becomė In order to accommodate anisotropic solutions we parameterize the induced fields on the radial slice Σ r in terms of fields compatible with the anisotropy. In particular, we decompose the induced metric γ ij and vector field A i as 8 where the indices a, b run from 1 to d and σ ab (r, t, x), n a (r, t, x), n(r, t, x), a(r, t, x) and A a (r, t, x) are the fields in terms of which we will parameterize the dynamics. In terms of the anisotropic fields the flow equations (4.2) take the form where we have used the leading asymptotic form of the flow equation for the Stückelberg field. The Lifshitz metric (3.18) implies that the most general asymptotic form of the fields n and n a compatible with locally Lif asymptotics is where n (0) (t, x), n (0)a (t, x), and g (0)ab (t, x) are arbitrary functions of the transverse coordinates and the constant β is to be determined. Since γ tt = −n 2 + n a n a , requiring that n a n a is at most divergent as n 2 imposes the restriction Inserting the asymptotic behaviors (4.5) in the flow equations (4.4) leads to a set of asymptotic conditions on the superpotential, namely (4.10) Using the inverse metric n a n 2 n a n 2 σ ab − n a n b n 2 , (4.11) where we have used (4.8) in the last step. Inserting this in (4.7) gives Moreover, using the leading form of the flow equation for ω to replaceȧ andȦ a withḃ andḂ a respectively in the vector flow equations, we see that the latter require that the time component, b, and the spatial component, B a , behave in the same way asymptotically, which we parameterize as where b (0) (t, x) and B (0)a (t, x) are arbitrary functions of the transverse coordinates and the exponent is as yet unspecified. Using this asymptotic form of B a in the vector flow equation together with (4.12) we find which is the asymptotic constraint (3.42) we found for the homogeneous solutions. Moreover, n a B a ∼ n (0) a B (0)a e (z−1−β+ )r , (4. 16) and so where, assuming B (0)a = 0, (4.18) However, (4.12) implies that, if B (0)a = 0, in order to satisfy (4.8) we must demand that which requires that either z < 1 or β < 0. The latter contradicts the above asymptotic conditions and so it is not an acceptable solution. Moreover, we have argued that z < 1 corresponds to the solutions I and II of the NEC in (1.8) and since θ ≥ d + z in those cases, there are no well-defined asymptotic expansions. A possible exception is the marginal case θ = d + z with 0 ≤ z < 1, but we will not consider this here. The only alternative, therefore, is to require Note that the inequality (4.19) need not hold in this case since (4.8) is automatically satisfied. Moreover, (4.10) determines (z − 1 − β)n (0)a = 0, (4.22) in this case, which can be solved by either setting β = z − 1 and leaving n (0)a (t, x) arbitrary, or by setting n (0)a (t, x) = 0 in which case β does not arise at all. Since we want to keep all possible sources compatible with Lif asymptotics, we set and keep n (0)a (t, x) unconstrained. To summarize, from this asymptotic analysis we have determined that locally Lifshitz boundary conditions amount to the gauge-invariant asymptotic constraint where n i = (n, 0) is the unit normal to the constant time surfaces and Y o (X) is defined in (4.15). This is a covariant way of writing the scalar constraint (4.15) and the spatial vector constraint (4.20). This covariant form of the asymptotic constraint allows us to obtain the corresponding asymptotic form of the covariant momenta which can be integrated to obtain the leading asymptotic from of the zero order solution of the Hamilton-Jacobi equation: The asymptotic form of the momentum conjugate to the Stückelberg field ω following from this HJ solution is which as we shall see shortly is subleading relative to the rest of the momenta in a precise sense that we will specify. In terms of the superpotential, the asymptotic conditions (4.25) imply the following conditions on the superpotential U (X, Y ) and its first derivatives: (4.28) Inserting these in the superpotential equation (3.10) one recovers the relations (3.19) between the various parameters. As we have seen from the homogeneous solutions in Section 3, there are additional constraints on the superpotential at subleading orders, coming from the consistency of the Taylor expansion in B i −B oi . Moreover, there are more sources appearing at subleading order due to the constraint (4.24). We will revisit these points later on, when we develop the recursive algorithm for determining the subleading terms of the HJ solution and when discussing the general Fefferman-Graham expansions.

Graded expansion in eigenfunctions of the derivative and gradation operators
A solution of the HJ equation of the form (4.1) captures all zero derivative terms. However, the general asymptotic solution of the HJ equation with spacetime-dependent sources contains asymptotically subleading terms with transverse derivatives acting on the induced fields. In order to account for these terms in a systematic way, and to consistently impose Lif boundary conditions, we are going to seek a solution in the form of a covariant expansion in eigenfunctions of a suitable functional operator. This is analogous to the expansion in the dilatation operator for asymptotically locally AdS spaces introduced in [44] or its generalization to asymptotically non AdS -but relativistic -backgrounds in [37]. The anisotropy introduced by the Lif boundary conditions, however, necessitates some generalization of the formalism. The dilatation operator method has been extended to Lifshitz backgrounds without a linear dilaton in the vielbein formalism [29] and in Lifshitz gravity [33]. However, the expansion we develop is both fully covariant and applicable in the presence of a linear dilaton, which is necessary in order to accommodate hvLf backgrounds.
The leading order solution of the Hamilton-Jacobi equation in this covariant expansion is of the form (4.1). Since the superpotential U (φ, B 2 ) depends on the choice of the potentials V (φ), Z(φ) and W (φ) in the Lagrangian, which we want to keep as general as possible at this stage, we demand that (4.1) be an eigenfunction of the functional operator we expand in for any choice of U (φ, B 2 ). There are two operators that satisfy this criterion, namely for which it is easy to check that and so S (0) is an eigenfunction of both δ and δ B , with respective eigenvalues d + 1 and 1, for any U (φ, B 2 ). Crucially, these operators commute which means that if S (2k) is an eigenfunction of δ, then so is δ B S (2k) with the same eigenvalue. This allows us to expand S covariantly in a double expansion. In order to construct the covariant expansion, we need to understand the structure of the eigenfunctions of δ and δ B . As we have argued, any function of B 2 (and trivially of φ) is automatically an eigenfunction of both operators. It therefore remains to understand how these operators act on terms with transverse derivatives, ∂ i . From the structure of the Hamiltonian constraint follows that any derivative expansion of the Hamilton-Jacobi functional will contain only even number of derivatives. Covariance then requires that for every pair of derivatives there is either an inverse metric, γ ij , or a factor of B i B j with which the two derivatives are contracted. A simple counting exercise then shows that δ counts the number of derivatives. Namely, any functional S (2k) containing 2k derivatives is an eigenfunction of δ with eigenvalue d + 1 − 2k, where d + 1 is the contribution of the volume element.
The eigenvalues of the operator δ B follow from the observation that it satisfies where is a projection operator: This implies that an eigenfunction S (2k) of δ with 2k derivatives can be split in a sum of up to k + 1 terms containing 0, 1, . . . , k powers of σ ij . This can be achieved systematically as follows. Terms in which all 2k derivatives are contracted with B i are eigenfunctions of δ B with eigenvalue 1 − 2k, since every factor of B i contributes −1 to the eigenvalue and the 1 comes from the volume element. Next, we consider terms where 2k − 2 derivatives are contracted with B i and 2 derivatives are contracted with γ ij . Such terms are not eigenfunctions of δ B but they can be written as a sum of two eigenfunctions of δ B with eigenvalues 1 − 2(k − 1) and 1 − 2k by writing where and S (0,0) = S (0) is given by (4.1). We will refer to the operator δ as the 'derivative operator' since it counts transverse derivatives, while δ B we will call the 'gradation operator'. It should be stressed, however, that there is an inherent assumption of locality for these expansions in local eigenfunctions of the operators δ and δ B to be meaningful. This assumption is of course not valid for the finite part of the solution of the HJ equation, i.e. the renormalized on-shell action. However, this is of no concern right now. Our strategy is to develop a recursive algorithm that determines iteratively increasingly asymptotically subleading terms in the solution of the HJ equation assuming locality. This recursive procedure breaks down exactly at the order where the finite contribution to the solution occurs. This finite part is required in order for the asymptotic solution of the HJ equation to qualify as a complete integral, and it is necessary for the derivation of the Fefferman-Graham expansions and the identification of the normalizable modes. As in the case of homogeneous solutions in Section 3, the finite non-local part must be addressed separately, and it will be the main subject of Section 5. Table 1. Action of the operators δ and δB on the canonical momenta.

Expansion of the canonical momenta
Since the canonical momenta are related to the solution of the Hamilton-Jacobi equation via (2.18), one might expect that the momenta defined via are also eigenfunctions of δ and δ B . This is in fact not true, and it should be emphasized that the subscripts in the momenta do not denote their eigenvalues under δ and δ B , since they are not eigenfunctions. The subscripts on the momenta instead indicate that they are gradients of the corresponding eigenfunctions S (2k,2 ) . The action of δ and δ B on these momenta can be obtained using the commutation relations (4.39) The results are summarized in Table 1. From the expressions in Table 1 the complete set of linearly independent eigenfunctions of both δ and δ B that are linear in the canonical momenta can be constructed. These eigenfunctions are listed in Table 2, along with their eigenvalues under δ and δ B . The eigenfunctions in Table 2 in turn allow us to decompose any quantity that involves the canonical momenta in terms of these eigenfunctions. For example, the metric and vector momenta can be decomposed in terms of eigenfunctions of δ and δ B as follows: where the quantity P (2k,2 ) i is defined in Table 2. For future reference we decompose all scalar quantities that are quadratic in the canonical momenta in terms of the eigenfunctions of these operators in Table 3. We will need these eigenfunctions in the next subsection in order to analyze the Hamiltonian constraint and to develop the recursion algorithm.

Expansion of the first class constraints
In order to develop a recursive algorithm for solving the Hamilton-Jacobi equations in terms of eigenfunctions of the derivative and gradation operators we must expand the first class constraints (2.16) in eigenfunctions of these operators. The momentum and U (1) gauge constraints are linear in the momenta and so they can be decomposed in eigenfunctions of δ and δ B using the eigenfunctions in Table 2. The Hamiltonian constraint, however, is quadratic in the momenta and the eigenfunctions in Table 3 are required instead. Let us consider each constraint in turn.
U (1) constraint: can be immediately decomposed in eigenfunctions of δ and δ B using the last eigenfunction in Table 2. Namely, k, and hence π ω (2k,2 ) = D i π (2k,2 ) i , ∀k, . Momentum constraint: Using the U (1) constraint we can write the momentum constraint in the form which can be expanded in eigenfunctions of δ so that for all k Using the decomposition of the momenta in eigenfunctions of both the derivative and gradation operators in (4.40), this can be in turn written as Matching terms of equal eigenvalues under δ B we obtain the two conditions for all 0 ≤ ≤ k. In particular, we note the special cases (4.48) Hamiltonian constraint: The Hamiltonian constraint in (2.16) is quadratic in the canonical momenta and it is the dynamical equation that determines the Hamilton-Jacobi function S. In particular, using the decomposition of the momenta in terms of the eigenfunctions of δ and δ B , we will turn the Hamiltonian constraint into a tower of linear equations for S (2k,2 ) , which can be solved iteratively.
Expanding the Hamiltonian constraint in eigenfunctions of δ and isolating terms with the same eigenvalue we obtain for k > 0 where We have written these constraints in the form of inhomogeneous linear equations for S (2k) by collecting all momenta coming from S (2k) on the LHS and grouping terms that originate in S (2k ) with k < k in the inhomogeneous term R (2k) . There is an exception to this, however, because as we have seen above the δ eigenvalue of 1 , and therefore, this term must be included in the source R (2k+2) . Inserting the the zero order momenta in these recursion relations we obtain (4.52) Finally, using Tables 2 and 3 these recursion relations can be expanded in eigenfunctions of δ B as for all k > 0 and 0 ≤ ≤ k. These recursion relations are the basis of our algorithm for systematically solving the Hamilton-Jacobi equation. We now explain how this can be achieved.

Recursion relations
We now turn to the question of how the recursion relations (4.53) can be utilized in order to determine the terms S (2k,2 ) of the Hamilton-Jacobi functional. A number of useful results that we will need in this section is presented in Appendix B. In particular, in the appendix we define the unintegrated versions of the functional operators δ and δ B , namely, Using these unintegrated operators we can rewrite (4.53) in the form This form of the recursion relations allows us to utilize the fact that S (2k,2 ) is a simultaneous eigenfunction of both δ and δ B . Some attention is required, however, in understanding the structure of various total derivative terms. Writing S (2k,2 ) =ˆd d+1 xL (2k,2 ) , (4.56) and using the results of Appendix B, we have as well as where we have invoked Lemma B.1 to deduce that L (2k,2 ) is an eigenfunction of δ, without any total derivative term. Combining these relations one can show that the operators δ and δ B act on the total derivative terms as follows: However, L (2k,2 ) is only defined up to a total derivative and so we are free to define Using the action of δ and δ B on the total derivative terms we now find where and it satisfies δ u (2k,2 ) More generally we define where λ is an arbitrary parameter, so that Inserting these expression in the recursion relation (4.55) we obtain where we have dropped the superscript λ in L λ (2k,2 ) . Provided the ratio of the functions Y U Y and α ξ U − 2(α ξ + d 2 ξ 2 )Y U Y + dξU X is constant, a suitable choice of the parameter λ eliminates the total derivative term. However, we will keep the total derivative term for the time being and proceed with solving these recursive equations. On the way we will determine the minimal condition the superpotential U (X, Y ) must satisfy so that this total derivative term can be eliminated.

Taylor expansion in the Lifshitz constraint
The expansion of the HJ functional in eigenfunctions of the commuting operators δ and δ B and the corresponding recursion relations (4.66) are not specific to Lif boundary conditions. In order to incorporate these we must impose the asymptotic constraint (4.24). This means that, in addition to the expansion in eigenfunctions of δ and δ B , the solution of the HJ equation must take the form of a Taylor expansion in B i − B oi . In particular, these two expansions must be consistent with each other, and so each term S (2k,2 ) in the graded covariant expansion must admit a Taylor expansion in B i − B oi . This Taylor expansion, except from imposing Lif boundary conditions, will allows us to eliminate the functional derivative with respect to B i in the recursion relations (4.66), leading to tractable linear functional differential equations in one variable.
The Taylor expansion in B i − B oi for the zero order solution S (0) can be immediately obtained from the Taylor expansion of the superpotential U (X, Y ) in Y − Y o in Section 3, using the identity However, since the operators δ and δ B depend on B i as well, they must also be Taylor expanded. Considering δ first, we evaluate and we have made use of the identity (C.8) in the third line. An analogous result holds for d B . This leads to the following identities where the operators are respectively the pullbacks of the operators d and d B on the constrained submanifold B i = B oi . Note that since B oi ∝ n i , the unit normal to the constant time slices, it follows that the pullback of the gradation operator, δ 0 B , counts time derivatives. Moreover, the pullback of the projection operator (4.33) becomes the spatial metric (see Table 12) The covariant expansion in simultaneous eigenfunctions of δ 0 and δ 0 B , therefore, is a derivative expansion with the number of derivatives given by the eigenvalue of δ 0 and graded according to the number of time derivatives, counted by the eigenvalue of δ 0 B .

Taylor expansion of the HJ equation
The HJ equation for the zero order solution L (0) is the superpotential equation (3.10). Since L (0) depends on B i only though Y = B i B i the Taylor expansion of the superpotential equation in B i − B oi is equivalent to the Taylor expansion in Y − Y o we discussed in the superpotential III part of Section 3. All the results there carry over, except that the flow equations must be generalized to account for components that were identically zero for homogeneous backgrounds. For now, we only need equations (3.93), (3.94) and (3.95), which follow from the Taylor expansion of the superpotential.
The HJ equations for L (2k,2 ) with k > 0 are the recursion relations (4.66). Inserting the expansion (4.68) and using the identity (C.7) the first two orders in B i − B oi give the following two equations: It must be stressed that with this definition of π 0ij (2k,2 ) and π 0 φ(2k,2 ) these quantities are not the O(B −B o ) 0 terms in the Taylor expansion of the corresponding momenta. In fact, using (C.8) and (C.7) we find (4.78) We will not present the equations for O(B − B o ) 2 and higher here, but note that provided B i − B oi sources a relevant operator, there is always some order at which the Taylor expansion can be truncated since higher order terms are subleading relative to the normalizable modes. At which order the Taylor expansion can be truncated depends on the leading asymptotic behavior of B − B o , which was discussed in Section 3. Moreover, we can identify some generic features that apply to the higher order equations as well. Firstly, recall that the Taylor  . Another generic feature of these equations is the structure of the total derivative terms. In particular, the relative coefficient of the two total derivative terms remains the same for any order. It follows that imposing a single condition on the functions u 0 (φ) and u 1 (φ), in addition to the three equations (3.98), ensures that the total derivative terms can be eliminated from all equations at any order. Namely, if holds for some constant c, then the total derivative terms can be eliminated by setting . (4.80) The constant c cannot take any value, however, since the asymptotic conditions (3.91) require that We will therefore restrict our attention to theories that satisfy in addition to (3.98). Using the third equation in (3.98), this condition (4.82) can alternatively be written Imposing this relation between u 1 and u 0 implies that the functions V (φ), W (φ) and Z(φ) are all parameterized in terms of one arbitrary function through (3.98). Note however, that (4.82) is automatically satisfied by the asymptotic form (3.91) of the functions u 1 and u 0 and so it imposes no additional constraint on the parameters of generic Lif solutions. It only constrains the structure of the subleading terms in u 1 and u 0 and in this sense it is a mild restriction. However, we believe that imposing this restriction is not essential in order to solve the equations (4.74) and (4.75), but we have found no alternative way to solve them in the generic case. Of course, in special cases one can use an ansatz to solve these equations, but besides being very inefficient, this approach cannot be applied to the general case. Incorporating the conditions (4.82) and ( which is defined provided u 0 + Z Z u 1 = 0. Moreover, the source B oj R 1j (2k,2 ) in the last two equations is given by  O(B − B o ). These equations provide a recursive algorithm that allows us to obtain the solution of the HJ equation at order k + 1 from the solution at order k. Namely, given the solution of the HJ at order k, the corresponding canonical momenta determine the inhomogeneous term in the linear equations for the order k + 1 solution. The main technical challenge in this algorithm is solving these recursion relations. Obtaining the canonical momenta from a given solution and constructing the inhomogeneous term for the next order can also be tedious, but it's straightforward. As we will see momentarily, the solution of the recursion relations can be streamlined using the integration technique developed in [37]. Solving the HJ equation then becomes entirely algorithmic and it is ideally suited for implementation in a symbolic computation package such as xAct [60]. The recursion relations (4.84), (4.85) and (4.86) are identical in form to the equations appearing in the recursive solution of the HJ equation for relativistic backgrounds [37] and exactly the same techniques can be applied here. Indeed, many of the results in [37] are directly relevant. Firstly, note that the solutions of (4.84), (4.85) and (4.86) are qualitatively different depending on whether u 0 + Z Z u 1 is zero or not. Using (3.91) we see that this quantity asymptotes to the constant parameter µ and so there are three cases to examine: i) µ = 0, ii) µ = 0 but u 0 + Z Z u 1 not identically zero, and iii) u 0 + Z Z u 1 = 0, at least up to normalizable modes. We will consider two examples of case iii) in Section 6. We will not discuss case ii) further here because it requires a specification the subleading terms in u 1 and u 0 that determine the asymptotic form of the scalar in this case. This can be easily done but would take us away from the generic case. In this section we will instead focus on case i), which is the generic situation.
Provided the parameter µ is not zero, all recursion relations (4.84), (4.85) and (4.86) admit a the homogeneous solution of the form which implies that the homogeneous solution is finite and so it corresponds to the usual renormalization scheme dependence. 9 The inhomogeneous solutions of (4.84), (4.85) and (4.86) can be written formally in the form As in Eq. (2.36)-(2.37) of [37], the expressions (4.94) for the inhomogeneous solutions are formal since the source terms, such as R 0 (2k,2 ) [γ,φ], generically contain derivatives of the scalar φ. In [37] these formal integrals were defined by systematically tabulating all possible derivative structures involving the scalar, Table 4. General integration identities for integrands that contain up to four derivatives on the scalars that were derived in [37]. The shorthand notation ffl φ k, ,m is defined in (4.96). R (2k,2 ) stands for any of the source terms on the RHS of (4.94), while L (2k,2 ) stands for any of the quantities on the LHS. The tensors t i 1 i 2 ...im and t ij are arbitrary totally symmetric tensors independent of φ, while t ijkl These formulas suffice for all terms appearing in R 0 (2,0) and R 0 (2,2) , but only for terms in R 0 (4,0) , R 0 (4,2) and R 0 (4,4) that are contracted with the particular tensors t ijkl 1 and t ijkl

2
. Although these tensors cover the most general 4-derivative terms in the relativistic case [37], this is not in general the case for the non-relativistic boundary conditions we impose here. However, the relevant integration formulas that generalize this table can be derived as in [37]. Moreover, as we will see in Section 6, these formulas are not required in the case of exponential potentials, since the integrals over the scalar can be evaluated in general independently of the tensor structure in that case.
up to four derivatives, and the corresponding integrals were evaluated generically. The results, adapted to the present problem, are summarized in Table 4. As in [37] we have introduced the shorthand notation φ k, ,m depending on which integral in (4.94) one considers. Using the map between integrands involving derivatives of the scalar and the corresponding integrals in Table 4, any integral containing zero or two derivatives of the scalar can be directly evaluated. Most integrals containing four derivatives on the scalar can be evaluated directly using this table as well, but there are few cases which require an extension of the results in Table 4 because only certain tensor structures at the four-derivative level were considered in [37]. It is straightforward to generalize these results to any tensor structure with four derivatives on the scalar following the procedure in Appendix A of [37]. However, we will not carry out this generalization here as we will not needed it explicitly.
We can now summarize the complete recursion algorithm. We start by organizing the source terms (4.50) into eigenfunctions R (2k,2 ) of the operator δ B , utilizing the results in Table 3. Taylor expanding these expressions in B i − B oi one obtains the source terms at each order of the Taylor expansion, which are eigenfunctions of δ 0 and δ 0 B . These eigenfunctions are then written in the form where the tensors T I k, contain only derivatives of the scalar φ, but are otherwise independent of φ. Using the identities in Table 4, the integrals in (4.94) can be evaluated to obtain L (2k,2 ) in the form This determines the complete solution of the HJ equation at order k up to linear order in B i − B oi . To obtain the solution at order k +1 we need to evaluate the momenta from the order k solution and substitute them in the source term (4.50) for the order k + 1 equation. We then proceed as before. This procedure is repeated in order to obtain the solution of the HJ equation up to the finite term, where the recursion procedure breaks down. We will discuss when precisely this happens and the significance of the finite part in Section 5.

Solution at order k = 1
In order to illustrate the recursion algorithm we now construct the general solution at order k = 1 and up to order O(B − B o ) in the Taylor expansion. The source term (4.50) for k = 1 and to lowest order in The first step in the algorithm is to decompose this into eigenfunctions of δ 0 B . The last term is an eigenfunction of δ 0 B with eigenvalue −1 and hence it belongs to R 0 (2,2) . This can be deduced by directly evaluating the action of δ 0 B on this term, or by invoking the last entry in Table 1 and noticing that σ i k B l π kl (0,0) = 0. The same result can also be read off the last entry in Table 3. The other three terms are not eigenfunctions of δ 0 B , but they can be decomposed into eigenfunctions of δ 0 B using the projection operator σ i j . For the scalar we have where both terms in this decomposition are eigenfunctions of δ 0 B with respective eigenvalues 1 and −1.
where the first term has δ 0 B eigenvalue 3 and the second 1. However, there cannot be any eigenfunction of δ 0 B with eigenvalue 3 when k = 1 and therefore σ ij σ kl F oik F ojl must vanish identically. Finally, the Ricci scalar can be decomposed into two eigenfunctions of δ 0 B with eigenvalues 1 and −1, but the decomposition is less trivial. Namely the naive decomposition is not correct in this case because these two terms are eigenfunctions of δ 0 B only up to total derivatives. In particular, where the first two eigenfunctions have eigenvalue 1 and the last two −1. Using the decomposition of the Ricci tensor in Table 12 it is easy to see why these particular combinations arise. In terms of anisotropic geometric quantities these become (4.106) which makes it manifest that the eigenfunction with eigenvalue 1 contains only spatial derivatives, while the one with eigenvalue −1 contains only time derivatives.
Next we need to write these terms in the form (4.98) by making explicit all the dependence on the scalar field φ. Since we have where f ij is defined in Table 12 in Appendix C. Hence, which confirms the conclusion we reached above that σ ij σ kl F oik F ojl must vanish identically based on its eigenvalue under δ 0 B . Moreover, and so Finally, Collecting all results, the source term R 0 (2) can be decomposed in terms of a convenient basis of eigenfunctions as described in Table 5, where we also introduce the linear operator The corresponding coefficients of the solutions L 0 (2,0) and L 0 (2,2) of the HJ equation, in the parameterization (4.99), are then obtained using the integration formulas in Table 4, which appear in the last column of Table 5.
Similarly we find that the O(B − B o ) source terms for k = 1 are Decomposing these in spatial and time components leads to the expressions presented in Table 6. In each case, the corresponding solutions of (4.94), obtained using Table 4, are listed in the last column. One must remember, however, that (4.114) do not provide the full source for B oj L 1j (2,2 ) given in (4.90). In particular, the full source for B oj L 1j (2,2 ) contains terms involving the momenta obtained from the O(1) solution in the Taylor expansion. Table 5. General solution of the first recursion relation in (4.94) at order k = 1. The second column from the right describes the source of the inhomogeneous equation in the form (4.98), while the last column gives the solution L 0 (2,0) and L 0 (2,2) in the parameterization (4.99). The shorthand notation used in the last column is defined in (4.96).

Computation of momenta at order k = 1
The general solution of the recursion relations (4.94) at order k = 1 and to the first two orders in the B i − B oi expansion is given in Tables 5 and 6. In order to proceed to the next order in k, we need to compute all the canonical momenta from the solution at order k = 1 by evaluating the corresponding functional derivatives. The identities It is useful to write these momenta entirely in terms of quantities that directly pertain to the geometry of the spatial surfaces and their embedding in the constant radial slices Σ r , rather than covariant variables with respect to Σ r diffeomorphisms, since these variables are best suited to facilitate the decomposition of the inhomogeneous term R 0 (2k) at the next order in k into eigenfunctions of δ 0 B . All these quantities and their geometric meaning is defined in Appendix C, where various useful identities are presented as well. Table 6. General solution of the second and third recursion relations in (4.94) at order k = 1. The second column from the right describes the sources σ i j R 1j (2,2 ) and BojR 1j (2,2 ) of the inhomogeneous equations in the form (4.98), while the last column gives the components σ i j L 1j (2,2 ) and BojL 1j (2,2 ) of the solution in the parameterization (4.99). The shorthand notation used in the last column is defined in (4.96). The results in this table can be extended to the full source Boj R 1j (2,2 ) in (4.90) once the canonical momenta at order O(1) in the Taylor expansion are evaluated.
In terms of the anisotropic variables the momenta following from the O(1) solution in Table 5 are π 0ij (2,0) p 1 (2,0) p 2 (2,0) p 3 (2,0) p 5 (2,2) p 1 (2,2) p 4 (2,2) p 5 The coefficients p I k, appearing in these expressions are given in the last column of Table 5. Finally, the vector momenta do not require functional differentiation since they are given directly by the solution of the last two equations in (4.94). Namely, from (4.77) we have

(4.122)
These expressions can be simplified by noticing that, based on the eigenvalues in Table 3, the following quantities must vanish: Since σ ij is asymptotically positive definite it follows that σ i k B l π kl (2,0) = 0, P i (2,2) = 0. These identities have derived abstractly using the eigenvalues of the derivative and gradation operators, but can be checked explicitly. The first of these identities is is easily seen to hold for the momenta (4.115). The second identity is less obvious at this point, but can be checked in the examples in Section 6.
Finally, using these identities, as well as (4.78) in order to properly isolate the O(1) part of R (4,0) , R (4,2) and R (4,4) , we can write the inhomogeneous terms at order k = 2 in the simpler form (2,2) (2,2) π 0kl (2,2) + 2 n i n j π 0ij , where we have defined and Moreover, the inhomogeneous term (4.90) can be written as Inserting the expressions for the canonical momenta from the order k = 1 solution in these inhomogeneous terms one can use Table 4 in order to obtain the corresponding solutions L 0 (4,0) , L 0 (4,2) and L 0 (4,4) of the recursion relations (4.84).

Asymptotic expansions, Ward identities & the holographic dictionary
So far we have concentrated on the algorithm for obtaining the general asymptotic solution of the radial Hamilton-Jacobi equation with Lifshitz or hyperscaling violating Lifshitz boundary conditions. The purpose of the current section is to point out certain generic features of this solution and to explain its relevance in the context of holography.

General structure of the solution, boundary counterterms & renormalized action
In the previous sections we have shown that this solution takes the form of a graded covariant expansion in simultaneous eigenfunctions of the operators δ and δ B , where each term in this expansion is a functional Taylor expansion in B i − B oi . Schematically, By construction, each term in this expansion has definite asymptotic behavior, which is counted by the dilatation operator, δ D , defined via the leading asymptotic behavior of the operator ∂ r [44]. In order to determine the form of the dilatation operator we need to identify which field components are allowed to have independent sources by the boundary conditions, as well as their asymptotic behavior. As we have seen in the Section 4, Lifshitz boundary conditions are equivalent to the covariant constraint (4.24) and so the covariant fields permitted to have independent sources are the metric γ ij , the scalar φ, and the time component of B i − B oi . More concretely, decomposing B i − B oi in timelike and spacelike components using the projection operator However, (4.24) implies that the source of σ j i B j must vanish for Lifshitz boundary conditions and therefore, since B oi is a function of γ ij and φ, the only independent source in B i − B oi is contained in the scalar field It follows that the dilatation operator can be identified with the asymptotic form of the operator The leading asymptotic form of γ ij and φ can be obtained immediately from (4.2) and (4.28), namelẏ The leading asymptotic behavior of ψ can be inferred from that of Y − Y o in (3.59), but it is instructive to derive it from first principles in the present more general setting. From (C.8) and (C.7) we obtaiṅ Combining this with (4.2) (ignoring transverse derivatives for now) yieldṡ  10) in complete agreement with the result (3.59) we obtained in Section 3. Moreover, (5.9) implies thaṫ ψ ∼ −∆ − ψ, (5.11) and therefore the dilatation operator takes the form Several comments are in order here. Firstly, it is clear from this form of the dilatation operator that every term in the expansion (5.1) has definite asymptotic behavior. Namely, 13) where recall that Secondly, we can now state more precisely why the dilatation operator is in general not a suitable operator in whose eigenfunctions to expand the solution of the HJ equation in the presence of a scalar field φ. Namely, each term in (5.1) is in general only an asymptotic eigenfunction of δ/δφ. However, an expansion in simultaneous eigenfunctions of δ 0 and δ 0 B allows us to determine the φ-dependence in closed form. Finally, note that in the relativistic limit z → 1 which is the operator used in [37] for the corresponding relativistic problem. The definite asymptotic form (5.13) of each term in the expansion (5.1) allows us to determine up to which order in k, and m we need to go. The criterion is that we need to determine all the terms for which When this quantity is positive the corresponding term in (5.1) clearly diverges in the UV and needs to be removed with a local counterterm. The terms for which the inequality is saturated (which can only happen for certain values of the parameters z, θ = −dµξ and ∆ − ) are also divergent, but only linearly in the radial UV cut-off r o . This follows from the fact that a term in the expansion (5.1) corresponding to the integers k, and m has a single factor of C k, + dµξ − m∆ − in the denominator. This can be seen directly from the recursion formulas (4.94). Terms corresponding to integers for which the above inequality is saturated (if there are any) consequently have poles. By the usual dimensional regularization trick [44] where the radial cut-off is defined via the pole is traded for explicit cut-off dependence. Such terms normally give rise to conformal anomalies since the explicit cut-off dependence breaks the invariance of the corresponding term under radial translations. In the absence of a linear dilaton, i.e. when µ = 0, this is the best one can do since there is no regularization scheme where full bulk diffeomorphism invariance is preserved. However, when µ = 0 the cut-off r o can be replaced with φ/µ, thus preserving complete diffeomorphism invariance [37]. The terms for which the above inequality is saturated, therefore, always require regularization but they only lead to conformal anomalies when µ = 0. This makes sense from the dual field theory point of view: for µ = 0 the theory has a running coupling in the UV. Irrespectively of whether there are integers for which the inequality (5.17) is saturated, there is always an independent solution of the HJ equation starting with dilatation weight zero and is therefore UV finite. Namely, the solution (5.1) takes the form where S reg is the lowest order term of this new independent solution and the dots stand for terms of negative dilatation weight that vanish in the UV. S reg satisfies δ D S reg = 0, (5.19) and can be parameterized as where the quantities π ij , π i and π φ correspond undetermined integration functions of the HJ equation, subject only to certain constraints that we will derive shortly. In particular they are not functions of the induced fields γ ij , B i and φ. As we have discussed in Section 3, a solution of the HJ equation that contains as many integration 'constants' as generalized coordinates is a complete integral of the HJ equation, meaning that it is a sufficiently general solution of the HJ equation to describe all solutions of the second order equations of motion. In particular, every solution of the second order equations corresponds to specific values for the integration constants π ij , π i and π φ . On the space of solutions of the equations of motion that have arbitrary sources for the fields γ ij , B i and φ (as allowed by the boundary conditions) and satisfy a certain regularity condition in the IR the quantities π ij , π i and π φ become non-local functionals of the sources. The significance of S reg stems from the fact that the solution, S, of the HJ equation is nothing but the on-shell action. More accurately, for every solution of the equations of motion, the corresponding on-shell action is exactly equal to a complete integral of the HJ equation, for a specific choice of the integration functions π ij , π i and π φ . The AdS/CFT dictionary identifies the on-shell action, and hence the complete integral S, with the generating function of connected correlation functions. The on-shell action is UV divergent, but its identification with the asymptotic complete integral (5.18) means that these UV divergences can removed by the local covariant counterterms defined by This means that S reg = S + S ct is identified with the regularized on-shell action, and therefore (by the AdS/CFT dictionary) with the regularized generating function of connected correlation functions. The renormalized on-shell action, or generating function, is given by the limit  More importantly, the flow equations allow us to identify generically the complete set of modes parameterizing the symplectic space of asymptotic solutions, without deriving the full form of these solutions. We have already identified a set of integration constants that parameterize S reg in the asymptotic complete integral (5.18) of the HJ equation. These integration constants enter in the flow equations aṡ n ∼ e zr n (0) (x), n a ∼ e 2r n (0)a (x), σ ab ∼ e 2r σ (0)ab (x), where n (0) (x), n (0)a (x), σ (0)ab (x), ω (0) (x), φ (0) (x) and ψ − (x) are arbitrary functions of the transverse coordinates, and the given asymptotic form of φ is valid for µ = 0. For µ = 0 the asymptotic form of φ depends on the subleading terms in the potentials that define the bulk theory. Note that the asymptotic behavior of the gauge field A i is completely determined in terms of these fields and does not contain any additional source allowed by the asymptotic Lifshitz condition (4.24), namely 10 The source ω (0) (x), therefore, corresponds to a pure gauge transformation. The radial dependence of the sources (5.25) allows us to determine the radial dependence of the modes π ij , π i and π φ parameterizing S reg . Since the only fields with independent sources are those in (5.25), (C.7) and (C.8) imply that so that This motivates us to define the following quantities: Note that the quantity σ i j π j couples to variations of B i orthogonal to B oi and hence it corresponds to the 1-point function of an irrelevant operator. Although Lifshitz boundary conditions do not allow for a source of this operator it can have a non-zero expectation value. In terms of these variables the general variation of S reg with respect to the sources becomes where δγ ij T ij = −2nδn T tt + 2δn a ( T ta + n a T tt ) + δσ ab ( T ab − n a n b T tt ). (5.31) The integration functions defined in (5.29) are the symplectic conjugate variables to the sources (5.25) (except for E i whose source is set to zero) and, therefore, they are identified via the holographic dictionary with the renormalized 1-point functions of the dual operators. The asymptotic form of these 1-point functions follows from the asymptotic form of the sources (5.25), together with the fact that S reg has dilatation weight zero. Namely, 1-point function source spatial stress tensor As we shall confirm shortly by deriving the Ward identities these modes satisfy, this is precisely the spectrum of the energy-momentum complex [29], plus the two additional scalar operators O φ (x) and O ψ (x). Note that the asymptotic form of the momentum density and the energy flux differ by a factor of e −r relative to the operators defined in [29], which reflects the fact that the indices of the corresponding operators in that reference are frame indices and not spacetime indices. The operators in [29] can be obtained by contracting our P i and E i with a spatial vielbein. However, the operators that enter the covariant Ward identities are P i and E i and not the ones with frame indices. Inverting the relations (5.29) and inserting the asymptotic behaviors (5.32) in (5.24) we obtain the dependence of the asymptotic expansions on the normalizable modes.

Holographic Ward identities
The holographic Ward identities follow directly from the first class constraints (2.16). The Hamiltonian constraint leads to the trace Ward identity, while the momentum and gauge constraints imply the anisotropic diffeomorphism Ward identities. However, the trace Ward identity can be derived much more easily from the invariance of the HJ solution under radial translations.

Diffeomorphism Ward identity
Combining the momentum and gauge constraints in (2.16) and applying them to S reg gives The leading asymptotic form of the vector field, B i ∼ (1 + ψ)B oi , implies that where we have assumed that ∆ − > 0 in the second step. The above constraint then takes the form Using the variables introduced in (5.29) and (5.32) we obtain the constraint Different components of this equation behave differently asymptotically. Isolating components with the same scaling behavior using the projection operator σ i j we arrive at the three anisotropic Ward identities with arbitrary sources When all sources are set to their background value for flat space these identities reduce to the Ward identities for the energy-momentum complex discussed in [29], plus conservation of the momentum density.

Trace Ward identity
The trace Ward identity can be derived by considering the transformation of S reg under an infinitesimal local radial translation r o → r o +δσ(x), which induces an anisotropic Weyl transformation on the boundary. Such a tranformation in general gives If there is no explicit dependence on the radial cut-off in the counterterms, this variation must vanish identically. If, however, there is an explicit dependence on the radial cut-off, then the counterterms are not invariant and hence there is an additional contribution from the coefficients of the radial cut-off in the counterterms, i.e. the conformal anomaly. In particular, where the conformal anomaly is given by As we pointed out earlier, there in no conformal anomaly when µ = 0 since in that case there is a regularization scheme that does not break radial translations.

Examples
In order to appreciate how the algorithm for solving the HJ equation recursively works in practice it is instructive to work through a few examples.
as well as for σ i j π 0j (2,2 ) : Differentiating the expressions in (6.6) with respect to the metric γ ij leads to the momenta (4.115) and (4.116), which now take the form π 0ij (2,0) p 1 (2,0) p 3 Using these expressions we obtain the solution to the remaining third recursion relation in (6.2): As an illustration let us consider the case d = z = 2 which has been discussed before e.g. in [33]. From (3.61) follows that in this case ∆ − = 0 and hence Y − Y o ∼ r as r → ∞ and so we must set this mode to zero to ensure asymptotically locally Lif boundary conditions [27]. The zero order solution of the Taylor expansion in B − B o therefore gives the full solution in this case. The terms that contribute to the UV divergences, therefore are S =ˆd d+1 x L 0 (0) + L 0 (2,0) + L 0 (2,2) + L 0 (4,0) , (6.19) where L 0 (0) was given in (3.43). The terms L 0 (2,2) and L 0 (4,0) have poles at d = z = 2 and therefore both contribute to the conformal anomaly. Setting z = 2 and where r o is the UV cut-off, these terms become where denotes equivalence up to total derivative terms and we have used the identities (see Table 12) D [i q j] = 0, D k q k = D k q k + q k q k , and R ij = 1 2 Rσ ij for d = 2. Using the fact that up to total derivative terms it is easy to check that L 0 (4,0) vanishes identically in agreement with [33].

Exponential potentials with µ = 0
A second interesting example is a generalization of the Einstein-Proca theory discussed above obtained by relaxing the condition that the scalar be constant and that ξ = 0. In particular, the scalar is not necessarily constant in this case and the potentials defining the Lagrangian take the form where The first three coefficients in the Taylor expansion of the superpotential correspondingly take the form where ∆ − now must be evaluated using the general expression (3.60) instead of (3.61). Note that as for the Einstein-Proca theory u 0 + Z Z u 1 = 0, (6.26) and therefore the recursion relations that determine the HJ solution are still algebraic and in fact identical to those of the Einstein-Proca theory given in (6.2). The source term (4.128) of the third recursion relation now takes the form whereū 2 := e ξφ u 2 , and from (4.127) we obtain The solutions to the three recursion relations (6.2) can therefore be written in the form   , (6.29) and again we have used (4.77). Note that in the limit ξ → 0 these expressions reduce to the corresponding ones in (6.5) for the Einstein-Proca theory. The source terms R 0 (2k,2 ) , σ i j R 1j (2k,2 ) and B oj R 1j (2k,2 ) for k = 1 are given in Tables 8 and 9.
Solution at order k = 1 The solution (6.29) at order k = 1 can be read off the last column of Tables 8 and 9. Namely, from Table 8 we see that the solution for L 0 (2,2 ) is: (6.30) Moreover, Table 9 gives for σ i j π 0j (2,2 ) : Table 8. General solution of the first recursion relation in (4.94) at order k = 1 for exponential potentials and µ = 0. The second column from the right describes the source of the inhomogeneous equation in the form (4.98), while the last column gives the solution L 0 (2,0) and L 0 (2,2) in the parameterization (4.99). Table 9. General solution of the second and third recursion relations in (4.94) at order k = 1 for exponential potentials and µ = 0. The second column from the right describes the sources σ i j R 1j (2,2 ) and BojR 1j (2,2 ) of the inhomogeneous equations in the form (4.98), while the last column gives the components σ i j L 1j (2,2 ) and BojL 1j (2,2 ) of the solution in the parameterization (4.99). Moreover, we have definedū2 := e ξφ u2. The results in this table can be extended to the full source Boj R 1j (2,2 ) in (4.90) once the canonical momenta at order O(1) in the Taylor expansion are evaluated.
As for the Einstein-Proca theory in the previous example, the solutions to the recursion relations (6.2) we obtained above determine the solution of the HJ equation up to and including order k = 1 and O(B − B o ). These suffice in order to determine the solution of the HJ equation at order k = 2 but only to order O(1) in the Taylor expansion in B − B o , corresponding to the solution of only the first recursion relation in (6.2) for k = 2. Again we will not write these solutions explicitly since they are too lengthy. But they can be evaluated straightforwardly with Mathematica by inserting the k = 1 results above into (4.125).
This example can be compared directly with the model discussed in [36], which corresponds to the following values of our parameters: Moreover, the two scalars in [36] are related to the scalar φ here as with φ → 0 in the UV. Dropping terms with derivatives on the scalar φ in this case we get the same result for L 0 (0) , L 0 (2,0) and L 0 (2,2) as in (6.21), but for L 0 (4,0) we now get where again denotes equality up to total derivative terms. This quantity is the only non-trivial conformal invariant with four spatial derivatives in d = 2 and for z = 2 [33]. Note that this model is related to the Einstein-Proca theory of the previous example only by a change of frame since ξ = 1/2 here. So the effect of going from the Einstein frame (where no purely spatial anomaly is generated) to a non-Einstein frame is to generate a non-zero coefficient for this conformal invariant in the anomaly. However, the expression for the anomaly given in [36] does not agree with our result. Namely, in our notation the purely spatial part of the expression in [36] is which is in fact not a conformal invariant. We have traced the discrepancy to the fact that the O(B − B o ) contribution to the 2-derivative momenta has not be taken into account in [36].
6.3 Exponential potentials with µ = 0 As a final example we consider a model with exponential potentials V ξ = V o , Z ξ = Z o e −2(ξ+ν)φ , W ξ = W = W o e −2(ξ+ν)φ , (6.43) corresponding to the first three superpotential coefficients (∆ − is again given by (3.60)) u 0 (φ) = (z − 1 + d(1 + µξ)) e −ξφ , u 1 (φ) = 1 2 (z − 1)e −ξφ , u 2 (φ) = 1 but without any restriction on the parameters that define the boundary conditions. In particular, the crucial difference in this example relative to the previous two is that µ = 0 and so the recursion relations (4.84), (4.85) and (4.86) are no longer algebraic. However, there is still some simplification due to the fact that the potentials are exactly -not merely asymptotically -exponentials. The inhomogeneous solutions (4.94) become where we have used the fact that for the present example Using this, together with A k, = 0, we see that the integrals in Table (4) reduce in this case to ordinary integrals over the exponential coefficients of any tensor structure involving derivatives on the scalar. In fact, since the overall exponential function of the scalar in the source terms in (4.94) is easily determined to be R 0 (2k,2 ) ∼ e dξφ , σ i j R 1j (2k,2 ) ∼ e (dµξ+z− )φ/µ , B oj R 1j (2k,2 ) ∼ e dξφ , (6.48) we can perform the integrals over the scalar fields generically without any reference to the explicit form of these source terms. The source term (4.128) of the third recursion relation in (4.94) can be written as while from (4.127) we get Q 0 (2k,2 ) := π 0 (2k,2 ) + dn i n j π 0ij (2k,2 ) − dν 2α Φ 0 (2k,2 ) , P 0 (2k,2 ) := B ok π 0k (2k,2 ) −  Table 10. General solution of the first recursion relation in (4.94) at order k = 1 for exponential potentials and µ = 0. The second column from the right describes the source of the inhomogeneous equation in the form (4.98), while the last column gives the solution L 0 (2,0) and L 0 (2,2) in the parameterization (4.99).
where ζ defined in (4.126) now becomes Performing the integrations over the scalar field in (4.94) we arrive at the solutions L 0 (2k,2 ) = − 1 C k, + dµξ R 0 (2k,2 ) , 2 ) [γ(x), φ(x); x ] − 2(z − 1)n k n l π 0kl (2k,2 ) (6.52) where again we have used (4.77). In the limit µ → 0 these expressions reduce to the corresponding ones in (6.29) of the previous example. The source terms R 0 (2k,2 ) , σ i j R 1j (2k,2 ) and B oj R 1j (2k,2 ) for k = 1 are given in Tables 10 and 11. Note that since the hyperscaling parameter θ in the Einstein frame is given by the combination −dµξ, we see that the denominators in these recursion relations are shifted by θ relative to the previous examples.   Table 11. General solution of the second and third recursion relations in (4.94) at order k = 1 for exponential potentials and µ = 0. The second column from the right describes the source terms σ i j R 1j (2,2 ) and BojR 1j (2,2 ) of the inhomogeneous equations in the form (4.98), while the last column gives the components σ i j L 1j (2,2 ) and BojL 1j (2,2 ) of the solution in the parameterization (4.99). The constantū2 ≡ e ξφ u2 has been introduced to simplify the expressions. The results in this table can be extended to the full source Boj R 1j (2,2 ) in (4.90) once the canonical momenta at order O(1) in the Taylor expansion are evaluated. and π 0 (2,2) = explicitly here. However, the results we have presented allow one to evaluate P 0 (2,0) and P 0 (2,2) easily by evaluating the last expression in (6.52) using Mathematica. The same holds for the solution at k = 2 and O(B − B o ) 0 , which can be obtained by inserting the k = 1 results in (4.125).

Concluding remarks
In this paper we have developed a general algorithm for constructing the holographic dictionary for a large class of theories that admit asymptotically locally Lifshitz and hyperscaling violating Lifshitz boundary conditions with arbitrary dynamical exponents. This dictionary only exists for θ ≤ d + z, z ≥ 1, since there are no well defined asymptotic expansions for θ > d + z and z ≤ 1.
The algorithm we developed relies entirely on the metric formulation of the dynamics and there is no need for the introduction of vielbeins at any point. The objective of the algorithm is the systematic construction of the most general asymptotic solution of the radial Hamilton-Jacobi equation subject to asymptotically locally Lifshitz and hyperscaling violating Lifshitz boundary conditions. This is achieved by expanding the solution of the Hamilton-Jacobi equation in simultaneous eigenfunctions of two commuting functional operators, which generalizes the standard expansion in eigenfunctions of the dilatation operator to non-relativistic and non-scale invariant boundary conditions. The resulting recursive procedure does not require any ansatz and it is entirely algorithmic. In future work we hope we will be able to implement this algorithm in a symbolic computation package.
The entire holographic dictionary can be derived from this asymptotic solution of the Hamilton-Jacobi equation as is shown in Section 5. In particular, the asymptotic Fefferman-Graham expansions, including the sources and 1-point functions, are derived directly from this asymptotic solution of the Hamilton-Jacobi equation, without any need for solving the second order equations of motion. In fact, the Hamilton-Jacobi equation leads to a much more efficient method for computing renormalized correlation functions as well [52,61,62]. Our method provides a solid basis for computing correlation functions in asymptotically Lifshitz and hyperscaling violating Lifshitz backgrounds, and we intend to explore this direction in future work. Another potential application of the present work is in the holographic computation of entanglement entropy.
Finally, we have shown that the unique non-trivial conformal invariant for z = 2 in 2 dimensions with four spatial derivatives appears in the conformal anomaly of an Einstein-Proca theory, provided the latter is coupled with a dilaton and one moves away from the Einstein frame. To our knowledge, this is the first example where this term is actually generated, implying that the detailed balance condition does not hold in this case [33]. More generally, the algorithm presented here provides a systematic tool for generating non-relativistic conformal invariants for any dimension and any value of the dynamical exponent z ≥ 1.