A 3+1 formulation of the 1/c expansion of General Relativity

Expanding General Relativity in the inverse speed of light, 1/c, leads to a nonrelativistic gravitational theory that extends the Post-Newtonian expansion by the inclusion of additional strong gravitational potentials. This theory has a fully covariant formulation in the language of Newton-Cartan geometry but we revisit it here in a 3+1 formulation. The appropriate 3+1 formulation of General Relativity is one first described by Kol and Smolkin (KS), rather than the better known Arnowitt-Deser-Misner (ADM) formalism. As we review, the KS formulation is dual to the ADM formulation in that the role of tangent and co-tangent spaces get interchanged. In this 3+1 formulation the 1/c expansion can be performed in a more systematic and efficient fashion, something we use to extend the computation of the effective Lagrangian beyond what was previously achieved and to make a number of new all order observations.


Introduction
The geometric formulation of nonrelativistic physics has experienced a resurgence of interest in the last few years, see e.g. [1][2][3][4] for a recent overview of the relevant geometries, symmetries and some applications. Formulating the nonrelativistic expansion of general relativity in a geometric fashion has the advantage that it keeps, by construction, general coordinate invariance manifest. Building on [5] this idea was applied to the Post-Newtonian (PN) expansion in [6], reformulating it in terms of Newton-Cartan geometry. The PN expansion, see e.g. [7] for an introduction, is an expansion around flat space-time and is for this reason a double expansion: nonrelativistic, i.e. an expansion in the inverse speed of light 1/c, as well as weakly coupled, i.e. an expansion in Newton's constant G N . From the geometric perspective the assumption of weak coupling was implemented in [5,6] as an assumption on the nonrelativistic connection. In [8] this assumption was relaxed, which extends the nonrelativistic theory to so called twistless torsional Newton-Cartan geometry. Physically speaking, this extended theory adds to Newtonian gravity an additional gravitational potential that describes strong gravitational time dilation effects. Rather than an expansion around flat space-time, it is an expansion, in powers of 1/c 2 , around an arbitrary static space-time [9] and as such extends the PN expansion. One of the advantages of the 1/c 2 expansion is that an expansion of the Einstein equations goes hand in hand with an expansion of the Einstein-Hilbert Lagrangian [10,11]. This makes it possible to formulate the expansion -truncated up to a given order -as a self-consistent nonrelativistic theory with a fully geometric/covariant action and associated variational principle.
Despite the elegance of the geometrically formulated 1/c 2 expansion [8,11], manifest space-time coordinate invariance comes with a cost, in that it introduces a number of fields which are pure gauge, and this increases the complexity of the expansion. Some work has been done to streamline this [12][13][14] or to organize it based on symmetry considerations [10,[15][16][17], but this has so far not led to higher order results.
In this paper we revisit the 1/c and 1/c 2 expansions of general relativity in a 3+1 formulation, where an explicit choice of time coordinate is made. Although the theory loses some of its elegance in such formalism, it keeps some key features: manifest invariance under spatial coordinate transformations as well as an action with well defined variational principle. Furthermore, it has the advantage that it makes the physical degrees of freedom more explicit and the expansion more transparent. This allowed us first of all to carry out the expansion to higher order than before. It also reveals that the 1/c expansion is equivalently an expansion in time derivatives around a quasi-stationary relativistic background, something which was pointed out at leading order in [18]. Also the relation between the 1/c expansion and the 1/c 2 expansion, a consistent truncation to even orders, becomes more clear so that we could use it to devise an algorithm that constructs the 1/c 2 expansion to power c −2N from the 1/c expansion to power c −N .
A 3+1 formulation of Newton-Cartan theory was recently considered in [19,20], for use in Newtonian cosmology. Rather than introducing the 3+1 formulation at the nonrelativistic level, we will start from General Relativity (GR) in 3+1 form and then expand. The most familiar 3+1 formulation of GR is that of Arnowitt, Deser and Misner (ADM) [21], see e.g. [22] for a detailed introduction. In ADM form the 4d relativistic metric is where the fields (N, N i , h ij ) depend on time t as well as the spatial coordinates x i . Mathematically speaking h ij is the induced metric on a spatial hypersurface with tangent vectors ∂ i . The basis in tangent space is completed by the vector , which is normal to the hypersurface. The factors of the speed of light c are naturally there for dimensional reasons. Now observe that the ADM decomposition becomes singular in the c → ∞ limit, since u is then no longer linearly independent of the ∂ i . This implies the ADM decomposition is not a suitable starting point for a 1/c expansion. (The same feature makes it 1 a good starting point for the ultra relativistic expansion in small c, see e.g. [23][24][25]). A 3+1 decomposition suitable for the 1/c expansion was introduced by Kol and Smolkin (KS) in [26]. In KS form the 4d relativistic metric reads where again the fields (M, C i , h ij ) depend on t and x i . This decomposition is based on a preferred basis for the cotangent space. The basis dx i for the spatial subspace gets extended to a basis of the total 4 dimensional cotangent space by the one-form n = M(cdt + C i dx i ). One sees that this basis remains generating in the limit c → ∞ since n remains linearly independent of the dx i . This indicates that the KS formulation is the appropriate 3+1 decomposition for the nonrelativistic expansion. Kol and Smolkin originally introduced the decomposition (1.2) in [27] in the context of the PN expansion. In that case, where one assumes M to be of the weak field form M = 1 − 2ϕ c 2 + O(c −4 ), the KS ansatz, also known as the Kaluza-Klein decomposition, greatly simplifies (part of) the PN expansion, see e.g. [28] for a review.
In this work, where we consider the 1/c expansion and make no weak field assumption, we will allow the fields (M, C i , h ij ) to be arbitrary analytic functions in 1/c, which in particular means that M = (0) M + (1) M c −1 + (2) M c −2 + O(c −4 ), with (k) M arbitrary functions of t, x i (only constrained by the expanded Einstein equations).
While we can relate (2) M to the Newtonian potential ϕ, M and (1) M are additional potentials not present in the PN expansion, capturing strong nonrelativistic gravitational effects. A similar observation holds for the expansion of C i and h ij . The precise relation to the PN expansion (up to 1PN order) is worked out in section 5.2, see also table 2.
Our paper is organized as follows. In section 2 we review both the ADM and KS decomposition of GR, we do so in a way that allows to treat both decompositions simultaneously and which emphasizes their dual nature. In section 3 we focus on the KS formulation and discuss the Einstein equations as well as the diffeomorphism symmetry and the associated Noether/Bianchi identities from the point of view of this 3+1 decomposition. At the end of the section, in subsections 3.3 and 3.4 we motivate and perform the convenient field redefinition M = e ψ 2 , h ij = e −ψ γ ij and make factors of c and time derivatives explicit. In section 4 we then come to the 1/c expansion. The new formalism allows us to make a number of new observations: we show how all subleading equations take the form of a sourced linear second order PDE and compute the linear differential operator that is the same at all orders. We also discuss how by a convenient choice of gauge all the subleading coefficients in the expansion of γ ij can be made trace-less. We then compute the Lagrangian of the expansion up to NNLO, i.e. up to order c −2 . Previously only the leading order, i.e. order c 0 , had been computed [18]. In section 5 we recall how the 1/c 2 expansion is a consistent truncation of the 1/c expansion to even orders. Here the use of the 3+1 formulation makes the connection to the 1/c expansion more transparent and allows us to derive an algorithm, that we call the shuffling algorithm, with which the 1/c 2 expansion up to order c −2N can be constructed from the 1/c expansion up to order c −N . At NLO, i.e. order c −2 , we recover the results of [8,11] and spell out the precise relation. The shuffling algorithm then allows us to compute the Lagrangian up to NNLO, i.e. order c −4 , for the first time. This is the order where the 1PN correction to Newtonian gravity finds itself and we show how indeed, upon setting to zero a number of additional potentials, the 1/c 2 expansion reproduces the well-known 1PN equations. We conclude in section 6 where we also point towards possible future directions. Finally there are three appendices. In appendix A we spell out two technical derivations of results used in the main text, in appendix B we provide some extra terms present outside the traceless gauge and the (un-expanded) equations of motion with powers of c and space-like time-like Table 1: We'll refer to a coordinate system (t, x i ) as a space-time split or d + 1 decomposition if it falls into either the ADM or KS class (the GN class is the intersection of both). They are defined by the properties of the coordinate bases of tangent and co-tangent spaces, as listed in this table.
time derivatives made explicit are written down in appendix C.

Two dual 3+1 formulations of GR
Let M be a d + 1 dimensional manifold equipped with a Lorentzian metric g µν . We'll denote the associated Levi-Civita connection with ∇.

Choice of coordinates
Intuitively a (local) d + 1 formulation would amount to a (local) choice of coordinates (x µ ) = (t, x i ) with t a 'time coordinate' and the x i 'spatial coordinates'. We could try to make this more precise by requiring ∂ t to be a time-like vector and ∂ i to be space-like vectors. Although imposing both these conditions is always possible, e.g. in Gaussian normal (GN) coordinates, this would restrict the allowed set of coordinates very much. A weaker choice, which we'll call Arnowitt-Deser-Misner (ADM) coordinates [21], only requires the ∂ i to be spatial. Let us however point out that this automatically implies that the one-form dt will be time-like. So in this case we see that it is a mix of the coordinate basis of tangent vectors and that of co-tangent vectors which is required to be space/time-like respectively. This then suggests a naturally dual alternative, which is to choose coordinates such that the dx i are space-like which then implies that ∂ t is time-like. We'll call such coordinates Kol-Smolkin (KS) coordinates [26]. For a summary of the definitions see table 2.1.
Elementary example Before discussing these classes in full generality in the next subsection, let us illustrate them in a simple example. Consider the Minkowski metric in inertial coordinates One easily verifies that (t, x i ) satisfy our definition of GN coordinates. If we define new coordinates (t,x i ), witht = t andx i = x i − v i t for some constants v i , then the metric takes the form While∂ i remains spatial,∂ t is no longer time-like when v i v i ≥ 1. But sincê g tt = −1 independent of v i , it follows that dt remains time-like. The coordinates (t,x i ) are thus an example of ADM coordinates. Alternatively we can work in coordinates (t,x i ), witht = t − w i x i andx i = x i for some constants w i . In these coordinates the metric reads One verifies that∂ t is always time-like, while∂ i is not space-like when w 2 i ≥ 1. But the dx i are always space-like and so (t,x i ) are an example of KS coordinates for Minkowski space, independent of the values of w i .
Note that indeed the ADM and KS classes are not mutually exclusive, their intersection is exactly the GN class.

Decomposition of the metric
Although our aim is to decompose the space-time metric g µν in either ADM or KS coordinates, it will be interesting to introduce a formalism that can treat both at the same time. This will then also allow us to decompose the Einstein-Hilbert action simultaneously for both cases in section 2.4. The decomposition of the metric we'll perform here will also shed light on the definitions of the previous subsection, which might have appeared a bit ad-hoc there.
We start by choosing 2 a time-like vector field u µ , which without loss of generality we can assume to be normalized: g µν u µ u ν = −1. Of course this is equivalent to the choice of a time-like one-form n µ = g µν u ν , which is again normalized. 2 Mathematically speaking a Lorentzian metric introduces an O(1, d)-structure on M. The addition of the vector field u µ refines this to an O(d)-structure, also known as an Aristotelian structure, which is a combination of compatible Galilean and Carrollian structures, see e.g. [29]. The tensors defining this Aristotelian structure are u µ , n µ = g µν u ν , h µν = g µν + u µ u ν and h µν = g µν + n µ n ν .
We can complete the time-like vector into a tangent frame (u, e i ), with each of the e i orthogonal to u: g µν u µ e ν i = 0. The relations e i µ e µ j = δ i j , u µ e i µ = 0 then provide a dual frame (n, e i ). In summary this amounts to Note that we do not require the frame e i to be orthonormal, instead we define It follows that h ij = g µν e i µ e j ν is the inverse of h ij . In terms of this frame and co-frame the metric and its inverse decompose as As we will now point out, there is a natural choice of frames of the form above associated to both ADM and KS coordinates.

ADM decomposition
In ADM coordinates, where ∂ i is spatial and dt time-like, a natural choice of frames satisfying (2.4) is Then, the metric (2.6) takes the well-known ADM form Note that the inverse metric takes the form (2.10)

KS decomposition
In KS coordinates it is dx i which is spatial while ∂ t is time-like, so in this case the natural choice of frames satisfying (2.4) is So the metric (2.6) takes the following form in KS coordinates: while the inverse metric is (2.14) Comparing the forms (2.9, 2.10) with (2.13, 2.14) one explicitly sees how they are dual, in the sense that the role of metric and inverse metric -i.e. tangent space and co-tangent space -get interchanged. We should emphasize that the spatial metric h ij in the KS form of the metric (2.13) is different from the spatial metric h ij in the ADM form of the metric (2.9). Indeed, in case the coordinates (t, x i ) are in the GN class, i.e. of both ADM and KS type, they are related as It is h ADM ij which has the geometrical interpretation of the pull-back of the Lorentzian metric g µν on the constant t hyper-surfaces. Although h KS ij does not have this natural geometric interpretation it is, by construction, a well-defined Riemannian metric on the constant t hyper-surfaces when (t, x i ) are in the KS class. Only if (t, x i ) are in the ADM class will the constant t hypersurfaces be spatial so that then h ADM ij is Riemannian.

Decomposition of the connection
Given a choice of frame and dual frame that satisfy (2.4) one has the following identities Our results below hold for more general frames than (2.7, 2.8) or (2.11, 2.12), as long as we assume them to satisfy the additional condition (2.20) in addition to (2.4).
The identities listed above are then sufficient to decompose the Levi-Civita connection of g µν in the independent components 3 The key point of introducing these objects is that, via (2. 16-2.19), they are sufficient to decompose the covariant derivatives of the frames. In turn that allows to decompose covariant derivatives of tensors as It is a bit tedious but rather obvious to see how this generalizes to tensors with any number of upper and lower 'spatial' indices, e.g. ∇ µ (e i ν e ρ j V i j ), so we will not write out these expressions explicitly.
In (2.25, 2.26) we introduced a modified covariant derivativeD t along t and a modified covariant derivativeD i along the x i . Their precise definitions on a 'spatial' tensor T i 1 ...im j 1 ...jn (t, x) arê Note thatD i differs from the standard covariant derivative D i with respect to the Levi-Civita connection of h ij , i.e.
Note that when combined with (2.17) the assumption (2.20) results in −e µ i ∇ ν e j µ = e j µ ∇ ν e µ i = e j µ Γ µ ρν e ρ i . This in turn guarantees thatΓ k ij is symmetric in ij since it can be written asΓ k ij = e µ i e ν j e k ρ Γ ρ µν .
The difference is two-fold, firstly inD i the standard partial derivative ∂ i is replaced by e µ i ∂ µ and secondly the 'spatial connection'Γ k ij is used instead of Γ k ij . Both these differences vanish in the ADM case, but not in the KS case -as we'll see below.
It is important to point out that the modified covariant derivative (2.28) remains compatible with h ij :D This is not so for the time-like covariant derivative, but one finds the elegant expressionD Let us remark that similarly to the decomposition (2.25) one computes that and finally we also point out the relationŝ (2.36)

The ADM case
When choosing the frames to be of the ADM type (2.7, 2.8) the hatted partial derivatives take the form The hatted time derivative, which in this case is hypersurface orthogonal, includes the well known shift by the vector field N i (see e.g. [22]), while the hatted spatial derivatives are simply partial derivatives. The connection components (2.21-2.24) become in this case First one recognizes in K ij the extrinsic curvature of the induced metric h ij , while Ω i is the so called Eulerian acceleration vector (see e.g. [22]). Secondly one sees that the connectionΓ k ij in this case equals the Levi-Civita connection (2.31), so thatD i = D i . Finally the object ∆ j i plays a role as a 'connection' in the time-like covariant derivative (2.27). Writing that out explicitly one gets the expression 42) where L N is the Lie derivative with respect to the vector field N i . Remark that this also allows one to rewrite (2.38) as

The KS case
Now we assume the frames to take the KS form (2.11, 2.12). The hatted partial derivatives now are∂ The situation is exactly opposite -or dual -to the ADM case (2.37), in that the hatted time derivative is simply (up to a prefactor) the partial derivative, while it is the spatial hatted derivatives that get a shift, this time by the one-form C i . When writing out the objects (2.21-2.24) one finds they are naturally expressed in terms of these hatted partial derivatives 4 : It might appear as if the introduction of the hatted partial derivatives is simply a notation to simplify the expressions (2. 45-2.48). But, as we will discuss in section 3, these hatted derivatives are the natural ones invariant under redefinitions of the time coordinate, which is why all invariant quantities involving derivatives can be expressed in terms of them. Let us now come to the expressions (2.45-2.48) themselves. While the symmetric part of K ij is analogous to the ADM expression (2.43) (see also (2.49)), it no longer has an interpretation as extrinsic curvature. More importantly, K ij now also has a non-vanishing anti-symmetric part (absent in the ADM case), which takes the form of a (generalized) curvature for the 1-form potential C i . The vector Ω i can also in this KS case be interpreted as an 'acceleration vector', in that it provides a relativistic version of Newton's gravitational field vector (see sections 3.3). Since ∆ i j vanishes we can conclude that On the contrary,Γ k ij is in this case not the Levi-Civita connection, rather As mentioned before, this implies that in the KS case, the spatial covariant derivative (2.28) is not the standard one, in that∂ i has a shift proportional to ∂ t , see (2.44), and the connection coefficients have the extra contribution ∆ k ij . When discussing general relativity in KS form the expression (2.28), together with (2.44) and (2.48) is most convenient to work with.

Decomposition of the Einstein-Hilbert action
To rewrite the Einstein-Hilbert action in d+1 form one needs to decompose the Ricci scalar. This can be done by first expressing the Ricci tensor as i )e i µ and re-expressing the covariant derivatives in terms of (2.21-2.24). Upon contraction with g µν one finds 5 Via (2.34) one recognizes the term in between brackets in (2.51) as a total derivative so that 6 It is gratifying to see that the Einstein-Hilbert action takes a simple and universal form for both ADM and KS type decompositions.

The ADM Lagrangian
The ADM Lagrangian [21] is well-known: Indeed it equals (2.53) when one takes the frames of the ADM form (2.7, 2.8), since then n t = N and -as discussed in the previous subsection -∂ i = ∂ i ,Γ k ij = Γ k ij , so thatR = R, the Ricci scalar of the Levi-Civita connection of h ij . Furthermore K ij has the standard ADM form (2.38).

The KS Lagrangian
Since for the KS frames (2.11, 2.12) n t = M we recover from (2.53) the KS Lagrangian [26]: Although this expression might appear almost identical to the ADM Lagrangian (2.54), it is rather different. This is since now, via (2.44) and (2.50),R has a number of additional contributions apart from R, the Ricci tensor of the Levi-Civita connection of h ij . Secondly the tensor K ij is also different from its ADM analog, see (2.45) versus (2.38), in particular it contains a non-trivial anti-symmetric part. When we write out the time derivatives explicitly in section 3.4, the Lagrangian (2.55) becomes (3.45, 3.57) and one sees some of the complexity that is elegantly packaged in the form (2.55).

Dynamics and symmetries in the KS formulation
General relativity in the ADM formulation has been widely discussed. This is much less so for the KS formulation, so in this section we shortly discuss the Einstein equations in the KS formulation, as well as the diffeomorphism invariance of the theory.

Einstein equations
The Einstein equations in KS form can be derived in two parallel and equivalent ways. The first would be to decompose the d+1-covariant Einstein equations using the frames (2.11, 2.12), the second to vary the KS Lagrangian with respect to the fields h ij , C i , M. Let us make the link between these two approaches explicit. The discussion in section 2.4 showed that Via the KS decomposition of the metric (2.13) and the definition of the KS frames (2.11, 2.12) it follows that This implies that if we define one gets the relations δgµν is the Einstein tensor of the metric g µν . In summary we see that indeed there are two ways to compute the G's: either via (3.4), or via (3.3) by a variation of L KS . One checks that both calculations match, with the result In the aboveĜ ij =R ij − 1 2 h ijR . The vacuum Einstein equations are thus equivalent to putting (3.5-3.7) to zero. It is interesting to point out that M appears algebraically in L KS (similarly to N in the ADM formulation), so that G 0 = 0 has the interpretation of a constraint. Another consequence is that the equations do not contain a second time derivative of M (while the second time derivative of both h ij and C i does appear). One obtains an equation for the second spatial derivatives of M, by taking the trace of G ij , we will discuss that equation a bit more in section 3.3.
One can couple matter to the theory by introducing a matter Lagrangian L mat . This then allows to define the KS energy-momentum tensors T : The Einstein equations then take the form Apart from computing the T 's by a variation of the matter Lagrangian wrt M, C i and h ij , one can also obtain them via a decomposition of the usual energy momentum tensor T µν = 1 √ −g δLmat δgµν . Via (3.2) one finds (3.10) By comparing (3.10) to (3.4) the equations (3.9) then follow immediately. One advantage of the form (3.10) is that it can also be used in cases where the matter has no Lagrangian description and is only defined in terms of an energy momentum tensor T µν . Clearly the conservation of energy momentum ∇ µ T µν = 0 will become equivalent to some equations for the T 's. These equations, and their origin in diffeomorphism invariance, will be discussed in the next subsection.

Diffeomorphism symmetries and Noether identities
Because a choice of KS coordinates amounts to a choice of preferred time coordinate, it breaks manifest diffeomorphism invariance. Since they leave the time coordinate invariant, spatial diffeomorphisms, i.e.x i (x j ), do remain manifest. But more general coordinate transformations map one choice of KS coordinates into another; this will leave the theory invariant, but in a less manifest way. One reason to have a closer look at diffeomorphism invariance in the KS formulation is that via the associated Noether identities it allows to find the equivalent of the Bianchi identity and energy-momentum conservation of the standard covariant formulation. Of course these equations can also be obtained in more direct fashion via the decompositions (3.4, 3.10). A second reason is that invariance under these symmetries will explain much of the structure of the Lagrangian and the modified derivatives appearing. It will be useful to split the d + 1 dimensional diffeomorphisms into two classes, which we will refer to as time redefinitions and time-dependent spatial diffeomorphisms. This split originates in the split of an infinitesimal diffeomorphism ξ µ into Λ = ξ µ n µ and ξ i = ξ µ e i µ , via the KS frames (2.11, 2.12). Although it is customary to discuss diffeomorphisms in the passive formulation, where one keeps the coordinates fixed and the action on tensors is via the Lie derivative, we find in this case the active formulation to be more insightful. For convenience we remind the reader of the relation between infinitesimal passive and active transformations and their action on an arbitrary tensor S: It might be relevant to note how derivatives transform differently, i.e. δ P ∂ µ = 0, so that The above generalizes straightforwardly to tensors with an arbitrary number of upper and lower indices. In the remainder of this section all infinitesimal transformations will be active and we drop the A subscript.

Time redefinitions
First we consider redefinitions of the coordinate t, i.e. diffeomorphisms of the form ξ µ n µ = −MΛ, ξ µ e i µ = 0. The coordinates transform as Demanding the line element ds 2 -see (2.13) -to be form invariant, leads to the transformations Note the strong similarity to the U(1) gauge transformations of Maxwell theory, as this will provide some intuition in the discussion of invariants below. We should stress however that this is but a similarity rather than a real equality, since contrary to standard U(1) gauge transformations here also the coordinates transform, the transformations are not abelian and it is the hatted rather than standard partial derivatives that appear in (3.16-3.18).
A short computation reveals that the following objects are invariant under these time redefinitions, i.e. δ Λ X = 0, for X any of First of all, one sees that although the partial derivatives are not invariant 7 , their hatted versions are. So we see that the appearance of these hatted partial derivatives is no coincidence, but fully dictated by the invariance under time redefinitions. SinceΓ k ij is built out of h ij and∂ i , it also follows that it is invariant, which makes it the natural connection to use. Similarly it isR ij , the modified Ricci tensor which is invariant rather than the standard Ricci tensor R ij . The symmetric part 2∂ t h ij is manifestly invariant while the anti-symmetric part takes a form similar to a U(1) curvature tensor Given the list of invariants (3.19), the KS Lagrangian (2.55) is manifestly invariant. The same is true for the equations of motion by inspection of (3.5-3.7). More generally speaking, time redefinition invariance guarantees that M, C i and derivatives can only appear through the objects (3.19).

Time dependent spatial diffeomorphisms
A second set of diffeomorphisms are those for which the generating vector field satisfies ξ µ n µ = 0. The remaining components, ξ i = e i µ ξ µ , define a spatial coordinate transformation, but one that is time dependent. Explicitly It will be useful to introduce a separate notation for the time derivative of this vector field:f Invariance of ds 2 then defines the transformations of the fields: Here we introduced the following linear operation on tensor components: The introduction of C ξ is useful when performing variations of objects build out of the basic fields, since this operation is nicely compatible with tensor multiplication and contraction. E.g.
The KS action (2.55) is not manifestly invariant under the transformations (3.23-3.25), a somewhat tedious calculation reveals that (3.30) This explicitly confirms that indeed general relativity in KS form is invariant also under time-dependent spatial diffeomorphisms, albeit not manifestly so. We stress that this is nothing but a consistency check, since the KS action is identical to the Einstein-Hilbert action and the time-dependent spatial diffeomorphisms are nothing but a subclass of d + 1 dimensional diffeomorphisms, so invariance is guaranteed by construction.

Noether identities
In any theory with a gauge or local symmetry the equations of motion are not all independent and the (differential) relations between them go under the name of Noether identities (sometimes also called classical Ward identities), see e.g. [30] for a pedagogic introduction. In the case of general relativity the Noether identities amount to the Bianchi identity for the Einstein tensor and the conservation equation for the energy-momentum tensor.
Given a local symmetry of the form 8 the associated Noether identity is Here E are the Euler-Lagrange equations, defined through a generic variation of the Lagrangian as We can then apply this to the Lagrangian L KS (M, C, h) (2.55) both for the time redefinitions and time dependent spatial diffeomorphisms discussed above. The resulting Noether identities are respectively equivalent tô Of course the same equations can also be obtained by decomposing the Bianchi identity ∇ µ G µν = 0 via the methods of the previous section. Through the Einstein equations in the form (3.9) one then directly obtains the energy-momentum conservation equations:D

Conformal redefinition and the relativistic Poisson equation
The field content of GR in KS formulation, as we introduced it above, is the triple (M, C i , h ij ). But since M appears algebraically in the Lagrangian its equation of motion, eqn (3.5), is a constraint and does not involve any derivatives on M.
Second derivatives on M do appear in eqn (3.7) but this is a tensorial equation instead of a scalar equation. This suggests that M is not a relativistic analog of the Newtonian potential, nor is (3.5) the analog of the Poisson equation. Indeed in a non-relativistic expansion, see section 5 or [11,26], it is the combination G 0 +G ij h ij that becomes the Newtonian Poisson equation. This particular linear combination appears naturally as the equation of motion E ψ of a scalar ψ introduced through the field redefinition 9 Indeed, defining 10 Let us point out that from here onward in the paper we set the number of spatial dimensions d = 3.
10 Throughout the remainder of the paper we will use the notation L = √ γL and E = √ γE to distinguish between tensorial quantities and the associated densities. a comparison to (3.3) reveals that This suggests ψ is a natural relativistic analog of the Newtonian potential, and indeed this change of variables has proven very fruitful in the nonrelativistic expansion; this was the key insight of [27] that motivated the KS formulation of [26] and it also appeared in the 1/c expansion in [11,18].
Expressed in terms of the alternative fields (3.38), and upon dropping some total derivative terms, the KS action (2.55) becomes Note that if one interprets ψ as the relativistic analog of the Newtonian potential then Ω i = 1 2∂ i ψ +Ω i is a natural analog of the gravitational force vector.

Making time derivatives explicit
General relativity in the KS formulation appears -at least to us -most elegant in the form (2.55), which is simple and very similar to the ADM form. In the previous subsection we made a change of variables that is more suitable for comparison to Newtonian gravity. This introduces a few extra terms in the Lagrangian, see (3.41), but still it retains some of its elegance due to the use of the quantitiesR,K and the hatted partial derivative∂ i (see (2.44)). These quantities each contain 'hidden' time derivatives, which we would now like to make explicit. This is motivated by our aim, in sections 4 and 5, to make an expansion in (inverse) powers of c, which naturally accompany each time derivative. The key point in this subsection is that in particular the time dependence in the term √ γR can be simplified quite a lot by discarding a total derivative term. We start by introducing the speed of light c in the original relativistic metric (2.13), which together with the redefinition (3.38) is then of the form (3.44) Since apart from an overall factor in front of the action this is the only place c appears in the purely metric sector of the theory it implies all factors of c can always be re-absorbed by simply rescaling t. Equivalently this implies that each time derivative has to be accompanied by a power c −1 . Since the Lagrangian is maximally second order in derivatives it implies we can write with L 0 containing no time derivatives, L 1 being first order in time derivatives and L 2 second order. We will discard total derivatives such that L 1 and L 2 will be quadratic in derivatives rather than contain second derivatives. L 0 will still contain second spatial derivatives. One can think of L 2 as the kinetic term, of L 1 as a Lorentz coupling to velocity and of L 0 as the potential term.
To simplify notation we will introduce a dot notation for time derivatives, i.e. = ∂ t . We make the important conventional choice that the time derivative will always act on tensors with lower indices. This implies that if one sees an object with upper indices and a dot, one should first lower the indices before interpreting the dot as a time derivative! For example: Note that under our convention we also havė One could proceed by simply inserting (3.44) into the Einstein Hilbert action, but we can use the results of the previous subsections to start from (3.41) instead. First we introduce the appropriate factors of c and split all the objects appearing in (3.41) in parts containing equal powers of c: In rewritingR[γ] it is useful to isolate a total derivative (the second line below), which will drop out of the action. One calculates that where D i and R are defined with respect to the Levi-Civita connection of γ ij and we additionally introduced Inserting (3.49-3.52) and (3.53) into (3.41) leads, upon dropping total derivatives, to Let us stress again that the Lagrangian (3.45, 3.57) equals the Einstein-Hilbert Lagrangian up to a total derivative. We refer to Appendix C for the equations of motion in a form with explicit time derivatives. These can be obtained by varying (3.45, 3.57) or equivalently by rewriting (3.5-3.7). We will refrain from an extensive discussion of the symmetries in this form, as we already discussed them in some detail in a more compact form in section 3. The discussion there also implies the transformations (3.58) and (3.59) do indeed leave the action for (3.45, 3.57) invariant. Still let us shortly mention a few key points and formulas. For use in the next sections we find it convenient to revert in this and the next sections to the passive formulation. In the passive formulation, using (3.11,3.12), and upon (3.38), the time redefinition symmetry (3.16-3.18) takes the while the time-dependent spatial diffeomorphisms (3.23-3.25) become One point of importance below is that the time-redefinition symmetry (3.58) implies a simpler scaling symmetry. Consider the special case Λ = c(t − t 0 )α, with t 0 an arbitrary constant and α a real parameter. Then observe that (3.38) implies that, when evaluated at t = t 0 , the transformations of the fields and their time derivatives take the scaling forms δ α X = [X]X, δ αẊ = α([X] + 1)Ẋ, where the scaling weights are (3.60) From this it follows that [ √ γdtd 3 x] = 2, which in turn implies that invariance of the action requires [L] = −2. One can indeed check that with the weights (3.60) all the terms in (3.57) carry weight −2. Note also the simple fact that since only ∂ t and C i have an odd weight, every term with a single time derivative should carry an odd number of C i factors, something which indeed can also be verified by inspection of (3.57). This feature will play a small but key role in the shuffling algorithm introduced in section 5.1.4 and it is thus interesting to point out it has its origin in the time redefinition symmetry of the theory. The time redefinition symmetry (3.58) will also play a key role in gauge fixing the trace of the subleading metric coefficients, as discussed in section 4.1.2. In section 5.2 we will shortly discuss how the symmetries (3.58, 3.59) upon expansion lead to the non-relativistic symmetries of [8,10,11].

The 1/c expansion in KS formalism
For a pedagogic summary and overview of earlier work on the 1/c expansion see [4,31]. All of this work has been in a form making 4-dimensional coordinate invariance manifest, and for that reason the analysis was based on Newton-Cartan geometry. Here we will revisit this expansion, but now in a 3+1 formulation, i.e. one that only keeps 3-dimensional coordinate invariance manifest and uses an explicit choice of time. This is not unnatural for two related reasons. First of all in the non-relativistic expansion there is a time direction that all non-relativistic observers agree upon. Secondly, by dimensional reasons, all time derivatives are accompanied by a power of 1/c and so the 1/c expansion is actually an expansion in time derivatives. As we argued in section 1, it is the KS formalism that provides the natural 3+1 split to use for the 1/c expansion. In this 3+1 form the general structure of the expansion will be more transparent, which allows us to make some new all order observations as well as push the expansion to higher order than was previously done. In section 5 we discuss the truncation of the 1/c expansion to the 1/c 2 expansion in detail.

Structure of the expansion
The starting point is the relativistic KS Lagrangian in the form (3.45, 3.57). To perform the expansion we assume 11 the fields to be analytic in 1/c : Note that the Lagrangian (3.45, 3.57) has some explicit dependence on 1/c as well.
In combination with the 1/c dependence of the fields (4.1) one gets an expansion γ ij , . . . , γ ij ]c −n . In the following subsections we will discuss some general features of this expansion. It will often be useful to collect all fields and their expansion coefficients in a single 'master field':

Actions and equations of motion
Note that one could either expand the Lagrangian L and then vary the expansion coefficients (n) L with respect to the various expanded fields to obtain equations 11 A priori one could consider a Laurent expansion in 1/c, a transseries including terms of the form e − A c or even more generic c dependence in the fields. Such more general ansatze remain largely unexplored and fall outside the scope of this paper. of motion, or alternatively one could expand the original equations of motion associated to the unexpanded Lagrangian. In [11] the interplay and compatibility between these two approaches was explained, for completeness we shortly revisit this in appendix A.1. The key relation derived there is Φ, k = 0, . . . , n. I.e. one could expand the action up to the preferred order n and then vary it with respect to all fields up to that order, to obtain all the relevant equations of motion: In practice this approach is however a rather contrived way to find, say, E , which one could also obtain by varying (1) S with respect to (0) Φ, a much shorter calculation. Indeed from a computational point of view it is more natural to compute order by order:

Gauge fixing the trace
The non-relativistic theory obtained through the 1/c expansion contains a tower (n) γ ij of symmetric tensor fields. It is natural to identify (0) γ ij as a metric defining the geometry underlying the expanded theory. Using this metric we can split the subleading tensor fields into a trace and traceless part 12 : The trace appears frequently in the expansion, both in the equations of motion and the Lagrangian, for example through the expansion of the determinant of γ ij . As we discuss in detail in appendix A.2, one can always go to a coordinate gauge where all traces are zero, i.e.
(n) γ = 0, ∀n ≥ 1. In such a gauge the equations of motion will simplify accordingly. Concerning the action one has to be a little more careful, since gauge fixing at the level of the action can lead to a loss of equations of motion 13 . This can however be circumvented in the 1/c expansion, in the following sense. Let us define With the above we mean that the gauge fixed collection of equations of motion is equivalent to the collection of equations of motion obtained from the action which has been gauge fixed. We should point out that the precise interpretation of (4.8) is somewhat subtle when the variation with respect to (0) γ ij is considered. One can treat (n) γ ij as an independent field that does not change under variation of (0) γ ij , which is consistent if one does not forget to replace all expressions of the formγ ij F ij byγ ijF ij before varying. Alternatively one can consider γ indicates the trace of (n) γ ij for n ≥ 1, we will use (0) γ to indicate the determinant of (0) γ ij . Although this might at first appear confusing this will not lead to any clash of notation as the trace of (0) γ ij will never appear (since γ , n ≥ 1, puts extra load on the notation, it prevents the proliferation of symbols, which is already quite large. 13 These are the so called constraints, such as the well known Gauss constraint of electromagnetism that needs to be supplemented to the Euler-Lagrange equations obtained from the Lagrangian in A 0 = 0 gauge.

The leading order as the stationary sector of GR
We have chosen notation and conventions in such a way that both the Lagrangian and the fields start at order c 0 , see (3.45) and (4.1, 4.2). It follows that the leading order of the 1/c expansion is where via (3.57) ψ . (4.11) In the above R[ C i , note that furthermore all indices are raised with (0) γ ij , which is the inverse of (0) γ ij .
Since the leading order fields γ ij will appear at all lower orders as well, it will be useful to introduce notation removing their superscript. By abuse of notation we will refer to the leading order fields in the expansion by exactly the same symbol as the un-expanded fields. This should not lead to confusion as it should be clear from the context whether we are discussing the 1/c expansion or not. I.e. from here onwards: (4.12) Similarly we will simply write √ γ, R and C ij for C ij respectively. Additionally indices will be raised at all orders of the expansion with γ ij , again abusive notation for γ ]. In this new notation then, we can summarize the Lagrangian at leading order as This form immediately reveals that this leading order Lagrangian (4.13) is nothing but the full KS Lagrangian (3.45, 3.57) with all time derivatives removed. In other words, for time independent field configurations the leading order 1/c expansion is exact. If one furthermore recalls the form of the metric (3.44) then one sees that such time independent fields correspond exactly to a 4 dimensional stationary (lorentzian) metric. We can thus conclude that the leading order of the 1/c expansion captures the full non-linear dynamics of the stationary sector of general relativity. The same conclusion was reached in [18], but the approach taken in this paper has the advantage that it is much more straightforward to arrive at the leading order Lagrangian (4.13) starting from the KS Lagrangian (3.45), instead of using the fully covariant approach of [18]. Note that stationary solutions to GR, in addition to solving the dynamics described by (4.13), will also not source any higher order corrections. One can however consider quasi-stationary solutions to GR, i.e. stationary solutions in which one makes the integration constants time dependent. Then such metrics will still solve the dynamics of (4.13), since this does not contain any time derivatives, but the time dependent integration constants will lead to non-zero time derivatives which in turn will source subleading corrections. We discuss the structure of these subleading corrections in the following subsection.

The universal linear part
As we discussed in the previous subsection, the leading order equations are nonlinear equations for the leading order fields (0) Φ. But once one considers the subleading equations, they are, as we will discuss in this subsection, linear equations determining the subleading fields (n) Φ, n ≥ 1. Furthermore these equations take the schematic form (4.14) Here D 2 is a second order linear differential operator that is universal, i.e. it is the same at each order n ≥ 1. The right hand side (n) S is a source term, built out of the fields of order lower than n, and their time derivatives. This gives the equations of motion in the 1/c expansion a hierarchic structure: at each order one determines (n) Φ by solving (4.14), and this field and its time derivatives then enter the source term in the equation for the higher order fields.
To arrive at the equations in the form (4.14) we first recall that the expansion of the Lagrangian takes the form: The second observation is then that the linear part is universal. One computes from (3.57) that 14 (4.19) The equations of motion then take the form (4.14) with L lin has the same form (4.19) at each order, D 2 (n) Φ needs to be computed only once. Doing so leads to 14 We remind the reader that we have now started using notation where ψ = where Contrary to D 2 , the source term (n) S is different at each order -indeed it becomes more intricate at each successive order -and thus it needs to be computed for each order separately via (4.20). We provide the results for (1) L source and (2) L source in section 4.2.

Explicit expansion up to NNLO
We discussed the leading order Lagrangian in section 4.1.3 and the all order linear part of the Lagrangian and equations of motion in section 4.1.4. Now we will present the complete Lagrangians at the first two subleading orders, in the traceless gauge -see section 4.1.2. The additional terms that appear without this gauge choice can be found in appendix B.
For notational convenience we have decided to remove the superscripts indicating the order and instead use different characters to label the different coefficients 15 : As before we split the Lagrangians into a linear part and a source part, and additionally we organize the source part by the number of time derivatives: (1) Although the linear part is the same at all orders and was already computed above in (4.19) we will add it here as well, to provide a complete reference.

NNLO
At the next order the split of the relativistic Lagrangian (3.45, 3.57) implies (2) L 0td = After an explicit calculation one finds (2) Previously only the leading order of an expansion of GR in 1/c including odd powers was computed [18], the results in sections 4.2.1 and 4.2.2 are new. They show the computational advantage of the 3+1 formulation since the already quite involved expressions we obtained here would become quite a bit more involved to derive in their fully covariant Newton-Cartan form.

The 1/c 2 expansion in KS formalism
In this section we review how the 1/c expansion contains the 1/c 2 expansion as a self-consistent sub-theory. In addition we spell out how the 1/c 2 expansion up to order c −2n can be obtained from the 1/c expansion up to order c −n by a reshuffling of the terms. We use this to compute up to order c −4 and compare to the results in the literature: in section 5.2.1 we rederive the results of [8,11] while in section 5.2.2 we show how a further truncation of our result reproduces the 1PN order of the PN expansion.

The even power ansatz and truncation
The explicit form of the relativistic metric (3.44) reveals that c naturally appears with an odd power. This can be circumvented by defining so that the metric takes the form It is important to stress that this is much more than a simple redefinition when combined with the assumption, crucial for an expansion in 1/c, that C i is analytic in 1/c. Then (5.1) implies that B i is analytic, but more importantly also that C i is 'subleading'. I.e. the relation (5.1) should be interpreted as the non-trivial assumption that lim c→∞ C i = 0. Inserting the ansatz (5.1) into the KS action in the form (3.45, 3.57) one gets where now In summary, after the redefinition (5.1) both the metric (5.2) and Lagrangian (5.3) contain only even powers of c. This in turn implies that if we assume the fields to be analytic in 1/c 2 , i.e.
. . (5.11) This follows from the fact that the action as expressed in terms of ψ, γ ij and B i contains only even powers, see (5.3, 5.4), and the same is thus true for the equations of motion E ψ = δS δψ and E ij = δS δγ ij when expressed in terms of the fields ψ, γ ij and B i . Note that via (5.1) E i = δS δC i = c δS δB i and so this equation will only contain odd powers when expressed in terms of ψ, γ ij and B i .
The discussion above shows in detail that the truncation (5.6) from the 1/c expansion to the 1/c 2 expansion is a consistent truncation. This means it can be performed at the level of the action and that variation of the truncated action will reproduce the truncated equations of motion. Equivalently it also shows that the even coefficients do not source the odd coefficients, in case all of those are set to zero (the reverse is not true). We should point out that our discussion of the even power truncation is restricted to the pure gravitational or vacuum sector. In the presence of non-trivial energy momentum one needs to perform a similar analysis of that sector as well.

The leading order as the static sector of GR
The leading order of the 1/c 2 expansion can be obtained directly by truncating the leading order of the 1/c expansion, which we discussed in section 4.1.3. This simply amounts to removing the field (0) γ ij . From now on, as we did previously, we will simply denote This Lagrangian coincides with the fully relativistic one (3.45, 3.57) where all time derivatives and the field C i have been put to zero. Comparing to the form of the relativistic metric (3.44) we see that extrema (ψ, γ ij ) of (5.12) coincide with quasi-static solutions of GR, i.e. static solutions with time dependent integration constants. This identification of the leading order of the 1/c 2 expansion with the static sector of GR was previously made in [9].
Note now that apart from a great simplification with respect to (4.21), the linear operator (5.14) is also block diagonal, so that the field This rather straightforward observation implies that the 1/c 2 expansion will contain exactly the same terms as the 1/c expansion with the replacement (5.16).
The key difference the replacement (5.16) makes is that it changes the order of the various terms. As we'll now discuss, this change of order has a little twist, leading to a shuffling of terms among orders. But keeping track of the order and the shuffling is not too hard. We start by associating a weight 17 to the relevant objects in each expansion: Let us now consider terms in the 1/c 2 expansion and trace back their origin in the 17 The weight is the power of c −1 with which this object appears. The only factors of c appearing are those associated to the weights as listed in (5.17, 5.18). Please be aware that this notion of weight as we use it in this section is unrelated to the weight as defined in (3.60) 1/c expansion through the replacement (5.16). Every term will be polynomial 18 in the expansion coefficients. We can group Ψ = (ψ, γ) since they have the same weight. Since spatial derivatives, as well as indices, have weight zero we can ignore them. We should keep track of time derivatives, but it is not relevant on which coefficient they act, so we will simply indicate the number of time derivatives at the beginning of the expression. So in this schematic fashion a term in the 1/c 2 expansion has the form Please note that ν a and µ a , as well as λ ≤ 2, indicate positive integer powers. Via (5.18) it follows that the weight of this term is Viewed through the replacement (5.16), the term T originated from a termT in the 1/c expansion, withT is integer as it should. This follows by inspection of (5.4) or (3.57) and observing that the terms with an even number of time derivatives come with even powers of B, respectively C, while the terms with one time derivative come with an odd power of B, respectively C. This is a consequence of invariance under the time redefinition symmetry (3.58), as discussed at the end of section 3.4.
Taking this into account, together with the fact that M is positive by definition, we can conclude whereT ♯td denotes terms with ♯ time derivatives and we split the case with two time derivatives in those that contain coefficients of C and those that do not. Formula (5.23) shows that the presence of the C i and time derivatives shuffles orders under the replacement (5.16) rather than simply relating them by a factor of two. The upshot of our discussion is however (5.24) which states that all terms up to order c −2N in the 1/c 2 expansion originate from a term up to order c −N in the 1/c expansion, the only exception being those terms containing two time derivatives and no coefficients (2k) B i , which originate from terms with two time derivatives and no (k) C i at order c −N −1 in the 1/c expansion. We can transform the above conclusion to the following algorithm, which we'll refer to as shuffling, to obtain the 1/c 2 expansion up to a given order c −2N .
Shuffling algorithm Compute all terms in the 1/c expansion up to order c −N . Make the replacement (5.16). Collect the resulting terms order by order by using the rule (5.18) to determine the order. Compute the c −N −1 'th order in the 1/c expansion of the two time-derivative terms without C i factors in (3.57). Make the replacement (5.16). Add the result to the collection of terms at order c −2N .
Note that under this procedure some terms will get a weight greater than 2N, these can be discarded as they do not appear in the 1/c 2 expansion up to order c −2N . For this reason the 1/c 2 expansion up to order c −2N is simpler than the 1/c expansion up to order c −N . The shuffling algorithm can be applied to the Lagrangian as well as the equations of motion.
Let us compare the three ways to compute the 1/c 2 expansion: • Direct approach Insert the ansatz (5.5) into the relativistic Lagrangian in the form (5.3, 5.4). Expand up to order c −2N .
• Shuffling Obtain the 1/c expansion up to order c −N , use the shuffling algorithm outlined above to get the 1/c 2 expansion up to order c −2N .
• Truncation Obtain the 1/c expansion up to order c −2N , put all odd power coefficients to zero to get the 1/c 2 expansion up to order c −2N .
In the absence of any results on the 1/c expansion the direct approach is the most efficient, since the 1/c 2 expansion contains less terms than the 1/c expansion, and so it is easier to directly expand in 1/c 2 up to order c −2N than it is to first expand in 1/c up to order c −N and then to shuffle 19 . In case there is a result of the 1/c expansion up to some order c −M available, then the most efficient way to obtain the 1/c 2 expansion up to order c −2⌊ M 2 ⌋ is of course to simply truncate. But shuffling is more powerful in that case, since it directly allows to reproduce the 1/c 2 expansion all the way up to order c −2(M −1) , and with the relatively little extra work of expanding the two time-derivative, no C i , part of the Lagrangian (3.57) to order c −M −1 it provides the 1/c 2 expansion up to order c −2M .
We will use this in practice in section 5.2. Since we computed the 1/c expansion of the Lagrangian up to order c −2 in section 4 we can use the shuffling algorithm to easily find the 1/c 2 expansion of the Lagrangian up to order c −4 , which goes beyond earlier results in the literature. At the next to leading order, i.e. order c −2 , one can explicitly see that the result obtained by shuffling matches with that obtained by truncation, as well as with results in the literature obtained by the direct approach.

Explicit expansion up to NNLO
We will now perform the 1/c 2 expansion explicitly up to NNLO, i.e. order c −4 . For notational convenience we indicate the various expansion coefficients with different symbols rather than with their superscript, as we did for the 1/c expansion. Furthermore we continue to use the same (abuse of) notation introduced in section 4.1.3, indicating the leading order fields with the same symbol as the full c dependent field. More precisely our expansion ansatz is With this notation the replacement (5.16), used in the shuffling algorithm, becomes We present the results for the expansion of the Lagrangian (5.3) as L e + c −2 (2) L e + c −4 (4) L e + O(c −6 )) .

(5.27)
Note that we will present these results in the traceless gauge, see section 4.1.2, which we indicate by putting a bar on the relevant expressions. The extra terms appearing outside this gauge can be found in appendix B.

NLO
The next to leading order, i.e. order c −2 , of the Lagrangian in the 1/c 2 expansion can be most easily obtained by simply truncating (2) L as obtained in the 1/c expansion, see (4.26, 4.30). Alternatively one can also find it by performing the shuffling algorithm of section 5.1.4 to (1) L, which is given in (4.25, 4.28). Both methods give the same result, which is The 1/c 2 expansion up to this order was first discussed at the level of the equations of motion in [8] and later at the level of the action in [10,11]. The Lagrangian L NRG as given in (3.29) in [11] equals our (2) L e up to (irrelevant) total derivatives: Here (2) L e is the Lagrangian density outside the traceless gauge, which via (B.14) is related to (5.28) as (2) One can verify the equality (5.29) via the identification of our variables with those of [11]: One of the interesting insights provided in [11], is that apart from obtaining the Lagrangian by expansion of the Einstein Hilbert action, it can also be constructed purely from symmetry considerations. Let us thus shortly comment on the relation between the symmetries as discussed in [8,11] and their shape in the KS formalism. The gauge parameters Λ and ξ µ in the transformations (3.58, 3.59) are a priori themselves functions of c −1 and should thus be expanded. Compatibility with the expansion ansatz (5.25) requires 20 Ignoring δ f , we have three type of transformations of the fields. We can start with the leading order diffeomorphisms ξ i : These correspond to spatial diffeomorphisms also present in [11], but note that the extra time derivative terms in the transformation of the subleading fields correspond to an additional Milne boost in the language of [8,11]. The precise boost parameter in the conventions of [11] (HHO) is For the λ transformation (5.32) one finds while the subleading diffeomorphisms parameterized by ζ i lead to The transformations (5.37, 5.38) are related to those of e.g. [11] via the field redefinition (5.31) and the following relation between the parameters

NNLO
We now proceed to the next order, i.e. order c −4 . As far as we are aware the gravitational action has not been previously expanded to this order, keeping all even power potentials as we do. As we will discuss below, upon further truncation of most of these potentials our result reproduces the post-Newtonian expansion at 1PN order. Given our computation of (1) L and (2) L, in the 1/c expansion, see sections 4.2.1, 4.2.2, it is a surprisingly short calculation to obtain (4) L e by the shuffling algorithm of section 5.1.4. The result is: As before we presented the Lagrangian in trace-fixed form, see section 4.1.2. The additional terms present outside this gauge choice are given in appendix B.
Let us now discuss how the above Lagrangian describes an extension of the post-Newtonian expansion up to 1PN order. The metric in the 1PN approximation, see e.g. [7], is The Ricci tensor of this metric is . (5.44) and the coefficients of the various powers of c −2 provide the vacuum gravitational equations up to 1PN order. If we compare (5.41) to our relativistic metric (5.2) and expansion ansatz (5.25) one finds the identification while the field ǫ ij remains undetermined by the 1PN metric ansatz (5.41).
One sees that most of the fields that a priori can be non-trivial in the 1/c 2 expansion are assumed to be zero in the 1PN expansion. This is because, by definition, the Post-Newtonian expansion is a non-relativistic expansion around flat space, making it a weak gravity expansion as well. The 1/c 2 expansion is however a non-relativistic expansion around an arbitrary quasi-static metric and includes certain non-linear, strong gravitational effects exactly.
Although it would still be quite an effort to derive the full equations of motion (4) E from (5.40), it becomes rather easy under the assumption (5.45) where most fields are set to trivial values. From (5.40) together with (5.12) and (5.28) we can find all equations of motion up to order c −4 , the result is 21 Here we see explicitly how the equations of motion in the 1/c 2 expansion up to order c −4 reproduce (a non-degenerate linear combination of) the 1PN equations (5.42-5.44). Note that we did not include the equation of motion (4) E ij e , since it involves the field ǫ ij , which is of higher than 1PN order.

Discussion
In this paper we revisited the 1/c expansion starting from the KS 3+1 formulation of GR, a lesser known dual version of the better known ADM decomposition. Although this 3+1 formulation renders space-time diffeomorphism invariance nonmanifest, it preserves manifest spatial diffeomorphism invariance and keeps the variational principle intact. It has the advantage it makes the degrees of freedom more explicit, with the 4 dimensional Lorentzian metric being parameterized by a scalar ψ, a vector field C i and a metric γ ij , with the last two carrying spatial indices only. These fields are assumed to be analytic in the inverse speed of light 1/c and the coefficients in their series expansion form the effective fields of the 1/c expansion. All effective fields up to order c −4 are listed at the top of table 2. The 3+1 formulation makes the structure of the expansion more transparent and that allowed us to compute the effective Lagrangians to higher order than before. We extended the computation of the Lagrangians from leading order [18] to next to next to leading order in case of the 1/c expansion including odd terms, which is described by the fields listed in the second table from the left in table 2. In addition we made some all order observations as well, with a computation of the universal linear part of the expanded equations of motion and an all order gauge fixing of the trace of the metric coefficients. The 3+1 formulation also clarifies the relation between the 1/c and 1/c 2 expansions, the latter being a consistent truncation to even powers of the former. This allowed us to formulate an algorithm to compute the 1/c 2 expansion up to order c −2N from the 1/c expansion to order c −N . Using this shuffling algorithm we obtained the Lagrangians describing the 1/c 2 expansion to order c −4 . At order c −2 this matches the earlier results of [8,11], at order c −4 the result is new. A further truncation simplifies the 1/c 2 expansion to the PN expansion -in table 2 one finds all potentials included in the 1/c 2 expansion up to order c −4 listed at the second right and the corresponding potentials in the PN expansion to 1PN order on the far right -and in that case our result simplifies to the standard expressions [7].
We believe that in the form introduced in this paper, the 1/c expansion is now finally ready to be applied. Rather than a push to ever increasing order, the priority in the near future should be to explore non-trivial solutions of the equations of motion, an interpretation of the physics they describe and to study if and how this might improve on results obtained by the PN expansion. That the 1/c expansion is not an empty theory is well established. Various example solutions are discussed in [8,11,18] but these all take the form of expansions of exact solutions of GR. Although such examples are useful to gain intuition into some features of the expansion they are not teaching us anything inherently new. Of greater interest would be solutions to the 1/c expansion that provide approximations to solutions of GR in situations where no exact solution is known. Natural situations to think of would be strong gravitational systems such as Neutron stars or gravitational dynamics close to merger, where the PN expansion would be put to the limit and the extra potentials included in the 1/c expansion could play a crucial role. We hope the work in this paper will have paved some of the way for future work in 1/c expansion 1/c expansion 1/c 2 expansion PN expansion alt notation up to NNLO up to NNLO up to 1PNO Table 2: On the top all potentials/fields in the 1/c expansion up to order c −4 are listed in a precise but somewhat pedantic notation. On the left of the second line the same potentials are listed in an alternative notation. Second to left are those potentials featuring in the 1/c expansion up to NNLO, while second to right are those potentials present in the truncation to the 1/c 2 expansion up to NNLO. All the way on the right we list those potentials present in the PN expansion up to first Post-Newtonian order. On the far left the order at which these fields enter the 1/c expansion of their relativistic relativistic counterparts is indicated.
this direction. Let us shortly return to the ADM and KS formulations of general relativity. Our discussion in section 2 showed how both can be implemented in a unified way, where they appear as each others dual, in that the first is based on a preferred frame while the second is base on the particular co-frame e i = dx i , n = −M(cdt + C i dx i ). The frame becomes degenerate in the c → ∞ limit, but remains generating in the c → 0 limit, while the opposite is true for the co-frame. This implies the KS formalism is the natural 3+1 decomposition to use in the non-relativistic or galilean limit, while the ADM formalism is well suited for the ultra-relativistic or carrolian limit. This interpretation suggests a link between KS/ADM duality and Galilei/Carroll duality [32][33][34][35][36]. In the ADM, respectively KS, formalism the Lorentzian metric takes the form (2.9), respectively (2.13) where the fields (N, N i , h ij ) respectively (M, C i , h ij ) are functions of both time t and space, x i . If one however assumes these fields to be time independent then both metrics (2.9) and (2.13) are stationary metrics. These two different forms of a stationary metric are known as the Zermelo and Randers (-Papapetrou) forms. The two forms can also be argued to be dual [37] and appear in carrolian and galilean fluid dynamics, see e.g. [38,39]. A better understanding of the relation between KS/ADM duality and Galilei/Carroll as well as Randers/Zermelo duality would be interesting on its own and might be applicable to the 1/c expansion or PN expansion as well, see e.g. [40].
Finally let us mention that in the last few years the study of nonrelativistic gravity per se -i.e. independent of a relativistic counterpart -has been very active, see e.g. [41][42][43][44][45][46][47][48][49][50][51][52][53], motivated by quantum gravity, holography and string theory as well as condensed matter applications. The 1/c 2 expansion has provided some concrete examples of nonrelativistic gravity theories and provides a rather generic technique to obtain them. This has influenced advancement in the wider field of nonrelativistic gravity as well and we hope the same might be true for the results and insights provided in this paper, in particular as it is the first development of the 1/c expansion beyond leading order keeping all odd power coefficients.
121C356. UZ is supported by TÜBİTAK -2218 National Postdoctoral Research Fellowship Program with grant number 118C512. DVdB was also partially supported by the Bilim Akademisi through a BAGEP award.

A.1 Expanded action and equations of motion
The relation between the expanded action 22 and expanded equations of motion is independent of the particular theory under question, so we will keep it general. Apart from setting some definitions and notation the main result derived in this appendix that is of use in the main text is (A.3).
Let S[Φ; ǫ] be an action functional of some fields Φ, that can have some explicit dependence on an expansion parameter ǫ. In the case of interest to this paper Φ = (ψ, C i , γ ij ) and ǫ = 1/c. We'll denote the Euler-Lagrange equations associated to this action as E[Φ; ǫ] = δS δΦ . We assume the fields themselves to be analytic functions of ǫ: Φ = ∞ Note that if we expand the theory up to order N then we would have N equations of motion, but also N actions that can be varied with respect to N fields. At first there might appear a mismatch, but the key insight [11] is that there is a large degeneracy between the variations of the expanded actions. Indeed, as is shown below, one has the following relations: To see this, first observe that:

A.2 Gauge fixing the traces: details
In this appendix we discuss in detail how the traces γ ij (n) γ ij , n ≥ 1, can be put to zero by a choice of coordinate gauge and then show the equivalence (4.9). We conclude with a summary of the result and provide a detailed description of the variational principle (4.8).

A.2.2 Equivalence to gauge fixed action
Here we provide the argument leading to the equivalence (4.9). Since the (n) γ are independent of (0) ψ and (0) The same is not true for the variation with respect to To understand more precisely how these two objects are related we split (m) γ ij in its trace and trace-less part, as defined in (4.7), so that 25  17) In practice this has the effect that while we can keep all (m) γ fixed while varying (0) γ ij , this is not the case for the (m) γ kl . We can however require them to vary in as 24 To ease notation we write simply (m) γ = 0, rather than (m) γ = 0, ∀m ≥ 1. 25 In (A.16) again we use some shorthand notation, where γ ij ] should be understood as γ ij ; (1) ψ, . . . , γ ij , . . . , γ ij is kept fixed. This is equivalent to where we introduced the following projector on the traceless part: This implies the Taking into account these subtleties, the variational principle forS can then be related to the variational principle for S via the standard procedure for changing variables. One finds the following relation between the two variational principles: (A.23) 26 Remark that (A.17) only constrains the trace of the variation of (m) γ ij , enforcing it to be non-zero. A priori one could consider adding a traceless part (in kl) to the right hand side of (A.18). Excluding this, as we do, amounts to taking (m) γ ij as independent as possible of (0) γ kl .

Inversely one has
After we dealt with these technicalities, it has now become straightforward to prove (4.9). Using the relations (A.22-A.26) one computes Summarizing we find (n) This relation, together with (A.13, A.14), then immediately implies (4.9).

A.2.3 Summary
The upshot is that the Euler-Lagrange equations { E ij = 0} in traceless gauge. We now summarize the precise definition of this alternative action and the related variational principle.
The gauge fixed actions We see that the extra terms, the F ij δ (n) γ ij originating from (A.35), seem to play a key role in removing the trace of F ij in the final result. But looking at the same calculation differently, we could observe that since (n) γ ij is traceless we could have from the beginning written δ( In this way of doing the calculation the extra terms, now δ (n) γ ijF ij , actually vanish. Of the two approaches the second might seem the most efficient, namely to project from the beginning everything contracted with the (n) γ ij on its traceless part so that one can ignore the extra contribution from δ (n) γ ij . On the other hand, the first approach has the advantage of being more straightforward. In particular when time derivatives are present the first approach seems more tractable. In that case one can use, by (A.35), that This illustrates that a number of further terms appear, but partial integration will only be needed for the last one.

B Extra trace terms
In the main text we have chosen to present a number of results in a gauge where we have put to zero the traces (n) γ of the subleading tensor fields (n) γ ij , n ≥ 1, see section 4.1.2 and appendix A.2 for details on how and why this is done. For completeness we provide in this appendix the additional terms that would appear in these results if they were not gauge-fixed. This could be useful in case one would like to study certain features in a gauge independent fashion, or in another gauge.
Additions to results in section 4.1.4 The linear part of the Lagrangian gets the additional trace terms in the 1/c expansion: L can be found in (4.13). The 2nd order linear differential operators appearing in the linear part in the 1/c expansion get the additional trace terms: are leading order equations of motion, and for this reason can be ignored on-shell.

Additions to results in section 4.2
The extra contributions to the linear parts of the Lagrangian are already given at all order in (B.1). In addition to those one has that (1) L 1td = (1) L 1td and (2) L 0td = (2) L 0td L 1td = (2) L 1td ᾱ kl →α kl L 2td = L e can be found in (5.12). The 2nd order linear differential operators appearing in the linear part in the 1/c 2 expansion get the additional trace terms: γ , (B.10) γ ij γ kl  The extra contributions to the linear parts of the even Lagrangian are already given at all order in (B.9). Outside the traceless gauge, the Lagrangians (5.28, 5.40) become (2) L e = (2) L e β kl →β kl L e = (4) L e β kl →β kl ǫ kl →ǫ kl

C Equations of motion
Here we present the equations of motion of the (un-expanded) KS action. As in (3.45) the Lagrangian can be split as with the L a given in (3.57). Similarly we can split the equations of motion (in the absence of matter) as We present the different parts of the equations of motion obtained from variation by each of the fields Φ = (ψ, C i , γ ij ) below.

Variation by ψ
Variation of the KS Lagrangian with ψ yields

Variation by C i
Next, one finds the equation of motion for the gauge field C i to be + e 2ψ 2ψC j + 2Ċ j + 1 2γ C j − C kγ jk C ij + C jĊ ij + C jγk i C jk , −E i 2 = − Γ ij kl C jγkl − 1 2 (C iγjkγ jk +γ ij C jγ ) + C kγijγ jk (C.8)

Variation by γ ij
It is a bit cumbersome but straightforward to obtain the contributions: