We give a detailed canonical analysis of the ndimensional fRiemann gravity correcting the earlier results in the literature. We also write the field equations in the FischerMarsden form which is amenable to identifying the nonstationary energy on a spacelike hypersurface. We give pure RSUP2SUP and  image  theories as examples.
['constraints', 'evolution', 'generic', 'gravity', 'General Relativity and Quantum Cosmology', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'High Energy Physics - Theory', 'Mathematical Physics']
Constraints and time evolution in generic fRiemann gravity
{'Constraints and Time Evolution in Generic f (Riemann) Gravity': 'Emel Altas 1, ∗ and Bayram Tekin b2, ‡  \n1 Department of Physics, Karamanoglu Mehmetbey University, 70100, Karaman, Turkey 2 Department of Physics, Middle East Technical University, 06800, Ankara, Turkey (Dated: September 10, 2024)  \nWe give a detailed canonical analysis of the n -dimensional f (Riemann) gravity, correcting the earlier results in the literature. We also write the field equations in the Fischer-Marsden form which is amenable to identifying the non-stationary energy on a spacelike hypersurface. We give pure R 2 and R µν R µν theories as examples.', 'I. INTRODUCTION': "Within the paradigm of effective field theories, General Relativity augmented with dark matter and dark energy is the lowest-order effective theory of gravity that works remarkably well at small and large scales. Of course, it is expected to be modified at extremely high energies, or extremely short distances, for example around the black hole or Big Bang singularities. The modifications can be computed from a microscopic theory, such as string theory, alternatively, one can study generic modifications consistent with symmetries (see, for example, [1] for more arguments of the raison d'être of modified gravity theories). A large subclass of modified gravity theories is described by the action of the form 1  \nS = 1 2 ˆ M d n x √ -g f ( R µνρσ ) , (1)  \nwhere f ( R µνρσ ) is a differentiable scalar invariant of the Riemann curvature tensor (we are taking the Newton's constant κ = 1). We consider the metric tensor to be the independent field. We shall work in the metric formulation and in n dimensions. Canonical analysis of this type of theory was given in the pioneering work [3], whose notation we shall adapt here. Our analysis mostly agrees with those of [3], however, we shall make some corrections and also recast the Hamiltonian formulation of the theory in the rather beautiful Fischer-Marsden form [4]. Such a construction easily allows one to define the Killing initial (KID) and the approximate KIDs that are used in the definition of non-stationary energy contained in a spatial hypersurface [5-9]. These computations would be relevant to identifying the initial gravitational wave content. Many papers are dedicated to various aspects of the f ( R µνρσ ) theory. For example in [10], the particle spectrum, the masses of the perturbative excitations of this generic theory with one massive spin-2, one massless spin-2, and a massive spin-0  \nparticle was given around any one of its constant curvature vacua. In [11], conserved charges of the theory, such as energy-momentum and angular momentum, were constructed.  \nThe layout of this paper is as follows: In section II, we introduce the action and the field equations of f (Riemann) theories. In section III, we summarize the ADM splitting [12] of the action, which yields the Hamiltonian, the constraint equations, and then the timeevolution equations of the initial data on the hypersurface. In section IV, we introduce the construction of the nonstationary energy for f (Riemann) theories. In section V, we consider the R 2 and R µν R µν theories as examples. The computations are long, therefore, we give most of the details of the calculations in Appendices.  \nHere one of our goals is to expound upon the Arnowitt-Deser-Misner (ADM) analysis [12] and give sufficient details of the computations, so that the reader can follow all the details rather easily. We also apply our results to pure R 2 and R µν R µν theories. We made a meticulous effort to write all the proofs of our statements in the appendices, so as not to significantly cut the flow of the discussion in the main text.", 'II. FIELD EQUATIONS OF f (RIEMANN) THEORIES': 'The field equations coming from variation of the action (1) in the presence of a minimally coupled matter field are 2  \n-1 2 g µν f -R ( µ γρσ ∂f ∂R ν ) γρσ -2 ∇ σ ∇ ρ ∂f ∂R σ ( µν ) ρ = T µν , (2)  \nwhere the round brackets denote symmetrization with a factor of 1 / 2. In Appendix A, we gave the ADM decompositions of the necessary spacetime tensor fields, and then the proof of (2) is given in Appendix B. This equation includes fourth-order derivatives of the metric. One introduces auxiliary variables to simplify the ensuing discussion and lower the number of derivatives. Following [3], let us consider the "mother action":  \nS = 1 2 ˆ M d n x √ -g ( f ( ρ µνρσ ) + ϕ µνρσ ( R µνρσ -ρ µνρσ ) ) , (3)  \nwhere two auxiliary fields ( ρ, ϕ ) have been introduced. These rank-4 tensors have the same symmetries as the Riemann tensor and are assumed to be independent of each other and the metric g µν tensor. Therefore, variation of the action with respect to all fields can be written as  \n1   \nδS = 2 ˆ M d n x  δ √ -g ( f ( ρ µνρσ ) + ϕ µνρσ ( R µνρσ -ρ µνρσ ) ) (4) + √ -g ( δf ( ρ µνρσ ) + δϕ µνρσ ( R µνρσ -ρ µνρσ ) + ϕ µνρσ ( δR µνρσ -δρ µνρσ ) )   .  \nOne has δf ( ρ µνρσ ) = ∂f ∂ρ µνρσ δρ µνρσ , and the variation of the Riemann tensor reads  \nδR µνρσ = δg µλ R λ νρσ + g µλ δR λ νρσ , = δg µλ R λ νρσ + ∇ ρ ∇ [ ν δg µ ] σ + ∇ [ ρ ∇ σ ] δg µν + ∇ σ ∇ [ µ δg ν ] ρ . (5)  \nDue to the symmetries of the tensor fields, one gets  \nϕ µνρσ ( δR µνρσ -δρ µνρσ ) = δg µν R ( µ λρσ ϕ ν ) λρσ +2 ϕ σ ( µν ) ρ ∇ ρ ∇ σ δg µν -ϕ µνρσ δρ µνρσ . (6)  \nUsing integration by parts and defining the following tensor field 3  \nE µν := -R ( µ λρσ ϕ ν ) λρσ -2 ∇ σ ∇ ρ ϕ σ ( µν ) ρ -1 2 g µν ( f ( ρ λγρσ ) + ϕ λγρσ ( R λγρσ -ρ λγρσ ) ) , (7)  \none can express the variation of the total action as  \nδS = 1 2 ˆ M d n x √ -g   -E µν δg µν + δϕ µνρσ ( R µνρσ -ρ µνρσ ) + ( ∂f ∂ρ µνρσ -ϕ µνρσ ) δρ µνρσ   . (8)  \nIn Appendix C, one can see more details of this computation including the variation of the "mother action". From the last expression, one gets a set of field equations.  \n- · Variation with respect to ϕ µνρσ gives  \nR µνρσ = ρ µνρσ . (9)  \n- · Variation with respect to ρ µνρσ yields  \nϕ µνρσ = ∂f ∂ρ µνρσ = ∂f ∂R µνρσ , (10)  \nwhere for the second equality we have used the previous result.  \n- · The metric variation in the action yields  \nE µν = T µν . (11)  \nSubstituting the additional equations in the explicit form of E µν , we recover the original fourth order derivative equations (2). We assumed minimal coupling of matter and the metric, and no direct coupling of the matter to the auxiliary fields ρ and ϕ .  \nFigure 1. Slicing of the spacetime in terms of co-dimension one spatial hypersurface.  \n<!-- image -->', 'III. ADM DECOMPOSITION OF THE f (RIEMANN) THEORY': "Let us assume that the topology of the spacetime manifold is M = R × Σ, with the first factor being the time dimension and Σ being a spacelike hypersurface. In Einstein's gravity, the initial data constitute the Riemannian metric γ and the extrinsic curvature K on Σ together with the initial matter. The connection on the hypersurface satisfies the metric compatibility condition: D i γ jk = 0. See Fig. 1 for the slicing of the spacetime.  \nThe metric in terms of the lapse function and the shift vector reads  \nds 2 = ( N i N i -N 2 ) dt 2 +2 N i dt dx i + γ ij dx i dx j . (12)  \nMore explicitly, the components of the metric and the inverse metric are  \ng 00 = -N 2 + N i N i , g 0 i = N i , g ij = γ ij , (13)  \ng 00 = -1 N 2 , g 0 i = 1 N 2 N i , g ij = γ ij -1 N 2 N i N j . (14)  \nWe choose the future-pointing unit normal vector n µ as  \nn µ = ( 1 N , -N i N ) , n µ = ( -N, /vector 0) , (15)  \nwhile the extrinsic curvature in terms of the unit-normal vector reads as  \nK ij = ∇ i n j = 1 2 N ( ˙ γ ij -D i N j -D j N i ) , (16)  \nwhere ˙ γ ij := ∂ 0 γ ij . To rewrite the action (3) in terms of the ADM fields, we start with the equality  \nϕ µνρσ ( R µνρσ -ρ µνρσ ) = ϕ ijkl ( R ijkl -ρ ijkl ) + 4 ϕ ijk 0 ( R ijk 0 -ρ ijk 0 ) + 4 ϕ i 0 j 0 ( R i 0 j 0 -ρ i 0 j 0 ) . (17)  \nWe introduce the ADM splitting of the corresponding components of the Riemann tensor in Appendix A in detail. Using these results, one gets  \nϕ µνρσ ( R µνρσ -ρ µνρσ ) = ϕ ijkl ( R ijkl -ρ ijkl ) + 4 ϕ ijk 0 ( N ( D i K jk -D j K ik ) + N l R ijkl -ρ ijk 0 ) +4 ϕ i 0 j 0 ( N k N l R ikjl + N 2 K ik K k j + NN k ( D i K jk + D j K ik -2 D k K ij ) + N ( D i D j N -˙ K ij + L N K ij ) -ρ i 0 j 0 ) . (18)  \nHere L N denotes the Lie derivative along the vector field N i , which is defined as  \nL N K ij = N k D k K ij + K ki D j N k + K kj D i N k . (19)  \nFollowing [3], let us introduce the once and twice hypersurface-projected spatial tensor fields, respectively as  \nΣ ϕ ijk := γ il γ jm γ kn n µ ϕ lmnµ ≡ ϕ ijk , (20)  \nΣ ψ ij := -2 γ ik γ jl n µ n ν ϕ kµlν ≡ ψ ij , (21)  \nwhich can equivalently be written as  \nϕ ijk = -Nϕ ijk 0 , (22)  \nψ ij = -2 N 2 ϕ i 0 j 0 . (23)  \nSimilarly, for the other auxiliary field, we introduce  \nΣ ρ ijk := n µ ρ ijkµ ≡ ρ ijk , (24)  \nΣ Ω ij := n µ n ν ρ iµjν ≡ Ω ij . (25)  \nThe last two expressions explicitly read  \nρ ijk = 1 N ρ ijk 0 -N l N ρ ijkl , (26)  \nΩ ij = 1 N 2 ρ i 0 j 0 -N k N l N 2 ρ jkil -N k N ( ρ ikj + ρ jki ) , (27)  \nwhere ρ ijkl and ϕ ijkl themselves are spatial tensor fields on the hypersurface by assumption: namely, one has ρ ijkl = Σ ρ ijkl and ϕ ijkl = Σ ϕ ijkl , and for the sake of brevity, we dropped the index Σ. In other words, we can express the spacetime tensor components ρ ijk 0 and ρ i 0 j 0 in terms of purely spatial tensor fields:  \nρ ijk 0 = N ρ ijk + N l ρ ijkl , ρ i 0 j 0 = N 2 Ω ij + N k N l ρ ikjl + NN k ( ρ ikj + ρ jki ) . (28)  \nijkl  \nSo then, the contracted term (18) becomes 4  \nϕ µνρσ ( R µνρσ -ρ µνρσ ) = ϕ ijkl ( R ijkl -ρ ijkl ) (29) -4 ϕ ijk ( D i K jk -D j K ik -ρ ijk + N l N ( R ijkl -ρ ijkl ) ) -2 ψ ij ( N k N l N 2 ( R ikjl -ρ ikjl ) + 1 N ( D i D j N -˙ K ij + L N K ij ) + K ik K k j -Ω ij + N k N ( D i K jk + D j K ik -2 D k K ij -ρ ikj -ρ jki ) ) .  \nInserting all of the variations to the generalized action (3), one arrives at  \nS = ˆ M d n x N √ γ   f ( ρ, Ω) 2 + 1 2 ϕ ijkl ( R ijkl -ρ ijkl ) (30) -2 ϕ ijk ( D i K jk -D j K ik -ρ ijk + N l N ( R ijkl -ρ ijkl ) ) -ψ ij ( N k N l N 2 ( R ikjl -ρ ikjl ) + 1 N ( D i D j N -˙ K ij + L N K ij ) + K ik K k j -Ω ij + N k N ( D i K jk + D j K ik -2 D k K ij -ρ ikj -ρ jki ) )   .  \nNow let us eliminate the spatial tensors ϕ ijkl and ϕ ijk , which are nondynamical. The field equation for ϕ ijkl is  \nρ  \n=  \nR  \nijkl  \n=  \nR  \nijkl  \n+  \nK  \nik  \nK  \njl  \n-  \nK  \nil  \nK  \njk  \n,  \n(31)  \nwhile the field equation for ϕ ijk is  \nD i K jk -D j K ik -ρ ijk + N l N ( R ijkl -ρ ijkl ) = 0 . (32)  \nThe last two terms inside the round brackets cancel each other due to the constraint and one arrives at the second constraint  \nρ ijk = D i K jk -D j K ik . (33)  \nUsing the constraints, the action (30) reduces to  \nS = ˆ M d n x N √ γ ( f ( ρ, Ω) 2 -ψ ij ( K ik K k j -Ω ij + 1 N D i D j N ) + 1 N ( ψ ij ˙ K ij -ψ ij L N K ij ) ) . (34)  \nUsing integration by parts, one can rewrite the last two terms, ignoring the boundary contributions, to arrive at the action S = ' M d n x L , where the Lagrangian density is  \nΣ  \nL = N √ γ   f ( ρ, Ω) 2 -ψ ij ( KK ij + K ik K k j -Ω ij + 1 N D i D j N ) -1 N K ij ( ˙ ψ ij -L N ψ ij )   . (35)  \nTherefore, with the help of the auxiliary fields, we have managed to recast the higher derivative action as a lower derivative one as desired. It is now easier to find the Hamiltonian of the theory from this lower-derivative action.", 'A. Hamiltonian of the theory': 'The dynamical fields are ( γ ij , ψ ij ), and hence one needs to introduce the two canonical momenta corresponding to these dynamical fields. These are  \nΠ ij := δ L δ∂ 0 ψ ij , (36)  \nand  \np ij := δ L δ∂ 0 γ ij = δ L δK lm δK lm δ ˙ γ ij . (37)  \nAfter a straightforward computation, one ends up with  \nΠ ij = -√ γK ij . (38)  \nThe constraints (31, 33) can now be recast in terms of the canonical momenta as  \nρ ijk = D j ( Π ik √ γ ) -D i ( Π jk √ γ ) , (39)  \nρ ijkl = R ijkl = Σ R ijkl + 1 γ (Π ik Π jl -Π il Π jk ) . (40)  \nSimilarly one arrives at  \np ij = √ γ 4 δf δK ij + √ γ 2 ( 1 N ( L N ψ ij -˙ ψ ij ) -γ ij ψ kl K kl -Kψ ij -ψ ik K j k -ψ jk K i k ) . (41)  \nThen the Hamiltonian density  \nH = p ij ˙ γ ij + π ij ˙ ψ ij -L , (42)  \nbecomes  \nH = 2 NK ij p ij -√ γK ij L N ψ ij +2 p ij D i N j + √ γN ( -f 2 + 1 N ψ ij D i D j N -ψ ij Ω ij + Kψ kl K kl + K ij ψ jk K i k ) . (43)', 'B. Constraint equations': 'Up to a boundary contribution, one can express the Hamiltonian density as a sum of constraint equations:  \nH = NΦ 0 + N i Φ i . (44)  \nHere Φ 0 denotes the Hamiltonian constraint and Φ i denotes the momentum constraints. Using (43), one can rewrite the Hamiltonian density as  \nH = N √ γ ( 2 √ γ K ij p ij -f 2 + D i D j ψ ij -ψ ij Ω ij + KK ij ψ ij + K ij K i k ψ jk ) + N i √ γ ( -2 D k ( p k i √ γ ) -K kl D i ψ kl -2 D k ( ψ kl K li ) ) . (45)  \nEquivalently, in terms of the canonical momenta one obtains  \nH = N √ γ ( -f 2 + D i D j ψ ij -ψ ij Ω ij ) + N √ γ ( -2Π ij p ij +Π Π ij ψ ij +Π ij Π i k ψ jk ) + N i ( -2 √ γD k ( p k i √ γ ) +Π kl D i ψ kl +2 √ γD k ( ψ kl Π li √ γ )) . (46)  \nTherefore the constraints are  \nΦ 0 = √ γ ( 2 √ γ K ij p ij -f 2 + D i D j ψ ij -ψ ij Ω ij + KK ij ψ ij + K ij K i k ψ jk ) , (47)  \nΦ i = √ γ ( -2 D k ( p k i √ γ ) -K kl D i ψ kl -2 D k ( ψ kl K li ) ) . (48)  \nThe constraints vanish in a vacuum; but if there is a non-zero energy momentum tensor, then they must be equal to the corresponding projection of the energy-momentum tensor onto the initial hypersurface:  \nΦ 0 = 2 T nn = 2 N 2 ( 2 N i T 0 i -T 00 -N i N j T ij ) , (49)  \nΦ i = 2 T ni = 2 N ( N j T ij -T 0 i ) . (50)  \nIn addition, the field equations of Ω ij also are constraints: δ H δ Ω ij = 0. Finally, let us write all the constraints. 5  \n- · The Hamiltonian constraint:  \nΦ 0 = √ γ ( -f 2 + D i D j ψ ij -Ω ij ψ ij ) + 1 √ γ ( -2Π ij p ij +ΠΠ ij ψ ij +Π ij Π i k ψ jk ) = 0 . (51)  \n- · The momentum constraint:  \nΦ i = -2 √ γD k ( p k i √ γ ) +Π kl D i ψ kl +2 √ γD k ( ψ kl Π li √ γ ) = 0 . (52)  \n- · The additional constraint of the auxiliary field:  \n2 ψ ij + δf δ Ω ij = 0 .', 'C. Time evolution equations': '1. The first set: ˙ γ ij , ˙ ψ ij  \nFrom now on we are going to construct the time evolution equations. 6 The phase space variables are ( γ ij , ψ ij , p ij , Π ij ). The canonical coordinates evolve via  \n˙ γ ij = δ H δp ij , ˙ ψ ij = δ H δ Π ij . (54)  \nThe definition of extrinsic curvature leads to  \n˙ γ ij = 2 NK ij + D i N j + D j N i . (55)  \nThe relation (38) additionally yields  \n˙ γ ij = -2 N √ γ Π ij + L N γ ij , (56)  \nand so we can write  \n˙ γ ij = 2 NK ij + D i N j + D j N i = -2 N √ γ Π ij + L N γ ij . (57)  \nOn the other hand, one has  \n˙ ψ ij = N √ γ ( -2 p ij + γ ij Π kl ψ kl +Π ψ ij +Π i k ψ jk +Π j k ψ ik ) + L N ψ ij -N √ γ 2 δf δ Π ij . (58)  \n(53)  \nand  \nand also  \n2. The second set: ˙ p ij , ˙ Π ij  \nNext, we find the time-evolution equations for the canonical momenta:  \n˙ p ij = -δ H δγ ij , ˙ Π ij = -δ H δψ ij . (59)  \nThe second one is easier to obtain since we only focus on the variations with respect to ψ ij . Using the Hamiltonian (46), we have  \n˙ Π ij = -N √ γ δ δψ ij ( D k D l ψ kl -ψ kl Ω kl ) -N √ γ δ δψ ij ( ΠΠ kl ψ kl +Π mn Π n k ψ mk ) -N m δ δψ ij ( Π kl D m ψ kl +2 √ γD k ( ψ kl Π lm √ γ )) , (60)  \nwhere  \nAfter ignoring the total derivative terms, we get  \nδψ kl δψ ij = 1 2 ( δ k i δ l j + δ l i δ k j ) . (61)  \nN δ δψ ij D k D l ψ kl = D i D j N, (62)  \nN m Π kl δ δψ ij D m ψ kl = -N m D m Π ij -Π ij D m N m , (63)  \nN m δ δψ ij D k ( ψ kl Π lm √ γ ) = -1 √ γ Π m ( i D j ) N m . (64)  \nSubstituting all of these pieces, we end up with  \n˙ Π ij = √ γ ( N Ω ij -D i D j N ) -N √ γ ( ΠΠ ij +Π ik Π k j ) + √ γ L N ( Π ij √ γ ) +Π ij D k N k . (65)  \nSimilarly we can find ˙ p ij . Clearly, we can express  \n˙ p ij = N √ γ ( γ ij ( -2Π kl p kl +ΠΠ kl ψ kl +Π km Π m n ψ nk ) +Π ij Π kl ψ kl +Π i l Π j k ψ kl ) -N √ γ δ δγ ij D k D l ψ kl + N 2 √ γ δf δγ ij -N m δ δγ ij ( -2 √ γγ mn D k ( p kn √ γ ) +Π kl D m ψ kl +2 √ γD k ( ψ kl Π lm √ γ )) , (66)  \nwhere the last two terms cancel each other because of the variation with respect to the spatial metric. Taking into account the covariant derivatives correctly, we have 7  \n˙ p ij = N √ γ   γ ij ( ΠΠ kl ψ kl +Π km Π m n ψ nk -2Π kl p kl ) +Π ij Π kl ψ kl +Π i l Π j k ψ kl   (67) + √ γ 2   D k ( ψ ij D k N -2 ψ k ( i D j ) N ) + γ ij ( ND k D l ψ kl -ψ kl D k D l N )   + √ γ L N ( p ij √ γ ) + N 2 √ γ δf δγ ij + p ij D k N k .  \nIn Appendix G, we gave the construction of the constraint and time evolution equations of General Relativity using the results we have obtained in this and in the previous sections.  \nAs explained in detail in [7], the Hamiltonian form of the Einstein-Hilbert action, when extremized, leads to the Fischer-Marsden form [4] of the field equations. For the generic f ( Riemann ) theory studied here, one can also recast the Hamiltonian flow in a concise form as  \nd dt      γ ψ p Π      = J · DΦ ∗ ( γ, ψ, π, Π)( N ) , J :=      0 0 0 1 0 0 1 0 0 -1 0 0 -1 0 0 0      , (68)  \nwhere the small circle represents the usual matrix product.  \nIn this matrix equation, DΦ ∗ ( γ, ψ, π, Π) is the formal adjoint of the linearized constraint map ( D Φ( γ, ψ, π, Π)). Why the adjoint map appears in the Hamiltonian flow can be understood from the discussion in [7]. Here N is the lapse-shift vector with components ( N,N i ). Observe that there is no time-evolution when DΦ ∗ ( γ, ψ, π, Π)( N ) = 0, and these points in the space of initial data yield Killing vectors in spacetime [13, 14]. Such a description of Killing symmetries is extremely useful in understanding the amount of non-stationary energy contained in a given initial data. Here is how: if N = ξ is a Killing vector, say a stationary Killing vector, then the time evolution is trivial. The failure of N to be a Killing vector field is given as  \nDΦ ∗ ( γ, ψ, π, Π)( N ) = J -1 · d dt      γ ψ p Π      . (69)  \nNext, we discuss the non-stationary energy in this generic theory based on the approximate Killing initial data.', 'IV. NON-STATIONARY ENERGY IN f (RIEMANN) THEORIES': 'Dain [5] introduced the concept of non-stationary energy for the time-symmetric initial data in General Relativity for vacuum asymptotically flat spacetimes. That definition is  \nof which the formal adjoint is  \nbased on the notion of approximate Killing initial data (KID), which is to be defined below. Dain\'s invariant was extended to the time-asymmetric case in [6], and for asymptotically non-flat spacetimes in [7], where another definition based on the time-evolution equations was given. In [9] the construction was extended to non-vacuum spacetimes.  \nLet us briefly recap Dain\'s construction as it is not widely known and involves several subtle steps. Let the constraint covector be Φ := (Φ 0 , Φ i ) and DΦ be its linearization about a given solution initial solution. Then, DΦ ∗ is the formal adjoint of the linearized constraint map that acts on the lapse and shift vector. A crucial tool in the construction of Dain\'s invariant is Bartnik\'s operator P defined as [15]  \nP := DΦ · ( 1 0 0 -D m ) , (70)  \nP ∗ ( N ) := ( 1 0 0 D m ) · DΦ ∗ ( N ) . (71)  \nIf one uses the densitized versions of the constraints, one must also rescale the Bartnik\'s operator as  \n˜ P ∗ ( N ) := ( γ -1 / 2 0 0 1 ) · P ∗ ( N ) . (72)  \nFinally, we can write the Dain\'s invariant, I ( ξ ), that quantifies the amount of non-stationary energy in the initial data that solves the constraint equations:  \nI ( ξ ) := ˆ Σ dV P ∗ ( ξ ) · P ∗ ( ξ ) , (73)  \nwhere ξ := ( N,N i ), and P ∗ ( ξ ) := P ∗ ( N N k ) . More explicitly, in (73) one has  \n( N N i ) · ( A B i ) := NA + N i B i . (74)  \nThe important step here is the following: in the integral (73), one considers only the lapse and shift functions that satisfy a fourth-order partial differential equation that arises in the integration by parts as  \nP · P ∗ ( ξ ) = 0 . (75)  \nThis is called (by Dain) the "approximate KID equation", which admit all the Killing initial data as solutions, but has more solutions than the Killing initial data.  \nOur formulation [7, 9] of Dain\'s invariant directly involves the time evolution equations since one can write the formal-adjoint of Bartnik\'s operator as  \nP ∗ ( N ) := ( 1 0 0 D m ) · DΦ ∗ ( N ) =      1 0 0 0 0 1 0 0 0 0 D m 0 0 0 0 D m      · J -1 · d dt      γ ψ p Π      . (76)  \nThen, Dain\'s invariant for generic f (Riemann) theories in the time evolution formulation reads as  \nI ( N ) = ˆ Σ dV ( | D m ˙ γ ij | 2 + | D m ˙ ψ ij | 2 + 1 γ ( | ˙ p ij | 2 + | ˙ Π ij | 2 ) ) , (77)  \nwhere the time derivatives of the phase space fields appear. One must also be careful with the notation as one has | D m ˙ ψ ij | 2 ≡ γ ik γ jl D m ˙ ψ ij D m ˙ ψ kl .  \nLet us remark on the possible use of the results of this section. Given initial data that solves the constraints, one can identify what fraction of that data will turn into gravitational waves using the expression (77). As a fully deterministic theory, this is what one expects in gravity. Unfortunately, it is generically hard to find analytical solutions to the constraints. Therefore, one needs to compute (77) for a numerical solution. Even in the simplest case, provided by Dain [5] for asymptotically flat time-symmetric initial data in Einstein\'s theory, a numerical evaluation of the related integral that gives the non-stationary energy has not been carried out. It is an outstanding problem. 8', 'A. The R 2 theory': 'From now on, we shall adapt our results to the R 2 theory. We consider the following function  \nf ( ρ µνρσ ) = ρ 2 = g µρ g νσ g αβ g γκ ρ µνρσ ρ αγβκ (78)  \nto represent the R 2 theory, where ρ = g µρ g νσ ρ µνρσ and so  \nρ = g µρ g vσ ρ µvρσ = g vσ ( g 0 ρ ρ 0 vρσ + g iρ ρ ivρσ ) = g vσ ( g 00 ρ 0 v 0 σ + g 0 i ρ 0 viσ + g i 0 ρ iv 0 σ + g ij ρ ivjσ ) , (79)  \nwhich yields  \nρ = g 0 σ ( g i 0 ρ i 00 σ + g ij ρ i 0 jσ ) + g kσ ( g 00 ρ 0 k 0 σ + g 0 i ρ 0 kiσ + g i 0 ρ ik 0 σ + g ij ρ ikjσ ) , (80)  \nand then  \nρ = g 00 g ij ρ i 0 j 0 + g 0 k ( g i 0 ρ i 00 k + g ij ρ i 0 jk ) + g k 0 ( g 0 i ρ 0 ki 0 + g ij ρ ikj 0 ) + g kl ( g 00 ρ 0 k 0 l + g 0 i ρ 0 kil + g i 0 ρ ik 0 l + g ij ρ ikjl ) . (81)  \nThe auxiliary field ρ µvρσ has the all symmetries of the Riemann tensor. Therefore, by renaming the indices we get  \nρ = g ij g kl ρ ikjl +2 ρ i 0 j 0 ( g 00 g ij -g 0 i g 0 j ) +4 g ik g 0 j ρ 0 ijk . (82)  \nInserting the corresponding components of the inverse spacetime metric, we arrive at  \nρ = γ ij γ kl ρ ikjl -2 γ ij 1 N 2 ( N k N l ρ ikjl + ρ i 0 j 0 -N k ( ρ 0 ikj + ρ 0 jki ) ) , (83)  \nwhere we have already introduced the hypersurface projected field Ω ij = n µ n v ρ iµjv with the future pointing unit normal vector n µ = (1 /N 1 -N i /N ). Hence we get  \nρ = γ ij γ kl ρ ikjl -2 γ ij Ω ij = Σ ρ -2Ω ,  \nwhere we have used ρ ikjl = Σ ρ ikjl and Ω ij = Σ Ω ij . Then  \nf = 4 γ ij γ kl Ω ij Ω kl -4 γ ij Ω ij γ kl γ mn Σ ρ kmln + γ kl γ mn Σ ρ kmln γ ps γ ij Σ ρ pisj . (84)  \nHere we have used γ ik γ jl ρ ijkl = Σ ρ , to make the difference clear between the trace with the spacetime metric, ρ = g µρ g νσ ρ µνρσ . Then, we write  \nf ( ρ µνρσ ) = ρ 2 = ( Σ ρ -2Ω ) 2 . (85)  \nTo construct the primary constraint, ∂f/∂ Ω ij = -2 ψ ij , we need to calculate ∂f/∂ Ω ij . It is easy to prove that  \n∂f ∂ Ω ij = 8Ω γ ij -4 γ ij Σ ρ. (86)  \nTherefore the primary constraint of the auxiliary field is  \nψ ij = γ ij ( 2 Σ ρ -4Ω ) . (87)  \nRecall that in General Relativity one has ψ ij = γ ij , and now we have γ, Σ ρ, Ω dependence in the hypersurface field ψ ij . The Hamiltonian constraint (51) reduces to  \n1 2  \nΦ 0 = 2 √ γ ( D k D k Σ ρ -2 D k D k Ω+Ω 2 -4 Σ ρ 2 ) + √ γ (( Π 2 ij +Π 2 ) ( Σ ρ -2Ω ) -p kl Π kl ) , (88)  \nand the momentum constraint (52) becomes  \nΦ i = -2 √ γD k ( p k i √ γ ) +2Π D i ( Σ ρ -2Ω ) +4 √ γD k ( Π k i √ γ ( Σ ρ -2Ω ) ) . (89)', 'B. The R µν R µν theory': "In this section, we are going to evaluate the R µν R µν theory as an example. To be able to do this, first, we have to compute the space and time decomposition of the contraction ρ µν ρ µν . We have  \nf ( ρ µνρσ ) = ρ µν ρ µν = ρ 00 ρ 00 + ρ 0 i ρ 0 i + ρ i 0 ρ i 0 + ρ ij ρ ij , (90)  \nand  \nf ( ρ µνρσ ) = ρ 00 ρ 00 +2 ρ 0 i ρ 0 i + ρ ij ρ ij . (91)  \nNow we should decompose the corresponding components into the ADM variables. We start with ρ µν , which can be obtained as follows  \nρ µν = g αβ ρ αµβν = g 0 β ρ 0 µβν + g iβ ρ iµβν = g 00 ρ 0 µ 0 ν + g 0 i ρ 0 µiν + g i 0 ρ iµ 0 ν + g ij ρ iµjν . (92)  \nUsing the symmetries of ρ µνgσ, one has  \nρ µν = g 00 ρ 0 µ 0 ν + g 0 i ( ρ νiµ 0 + ρ µiν 0 ) + g ij ρ iµjν , (93)  \nwhich yields  \nρ 00 = g ij ρ i 0 j 0 . (94)  \nRecall that ρ 0000 and ρ 0 i 00 automatically vanish because of the symmetries. Similarly ρ 0i reads  \nρ 0 i = -g 0 k ρ i 0 k 0 + g kl ρ kil 0 . (95)  \nMoreover, the spatial component ρ ij can be written as  \nρ ij = g 00 ρ 0 i 0 j + g 0 k ( ρ ikj 0 + ρ jki 0 ) + g kl ρ kilj . (96)  \nRecall that, we have already introduced the hypersurface projected tensor fields ρ ijk and Ω ij via (28) and the inverse metric components. Let us reexpress ρ 00 . Using (94) we can write  \nρ 00 = ( γ mn -N m N n N 2 )( N 2 Ω mn + N k N l ρ mknl + NN k ( ρ mkn + ρ nkm ) ) , (97)  \nwhich yields  \nρ 00 = N 2 γ mn Ω mn + N k N l γ mn ρ mknl +2 NN k γ mn ρ mkn -N m N n Ω mn -N m N n N 2 N k N l ρ mknl -N m N n N N k ρ mkn , (98)  \nwhere the last two terms vanish because of the symmetries. For simplicity, we introduce  \nΩ ≡ γ mn Ω mn , Σ ρ kl ≡ γ mn ρ mknl , Σ ρ k ≡ γ mn ρ mkn ,  \nwhere ρ mknl and ρ mkn are purely spatial by assumption and therefore we removed the over Σ on these fields. Then, ρ 00 reduces to  \nρ 00 = N 2 Ω+ N k N l ( Σ ρ kl -Ω kl ) +2 NN k Σ ρ k . (99)  \nNow let us compute ρ 0i . One has (95), which yields  \nρ oi = -N m N 2 ( N 2 Ω mi + N k N l ρ mkil + NN k ( ρ mki + ρ ikm ) ) + ( γ mn -N m N n N 2 ) ( Nρ nim + N l ρ niml ) , (100)  \nand  \nThen it reduces to  \nρ 0 i = -N m Ω mi -N m N k N l N 2 ρ mkil -N m N k N ( ρ mki + ρ ikm ) + Nρ ni n + N l ρ ni n l -N m N n ρ nim N N n -N m N n N l ρ niml N 2 . (101)  \nρ 0 i = -N m Ω mi -N m N k ρ ikm N 2 + N Σ ρ i + N m Σ ρ im -N m N k ρ kim N , (102)  \nwhere ρ ikm = -ρ kim . Then, one ends up with  \nρ 0 i = N Σ ρ i + N k ( Σ ρ ik -Ω ik ) . (103)  \nSimilarly we can compute ρ ij . Using (96), one has  \nρ ij = -1 N 2 ( N 2 Ω ij + N k N l ρ ikjl + NN k ( ρ ikj + ρ jki ) ) + N m N 2 ( N ( ρ jmi + ρ imj ) + N l ( ρ jmil + ρ imjl )) + ( γ mn -N m N n N 2 ) ρ minj . (104)  \nThen, we obtain  \nρ ij = -Ω ij -N k N l N 2 ρ ikjl -N k N ( ρ ikj + ρ jki ) + N m N ( ρ jmi + ρ imj ) + N k N l N 2 ( ρ jkil + ρ ikjl ) + ρ mi m j -N k N l N 2 ρ kilj , (105)  \nwhich reduces to the following compact form  \nρ ij = Σ ρ ij -Ω ij . (106)  \nTo compute the contraction ρ µν ρ µν , we have to compute the higher indices versions of the components. Let us start with ρ 00 . One has  \nρ 00 = g µ 0 g 0 ν ρ µν = g 00 g 00 ρ 00 +2 g m 0 g 00 ρ m 0 + g n 0 g m 0 ρ nm . (107)  \nInserting the inverse spacetime metric components, we have  \nρ 00 = 1 N 4 ρ 00 -2 N m N 4 ρ 0 m + N n N m ρ nm N 4 , (108)  \nand making use of (99, 103, 106) one arrives at  \nρ 00 = 1 N 4 ( N 2 Ω+ N k N l ( Σ ρ kl -Ω kl ) +2 NN k Σ ρ k -2 NN m Σ ρ m -2 N m N k ( Σ ρ km -Ω km ) + N m N n ( Σ ρ nm -Ω mn ) ) , (109)  \nwhere all the terms, except the first one on the right-hand side of the last equation cancels each other. Therefore we arrive at a simple result  \nρ 00 = Ω N 2 . (110)  \nNow we can compute the first piece in (91), that is ρ 00 ρ 00 . One has the following  \nρ 00 ρ 00 = Ω 2 + N k N l N Ω ( Σ ρ kl -Ω kl ) +2 N k N Ω Σ ρ k . (111)  \nSimilarly, we can evaluate ρ 0 i :  \nρ 0 i = g µ 0 g νi ρ µν = g 00 g 0 i ρ 00 + g m 0 g 0 i ρ m 0 + g 00 g mi ρ 0 m + g n 0 g mi ρ nm , (112)  \nand we can write  \nρ 0 i = g 00 g 0 i ρ 00 + ρ m 0 ( g 00 g mi + g m 0 g 0 i ) + g n 0 g mi ρ nm . (113)  \nMore explicitly, one obtains  \nρ 0 i = g 00 g 0 i N 2 Ω + ( Σ ρ kl -Ω kl ) ( N k N l g 00 g 0 i + g k 0 g li + N k ( g 00 g li + g l 0 g 0 i ) ) + Σ ρ k ( 2 NN k g 00 g 0 i + N ( g 00 g ki + g k 0 g 0 i )) . (114)  \nInserting the inverse metric components, one ends up with  \nρ 0 i = -1 N ( N i Ω+ Σ ρ i ) . (115)  \nThen, the second piece in (91) becomes  \nρ 0 i ρ 0 i = -Σ ρ i Σ ρ i -1 N N i Σ ρ i Ω -N i N k N 2 Ω ( Σ ρ ik -Ω ik ) -N k N Σ ρ i ( Σ ρ ik -Ω ik ) . (116)  \n/negationslash  \nNote that ρ ij = Σ ρ ij . Let's continue with ρ ij . One can express  \nρ ij = g iµ g jν ρ µν = g i 0 g j 0 ρ 00 + ρ m 0 ( g im g j 0 + g i 0 g jm ) + g in g jm ρ nm . (117)  \nSince we will compute the contraction ρ ij ρ ij , we may use the symmetries of the indices at this step to simplify the construction from now on. Inserting the results (99, 103, 106), we obtain  \nρ ij = N 2 g i 0 g j 0 Ω+2 N Σ ρ k ( N k g i 0 g j 0 + g jk g i 0 ) + ( Σ ρ kl -Ω kl ) ( g i 0 g j 0 N k N l +2 N k g i 0 g jl + g ik g jl ) . (118)  \nAfter using the inverse metric components, the last equation reduces to  \nρ ij = N i N j N 2 Ω+ Σ ρ ij -Ω ij + 2 N N i Σ ρ j . (119)  \nThe last term in (91) is easy to construct. We can easily obtain  \nρ ij ρ ij = ( Σ ρ ij -Ω ij ) 2 + N i N j N 2 Ω ( Σ ρ ij -Ω ij ) + 2 N N i Σ ρ j ( Σ ρ ij -Ω ij ) ,  \nwhere  \n( Σ ρ ij -Ω ij ) 2 = ( Σ ρ ij -Ω ij )( Σ ρ ij -Ω ij ) . (120)  \nCollecting the pieces, ρ µν ρ µν becomes  \nf ( ρ µνρσ ) = ρ µν ρ µν = Ω 2 -2 Σ ρ i Σ ρ i + Σ ρ ij Σ ρ ij -2Ω ij Σ ρ ij +Ω ij Ω ij . (121)  \nRecall that the constraint on the auxiliary field ψ was given in (53). In our case, differentiation of f with respect to Ω ij yields  \n∂f ∂ Ω ij = ∂ ∂ Ω ij ( γ kl γ mn Ω k 1 Ω mn -2 Σ ρ i Σ ρ i + Σ ρ ij Σ ρ ij -2 Σ ρ mn Ω mn +Ω mn Ω kl γ km γ ln ) . (122)  \nWorking out the details, one has  \n∂f ∂ Ω ij = γ kl γ mn ( Ω kl ∂ Ω mn ∂ Ω ij +Ω mn ∂ Ω k 1 ∂ Ω ij ) -2 Σ ρ mn ∂ Ω mn ∂ Ω ij + γ km γ ln ( ∂ Ω mn ∂ Ω ij Ω kl + ∂ Ω kl ∂ Ω ij Ω mn ) , (123)  \nand using  \none arrives at  \n∂ Ω mn ∂ Ω ij = 1 2 ( δ i m δ j n + δ i n δ j m ) , (124)  \n∂f ∂ Ω ij = 2 ( Ω ij +Ω γ ij -Σ ρ ij ) . (125)  \nThen, the constraint equation of the auxiliary field (53) yields  \nψ ij = Σ ρ ij -Ω γ ij -Ω ij . (126)  \nWe have already introduced the Hamiltonian and the momentum constraint equations of generic f (Riemann) theories in (51, 52). Using these expressions, the momentum constraint reduces to  \nΦ i = -2 √ γD k ( P i k √ γ ) +Π kl D i ( Σ ρ kl -Ω kl ) -Π D i Ω (127) +2 √ γD k ( Π li √ γ ( Σ ρ kl -Ω kl ) ) -2 √ γD k ( ΩΠ ki √ γ ) .  \nSimilarly, the Hamiltonian constraint becomes  \nΦ 0 = √ γ 2 ( Ω 2 +2 Σ ρ i Σ ρ i -Σ ρ 2 ij +Ω 2 ij +2 D i D j ( Σ ρ ij -Ω ij ) -2 D i D i Ω ) + 1 √ γ ( -2Π ij p ij + ( Σ ρ ij -Ω ij ) ( Π ij +Π ik Π k j ) -Ω ( Π 2 +Π 2 ij ) ) . (128)", 'VI. CONCLUSIONS': "We studied the time evolution and the constraint structure of f ( Riemann )-type theories using the auxiliary fields as was done in [3] and recast the Hamiltonian flow in the compact Fischer-Marsden form [4]. This form of Einstein's equations can be considered to be a failure of the initial data to possess an exact time translation symmetry, a vantage point that led to a definition of non-stationary energy or Dain's invariant [5]. The type of theories we studied here represent a large class of theories that can be handled with two auxiliary fields, going beyond these and including the derivatives of the Riemann tensor is somewhat challenging which we shall do in a separate work. One of our motivations for this work was to give a detailed account of the computations leading to the final constraints and time evolution expressions as there are several important mistakes and omissions in the existing literature. We also gave two concrete examples: the R 2 and R µν R µν theories.", 'ACKNOWLEDGMENTS': 'E. Altas is supported by the TUBITAK Grant No. 123F353.', 'VII. APPENDICES': 'Appendix A: ADM decomposition', '1. The metric and the connection': "For the sake of completeness let us give here the ADM split of the Einstein's equations and all the relevant tensors. Using the ( n -1) + 1 dimensional splitting of the spacetime metric (12), we can express  \ng 00 = -N 2 + N i N i , g 0 i = N i , g ij = γ ij , (A1)  \nand the inverse metric components are  \ng 00 = -1 N 2 , g 0 i = 1 N 2 N i , g ij = γ ij -1 N 2 N i N j . (A2)  \nIn generic n dimensions, the spacetime metric in a matrix form reads  \ng µν = ( g 00 g 0 i g 0 i g ij ) = ( N i N i -N 2 N i N i γ ij ) . (A3)  \nTaking the determinant of the spacetime metric, we can relate the determinants of the spacetime metric and the spatial metric as  \n√ -g = N √ γ, (A4)  \nwhere we have used g = det g µν and similarly γ = det γ ij . Let Γ µ νρ denote the Christoffel symbol of the n dimensional spacetime  \nΓ µ νρ = 1 2 g µσ ( ∂ ν g ρσ + ∂ ρ g νσ -∂ σ g νρ ) , (A5)  \nand let Σ Γ k ij denote the Christoffel symbol of the n -1 dimensional hypersurface, which is compatible with the spatial metric γ ij , D k γ ij = 0, as  \nΣ Γ k ij = 1 2 γ kl ( ∂ i γ jl + ∂ j γ il -∂ l γ ij ) . (A6)  \nThen a simple computation gives one the following components of the connection  \nΓ 0 00 = 1 N ( ˙ N + N k ( ∂ k N + N i K ik ) ) , (A7)  \nand  \nand  \nand also  \nΓ i 00 = -N i N ( ˙ N + N k ( ∂ k N + N l K kl )) + N ( ∂ i N +2 N k K k i ) + ˙ N i + N k D k N i . (A10)  \nTo compute the decomposition of the field equations, we need to express the additional tensor quantities such as the Riemann and the Ricci tensor components, the scalar curvature.", '2. ADM splitting of the Riemann tensor': 'The Riemann tensor of the spacetime is defined as  \nR µ νρσ = ∂ ρ Γ µ νσ -∂ σ Γ µ νρ +Γ µ ργ Γ γ νσ -Γ µ σγ Γ γ νρ . (A11)  \nΓ 0 0 i = 1 N ( ∂ i N + N k K ik ) , Γ 0 ij = 1 N K ij , Γ k ij = Σ Γ k ij -N k N K ij , (A8)  \nΓ i 0 j = -1 N N i ( ∂ j N + K kj N k ) + NK j i + D j N i , (A9)  \nSo, it is straightforward to compute the components given below  \nR m jkl = Σ R m jkl + K jl K m k -K jk K m l + N m N ( D l K jk -D k K jl ) , (A12)  \nR 0 jkl = 1 N ( D k K jl -D j K kl ) , (A13)  \nR 0 j 0 i = 1 N ( ˙ K ij -D i D j N -N l D i K jl -2 K l ( i D j ) N l ) -K jk K k i , (A14)  \nwhere Σ R m jkl is the Riemann tensor of the hypersurface and it explicitly reads  \nΣ R m jkl = ∂ k Σ Γ m jl -∂ l Σ Γ m jk + Σ Γ m ks Σ Γ s jl -Σ Γ m ls Σ Γ s jk . (A15)  \nAlso, we need to compute R i 0 j 0 . Using the above results it becomes  \nR i 0 j 0 = N k N l Σ R ikjl + NN k ( D i K jk + D j K ik -2 D k K ij ) + N ( D i D j N -˙ K ij + L N K ij ) + N 2 K ik K k j , (A16)  \nwhere L N denotes the Lie differentiation along the shift vector N i and when it acts on the extrinsic curvature, one has  \nL N K ij = N k D k K ij + K ki D j N k + K kj D i N k . (A17)  \nMoreover, we can introduce the hypersurface projected components of the Riemann tensor as R ijkl , R ijk/vectorn and R i/vectornj/vectorn . Below we will prove the following three statements:  \nR ijkl = Σ R ijkl + K ik K jl -K il K jk , (A18)  \nR ijk/vectorn = n µ R ijkµ = D i K jk -D j K ik , (A19)  \nR i/vectornj/vectorn = n µ n ν R iµjν = 1 N ( L N K ij + D i D j N -˙ K ij ) + K ik K k j . (A20)  \nLet us start with the proof of the first statement (A18). We have  \nR ijkl = g iµ R µ jkl = g i 0 R 0 jkl + g im R m jkl = N i R 0 jkl + γ im R m jkl , (A21)  \nwhere one can express  \nR 0 jkl = ∂ k Γ 0 jl -∂ l Γ 0 jk +Γ 0 kµ Γ µ jl -Γ 0 lµ Γ µ jk = ∂ k ( 1 N K jl ) -∂ l ( 1 N K jk ) + 1 N 2 K ˙ jl ( ∂ k N + N m K mk ) + 1 N K km ( Σ Γ m jl -N m N K ˙ jl ) -1 N 2 K ˙ jk ( ∂ l N + N m K ml ) -1 N K lm ( Σ Γ m ˙ jk -N m N K jk ) , (A22)  \nwhich yields  \nR 0 jkl = 1 N ( ∂ k K jl -∂ l K jk + K km Σ Γ m jl -K lm Σ Γ m jk ) , (A23)  \nwhere  \nwhich reads  \nR ij = ∂ 0 Γ 0 ij + ∂ k Γ k ij -∂ i (Γ 0 0 j +Γ k kj ) + Γ 0 ij (Γ 0 00 +Γ k k 0 ) +Γ k ij Γ 0 0 k +Γ k kl Γ l ij -Γ 0 0 j Γ 0 0 i -Γ 0 kj Γ k 0 i -Γ 0 ki Γ k 0 j -Γ k jl Γ l ki ,  \nor in terms of the hypersurface covariant derivatives  \nR 0 jkl = 1 N ( D k K jl -D l K jk ) . (A24)  \nSimilarly we can compute R m jkl . By definition we have  \nR m jkl = ∂ k Γ m ˙ jl -∂ l Γ m ˙ jk +Γ m kµ Γ µ jl -Γ m lµ Γ µ jk , (A25)  \nwhich reduces to  \nR m jkl = Σ R m jkl + K jl K k m -K jk K m l + N m N ( D l K jk -D k K jl ) . (A26)  \nCollecting the pieces we end up with (A18).  \nLet us prove (A19). By definition, projection once can be written as  \nR ijk/vectorn = n µ R ijkµ = n 0 R ijk 0 + n l R ijkl = 1 N R ijk 0 -N l N R ijkl , (A27)  \nR ijk 0 = g 00 R 0 kji + g 0 l R l kji , (A28)  \nR ijk 0 = N l Σ R lkji + N ( D i K jk -D j K ik ) + N l ( K ki K jl -K kj K il ) . (A29)  \nUsing the last equation in (A27), one arrives at the desired result (A19).  \nWe also need to construct R i -→ nj -→ n , which reads  \nR i/vectornj/vectorn = n µ n ν R iµjν = -R 0 j 0 i + N k R 0 jki , (A30)  \nwhere the nonvanishing components of the first term are  \nR 0 j 0 i = 1 N [ ˙ K ij -D i D j N -N l D i K lj -K lj D i N l -K ik D j N k ] -K ik K k j . (A31)  \nSubstituting the results, we obtain the desired result (A20).', '3. ADM splitting of the Ricci tensor and the scalar curvature': 'Starting with the definition of the spacetime Ricci tensor  \nR ρσ = ∂ µ Γ µ ρσ -∂ ρ Γ µ µσ +Γ µ µν Γ ν ρσ -Γ µ σν Γ ν µρ , (A32)  \none can express  \nwhich yields  \nR ij = Σ R ij + KK ij -2 K ik K k j + 1 N ( ˙ K ij -N k D k K ij -D i D j N -K ki D j N k -K kj D i N k ) , (A33)  \nwhere Σ R ij denotes the ij component of the Ricci tensor on the hypersurface  \nΣ R ij = ∂ k Σ Γ k ij -∂ i Σ Γ k kj + Σ Γ k kl Σ Γ l ij -Σ Γ k jl Σ Γ l ki . (A34)  \nThe 0 i component can be written as  \nR 0 i = ∂ 0 Γ 0 0 i + ∂ k Γ k 0 i -∂ i (Γ 0 00 +Γ k k 0 ) + Γ 0 0 i Γ k k 0 +Γ k kl Γ l i 0 -Γ k 00 Γ 0 ki -Γ k 0 l Γ l ki , (A35)  \nand this expression gives us the following simple result  \nR 0 i = N j R ij + N ( D m K m i -D i K ) . (A36)  \nSimilarly, the 00 component  \nR 00 = ∂ k Γ k 00 -∂ 0 Γ k 0 k +Γ 0 00 Γ k k 0 +Γ k kl Γ l 00 -Γ k 00 Γ 0 k 0 -Γ k 0 l Γ l k 0 , (A37)  \ncan be written in a compact form as  \nR 00 = N i N j R ij -N 2 K ij K ij + N ( D k D k N -˙ K -N k D k K +2 N k D m K m k ) . (A38)  \nThen, the scalar curvature of the spacetime, R = g µν R µν , can be expressed in terms of the scalar curvature of the spatial hypersurface, Σ R = γ ij Σ R ij , as  \nR = Σ R + K 2 + K ij K ij + 2 N ( ˙ K -D k D k N -N k D k K ) . (A39)', 'Appendix B: Field equations of the f (Riemann) theory': 'Let us start with the action  \nS [ g µv ] = 1 2 ˆ M d n x √ -g f ( R µvρσ ) , (B1)  \nof which the first-order variation is  \nδS [ g µv ] = 1 2 ˆ M d n x ( δ √ -g f ( R µvρσ ) + √ -g δf ( R µvρσ ) ) , (B2)  \nwhere  \nand  \nδ √ -g = 1 2 √ -gg µν δg µν = -1 2 √ -gg µν δg µν , (B3)  \nδf ( R µvρσ ) = ∂f ∂R λγρσ δ ( g λτ R τ γρσ ) = ∂f ∂R λγρσ ( R τ γρσ δg λτ + g λτ δR τ γρσ ) . (B4)  \nHere, the linear order variation of the Riemann tensor is δR τ γρσ = ∇ ρ δ Γ τ γσ -∇ σ δ Γ τ γρ and the variation of the spacetime connection is  \nδ Γ σ µv = 1 2 g σλ ( ∇ µ δg vλ + ∇ v δg µλ -∇ λ δg µv ) . (B5)  \nSo we have  \nδf ( R µνρσ ) = ∂f ∂R λγρσ ( R τ γρσ δg λτ + g λτ ( ∇ ρ δ Γ τ γσ -∇ σ δ Γ τ γρ )) = ∂f ∂R λγρσ ( R τ γρσ δg λτ + ∇ ρ ∇ [ γ δg λ ] σ + ∇ σ ∇ [ λ δg γ ] ρ ) . (B6)  \nBy renaming the indices and using the antisymmetry of the Riemann tensor, the last equation can be written as  \nδf ( R µνρσ ) = ∂f ∂R λγρσ ( R τ γρσ δg λτ +2 ∇ ρ ∇ γ δg σλ ) . (B7)  \nCollecting the pieces, we arrive at  \nδS [ g µν ] = 1 2 ˆ M d n x √ -g ( 1 2 g µν δg µν f ( R λγρσ ) + ∂f ∂R λγρσ ( R τ γρσ δg λτ +2 ∇ ρ ∇ γ δg σλ ) ) . (B8)  \nUsing integration by parts and ignoring the boundary terms, one has  \nδS [ g µν ] = 4 ˆ M d n x √ -g  g µν δg µν f ( R λγρσ )  \n1  (B9) + δg µν ( ∂f ∂R µγρσ R ν γρσ + ∂f ∂R νγρσ R µ γρσ ) +2 δg µν ∇ σ ∇ ρ ( ∂f ∂R µσρν + ∂f ∂R νσρµ )   ,  \nand in a compact form it reads  \nδS [ g µν ] = 1 2 ˆ M d n x √ -gδg µν ( 1 2 g µν f ( R λγρσ ) + R ( µ γρσ ∂f ∂R ν ) γρσ +2 ∇ σ ∇ ρ ∂f ∂R σ ( µν ) ρ ) , (B10)  \nwhich yields the field equations  \n-1 2 g µν f ( R λγρσ ) -R ( µ γρσ ∂f ∂R ν ) γρσ -2 ∇ σ ∇ ρ ∂f ∂R σ ( µν ) ρ = T µν , (B11)  \nwhere  \nT µν = -2 √ -g δS matter δg µν . (B12)', 'Appendix C: Introducing auxiliary fields': 'To turn the field equations (B11) into a set of first-order differential equations, we start with the augmented action  \nS [ g µν , ρ µνρσ , ϕ µνρσ ] = 1 2 ˆ M d n x √ -g ( f ( ρ µνρσ ) + ϕ µνρσ ( R µνρσ -ρ µνρσ ) ) , (C1)  \nwhere the two auxiliary fields ρ µvρσ and ϕ µvρσ have all the algebraic symmetries of the Riemann tensor R µvρσ .  \nAssuming that the matter couples minimally to the metric and not to the auxiliary fields, the variation of the action reads  \nδS [ g, ρ, ϕ ] = 1 2 ˆ M d n x ( δ √ -g ( f ( ρ µνρσ ) + ϕ µνρσ ( R µνρσ -ρ µνρσ ) ) (C2) + √ -g ( δf ( ρ µνρσ ) + δϕ µνρσ ( R µνρσ -ρ µνρσ ) + ϕ µνσ ( δR µνρσ -δρ µνρσ ) ) ) ,  \nwhere  \nδR µνρσ = R λ νρσ δg µλ + 1 2 ( ∇ ρ ∇ ν δg σµ + ∇ ρ ∇ σ δg νµ -∇ ρ ∇ µ δg νσ -∇ σ ∇ ν δg ρµ -∇ σ ∇ ρ δg νµ + ∇ σ ∇ µ δg νρ ) . (C3)  \nUsing the symmetries of the fields, we have  \nϕ µνρσ ( δR µνρσ -δρ µνρσ ) = ϕ µνρσ R λ νρσ δg µλ -ϕ µνρσ δρ µνρσ (C4) + 1 2 ϕ µνρσ ( ∇ ρ ∇ ν δg σµ -∇ ρ ∇ µ δg νσ -∇ σ ∇ ν δg ρµ + ∇ σ ∇ µ δg νρ ) ,  \nwhich can be rewritten as  \nϕ µνρσ ( δR µνρσ -δρ µνρσ ) = δg µν R ( µ λρσ ϕ ν ) λρσ +2 ϕ σ ( µν ) ρ ∇ ρ ∇ σ δg µν -ϕ µνρσ δρ µνρσ . (C5)  \nInserting these results in (C2), and integrating by parts, one ends up with  \nδS [ g, ρ, ϕ ] = 1 2 ˆ M d n x √ -g   δg µν [ R ( µ λρσ ϕ ν ) λρσ +2 ∇ σ ∇ ρ ϕ σ ( µν ) ρ (C6) + 1 2 g µν ( f ( ρ λγρσ ) + ϕ λγρσ ( R λγρσ -ρ λγρσ ) ) ] + ( ∂f ∂ρ µνρσ -ϕ µνρσ ) δρ µνρσ +( R µνρσ -ρ µνρσ ) δϕ µνρσ   + I Boundary , (C7)  \nwhere the boundary terms read  \nI Boundary = ˆ M d n x √ -g ( ∇ ρ ( ϕ σ ( µν ) ρ ∇ σ δg µν ) -∇ σ ( δg µν ∇ ρ ϕ σ ( µν ) ρ )) . (C8)  \nIntroducing  \nE µν := -R ( µ λρσ ϕ ν ) λρσ -2 ∇ σ ∇ ρ ϕ σ ( µν ) ρ -1 2 g µν ( f ( ρ λγρσ ) + ϕ λγρσ ( R λγρσ -ρ λγρσ ) ) , (C9)  \nand dropping the boundary terms, one arrives at  \nδS = 1 2 ˆ d n x √ -g ( -E µν δg µν +( R µνρσ -ρ µνρσ ) δϕ µνρσ + ( ∂f ∂ρ µνρσ -ϕ µνρσ ) δρ µνρσ ) . (C10)  \nThe field equations given in Section II follow from the above variation. One can show that using the field equations for the auxiliary fields in the E µν = T µν equation, one gets back the correct second-order field equations (2), hence the consistency.', 'Appendix D: ADM splitting of the auxiliary fields': 'Let us give some details of the computations leading to the action (30). One has  \nϕ µνρσ ( R µνρσ -ρ µνρσ ) = ϕ 0 νρσ ( R 0 νρσ -ρ 0 νρσ ) + ϕ iνρσ ( R iνρσ -ρ iνρσ ) (D1) = ϕ 0 iρσ ( R 0 iρσ -ρ 0 iρσ ) + ϕ i 0 ρσ ( R i 0 ρσ -ρ i 0 ρσ ) + ϕ ijρσ ( R ijρσ -ρ ijρσ ) .  \nDue to the symmetries, it can be written as  \nϕ µνρσ ( R µνρσ -ρ µνρσ ) = 2 ϕ 0 iρσ ( R 0 iρσ -ρ 0 iρσ ) + ϕ ijρσ ( R ijρσ -ρ ijρσ ) , (D2)  \nwhich yields  \nϕ µνρσ ( R µνρσ -ρ µνρσ ) = 2 ϕ 0 i 0 σ ( R 0 i 0 σ -ρ 0 i 0 σ ) + 2 ϕ 0 ijσ ( R 0 ijσ -ρ 0 ijσ ) + ϕ ij 0 σ ( R ij 0 σ -ρ ij 0 σ ) + ϕ ijkσ ( R ijkσ -ρ ijkσ ) , (D3)  \nand then  \nϕ µνρσ ( R µνρσ -ρ µνρσ ) = 2 ϕ 0 i 0 j ( R 0 i 0 j -ρ 0 i 0 j ) + 2 ϕ 0 ij 0 ( R 0 ij 0 -ρ 0 ij 0 ) + 2 ϕ 0 ijk ( R 0 ijk -ρ 0 ijk ) + ϕ ij 0 k ( R ij 0 k -ρ ij 0 k ) + ϕ ijk 0 ( R ijk 0 -ρ ijk 0 ) + ϕ ijkl ( R ijkl -ρ ijkl ) . (D4)  \nUsing the symmetries, one ends up with  \nϕ µνρσ ( R µνρσ -ρ µνρσ ) = ϕ ijkl ( R ijkl -ρ ijkl ) + 4 ϕ ijk 0 ( R ijk 0 -ρ ijk 0 ) + 4 ϕ i 0 j 0 ( R i 0 j 0 -ρ i 0 j 0 ) .  \nMoreover, using the decomposition of the components of the Riemann tensor one obtains  \nµνρσ  \nϕ  \n(  \nR  \nDefining  \n-  \nρ  \nϕ  \nijkl  \n+4  \n+4  \n(  \nϕ  \nϕ  \nijkl  \nR  \nijk  \n0  \ni  \n0  \nj  \n0  \n(  \n-  \nρ  \nN  \n(  \nijkl  \nD  \n(  \nNN  \ni  \n)  \n(D5)  \nK  \njk  \nD  \n(  \ni  \n-  \nD  \nK  \njk  \nj  \n+  \nK  \nik  \nD  \nj  \n) +  \nK  \nik  \nN  \n-  \nl  \nR  \nijkl  \nD  \n2  \nk  \n-  \nρ  \nK  \nij  \n)  \n+ N ( -˙ K ij + L N K ij + D i D j N ) + N 2 K ik K k j + N k N l R ikjl -ρ i 0 j 0 ) .  \nϕ ijk ≡ γ il γ jm γ kn n µ ϕ lmnµ , ψ ij ≡ -2 γ ik γ jl n µ n v ϕ kµlν , ρ ijk ≡ n µ ρ ijkµ, Ω ij ≡ n µ n ν ρ iµjν , (D6)  \n) =  \nµνρσ  \nµνρσ  \nk  \nijk  \n0  \n)  \none can express the corresponding components of the auxiliary field ϕ µνρσ in terms of the spatial tensor fields as follows  \nϕ ijk 0 = -ϕ ijk N , ϕ i 0 j 0 = -ψ ij 2 N 2 , (D7)  \nand similarly, in terms of the spatial tensors, we can write the components of ρ µνρσ as  \nρ ijk 0 = Nρ ijk + N l ρ ijkl , (D8) ρ ijk 0 = N 2 Ω ij + N k N l ρ ikjl + NN k ( ρ ikj + ρ jki ) .  \nThen, we arrive at  \nϕ µνρσ ( R µνρσ -ρ µνρσ ) = ϕ ijkl ( R ijkl -ρ ijkl ) -4 ϕ ijk ( D i K jk -D j K ik -ρ ijk + N l N ( R ijkl -ρ ijkl ) ) -2 ψ ij ( N k N l N 2 ( R ikjl -ρ ikjl ) + K ik K k j -Ω ij + N k N ( D i K jk + D j K ik -2 D k K ij -ρ ikj -ρ jki ) + 1 N ( -˙ K ij + L N K ij + D i D j N ) ) , (D9)  \nfrom which one obtains (30).', 'Appendix E: The constraints': 'By definition, the canonical momenta are defined as  \nΠ ij := δ L δ∂ 0 ψ ij , p ij := δ L δ∂ 0 γ ij = δ L δK lm δK lm δ ˙ γ ij , (E1)  \nwhere the Lagrangian density is given in (35), and so we can express  \nΠ ij = -√ γK ij . (E2)  \nIn the computation of p ij , one needs to compute the variation of the extrinsic curvature with respect to ˙ γ ij . By definition, we have  \nδK lm δ ˙ γ ij = 1 2 N δ δ ˙ γ ij ( ˙ γ lm -D l N m -D m N l ) = 1 2 N δ i l δ j m . (E3)  \nOne also uses the explicit form of the Lagrangian (35) to arrive at  \nδ L δK lm = N √ γ 2 δf δK lm + N √ γ ( 1 N ( L N ψ lm -˙ ψ lm ) -γ lm ψ ps K ps -Kψ lm -ψ ln K m n -ψ mn K l n ) .  \nCollecting the pieces we express the conjugate momenta as  \np ij = √ γ 4 δf δK ij + √ γ 2 ( 1 N ( L N ψ ij -˙ ψ ij ) -γ ij ψ kl K kl -Kψ ij -ψ ik K k j -ψ jk K i k ) .  \nThe last expression directly yields the velocities as  \n˙ ψ ij = N 2 δf δK ij -2 N √ γ p ij + L N ψ ij + N ( -γ ij ψ kl K kl -Kψ ij -ψ ik K j k -ψ jk K i k ) . (E4)  \nThe Hamiltonian density reads  \nH = p ij ˙ γ ij +Π ij ˙ ψ ij -L , (E5)  \nand inserting the velocities together with the Lagrangian density (35) it can be expressed as follows  \nH =2 NK ij p ij -√ γK ij L N ψ ij +2 p ij D i N j (E6) + √ γN ( -f 2 + 1 N ψ ij D i D j N -ψ ij Ω ij + Kψ kl K kl + K ij ψ ik K j k ) .  \nUp to a boundary term, the Hamiltonian density can be written as a sum of the constraint equations. Namely, one has  \nH = N Φ 0 + N i Φ i . (E7)  \nAfter a straightforward computation, one obtains  \nH = N √ γ ( 2 √ γ K ij p ij -1 2 f ( ρ µvρσ ) + D i D j ψ ij -ψ ij Ω ij + Kψ kl K kl + K ij K j k ψ ik ) + N i √ γ ( -2 D k ( p k i √ γ ) -K kl D i ψ kl -2 D k ( ψ kl K li ) ) . (E8)  \nSince we have already obtained the relation between the extrinsic curvature and the conjugate momenta, K ij = -Π ij / √ γ , we can equivalently write  \nH = N √ γ ( D i D j ψ ij -ψ ij Ω ij -1 2 f ( ρ µvρσ ) ) + N √ γ ( -2Π ij p ij +ΠΠ ij ψ ij +Π ij Π j k ψ ik ) + N i ( -2 √ γD k ( p k i √ γ ) +Π kl D i ψ kl +2 √ γD k ( ψ kl Π li √ γ )) , (E9)  \nwhich yields the Hamiltonian and the momentum constraints as  \nΦ 0 = √ γ ( D i D j ψ ij -ψ ij Ω ij -1 2 f ( ρ µvρσ ) ) + 1 √ γ ( -2Π ij p ij +Π ij ψ ij +Π ij Π j k ψ ik ) , (E10)  \nΦ i = -2 √ γD k ( p k i √ γ ) +Π kl D i ψ kl +2 √ γD k ( ψ kl Π li √ γ ) . (E11)', 'Appendix F: Time evolution equations': 'As for the dynamical equations, the first set of the evolution equations is  \n˙ γ ij = δ H δp ij , ˙ ψ ij = δ H δ Π ij . (F1)  \nThis set gives the velocities in terms of canonical variables. Ignoring the total derivative terms, the Hamiltonian density can be expressed as  \nH = N √ γ ( D k D l ψ kl -ψ kl Ω kl -1 2 f ( ρ µvρσ ) ) + N √ γ ( -2Π kl p kl +Π Π kl ψ kl +Π kl Π k m ψ ml ) + N m ( -2 √ γD k ( p k m √ γ ) +Π kl D m ψ kl +2 √ γD k ( ψ kl Π lm √ γ )) . (F2)  \nThen automatically we obtain  \n˙ γ ij = -2 N √ γ Π ij + L N γ ij = 2 NK ij + D i N j + D j N i , (F3)  \nwhich is the expected result (57) and by definition Lie derivative yields L N γ ij = D i N j + D j N i . Similarly ˙ ψ ij can be written as  \n˙ ψ ij = N √ γ ( -2 p ij + γ ij Π kl ψ kl +Π ψ ij +Π i k ψ kj +Π j k ψ ki ) + L N ψ ij -N √ γ 2 δf δ Π ij , (F4)  \nwhere one can evaluate the last term easily for a given f .  \nNow, we can continue with the second set of the evolution equations. One has the following relations  \n˙ p ij = -δ H δγ ij , ˙ Π ij = -δ H δψ ij . (F5)  \nUsing the Hamiltonian density again let us construct ˙ Π ij . We have  \n˙ Π ij = N √ γ Ω ij -N √ γ Π Π ij -N √ γ Π ik Π k j -N √ γ δ δψ ij ( D k D 1 ψ kl ) -N m Π kl δ δψ ij ( D m ψ kl ) -2 N m δ δψ ij ( D k ( ψ kl Π lm ) ) , (F6)  \nand we have to compute the last three terms. We have  \nD m ψ kl = ∂ m ψ kl + Σ Γ k mn ψ nl + Σ Γ l mn ψ kn , (F7)  \nwhich yields the following variation  \nδD m ψ kl = ∂ m δψ kl + Σ Γ k mn δψ n 1 + Σ Γ l mn δψ kn + ψ nl δ Σ Γ k mn + ψ kn δ Σ Γ l mn , (F8)  \nand it can be written in a more compact form as  \nδD m ψ kl = D m δψ kl + ψ nl δ Σ Γ k mn + ψ kn δ Σ Γ l mn . (F9)  \nThe variation of the hypersurface connection can be expressed as  \nδ Σ Γ k mn = 1 2 γ kp ( D m δγ np + D n δγ mp -D p δγ mn ) . (F10)  \nTherefore in a more explicit form, one obtains  \nδD m ψ kl = D m δψ kl + 1 2 ψ nl γ kp ( D m δγ np + D n δγ mp -D p δγ mn ) + 1 2 ψ kn γ 1 p ( D m δγ np + D n δγ mp -D p δγ mn ) , (F11)  \nwhich can be reexpressed as  \nδD m ψ kl = D m δψ kl + ψ n ( k γ l ) p ( D m δγ np + D n δγ mp -D p δγ mn ) , (F12)  \nand directly yields variation with respect to ψ ij as  \nN m Π kl δ ( D m ψ kl ) δψ ij = -N m D m Π ij -Π ij D m N m . (68)  \nNote that there is no contribution coming from the variations of the connection in the last expression. Similarly we can compute δD k D l ψ kl as  \nD k D l ψ kl = ∂ k ( ∂ l ψ kl + Σ Γ k lm ψ ml + Σ Γ l lm ψ km ) + Σ Γ k km ( ∂ l ψ ml + Σ Γ m ln ψ nl + Σ Γ l ln ψ mn ) . (F13)  \nTaking the variation, we write  \nδD k D l ψ kl = D k D l δψ kl + D p ( ψ ml D l δγ mp ) -1 2 D k ( ψ ml D k δγ lm ) + 1 2 D k ( ψ km γ lp D m δγ lp ) + 1 2 γ kp D l ψ ml D m δγ kp . (F14)  \nSince we focus on variation of the spatial field ψ , we only consider the first term on the right-hand side. Hence, ignoring the total derivative terms, we get  \nN √ γ δ δψ ij ( D k D l ψ kl ) = √ γD i D j N. (F15)  \nNow we should compute the last term: δD k ( ψ kl Π 1 m ) . The variation of this term gives us  \nδD k ( ψ kl Π lm ) = D k δ ( ψ kl Π lm ) -ψ kl Π ln δ Σ Γ n km . (71)  \nUp to a boundary term we get  \nN m δ δψ ij D k ( ψ kl Π lm ) = -Π m ( i D j ) N m . (72)  \nInserting them in ˙ Π ij , we end up with the desiring evolution equation  \n˙ Π ij = √ γ ( N Ω ij -D i D j N ) -N √ γ ( ΠΠ ij +Π ik Π k j ) + L N Π ij +Π ij D k N k . (F16)  \nThen, we can construct ˙ p ij = -δ H /δγ ij . Since we have the variations  \nδ √ γ δγ ij = 1 2 √ γγ ij , δγ -1 / 2 δγ ij = -1 2 γ -1 / 2 γ ij , (F17)  \nfrom the Hamiltonian density (46) we directly obtain  \n˙ p ij = -N 2 √ γγ ij ( D k D l ψ kl -ψ kl Ω kl -1 2 f ( ρ µvρσ ) ) -N √ γ δ δγ ij ( D k D l ψ kl -1 2 f ( ρ µvρσ ) ) + N 2 √ γ γ ij ( -2Π kl p kl +Π Π kl ψ kl +Π kl Π k m ψ ml ) -N √ γ δ δγ ij ( Π mn γ mn Π kl ψ kl +Π kl Π mn γ kn ψ ml ) -N m δ δγ ij ( -2 γ mn √ γD k ( p kn √ γ ) +Π kl D m ψ kl +2 √ γD k ( ψ kl Π lm √ γ )) . (F18)  \nTo simplify the last equation, we use the Hamiltonian constraint (51) together with  \nδγ mn δγ ij = -1 2 ( γ mi γ jn + γ mj γ in ) , (F19)  \nin (F18) to arrive at  \n˙ p ij = N √ γ   γ ij ( Π Π kl ψ kl +Π kl Π k m ψ ml -2Π kl p kl ) +Π ij Π kl ψ kl +Π i l Π j k ψ kl   -N √ γ δ δγ ij ( D k D l ψ kl ) + N 2 √ γ δ δγ ij f ( ρ µvρσ ) -N m δ δγ ij ( -2 γ mn √ γD k ( p kn √ γ ) +Π kl D m ψ kl +2 √ γD k ( ψ kl Π lm √ γ )) . (F20)  \nA straightforward calculation gives us  \nand  \nand also  \nN √ γ δ δγ ij ( D k D l ψ kl ) = -√ γ 2 ( D k ( ψ ij D k N -2 ψ k ( i D j ) N ) + γ ij ( ND k D l ψ kl -ψ kl D k D l N ) ) , (F21)  \n2 N m δ δγ ij ( γ mn √ γD k ( p kn √ γ )) = √ γ L N ( p ij √ γ ) + p ij D k N k , (F22)  \nN m δ δγ ij ( Π kl D m ψ kl +2 √ γD k ( ψ kl Π lm √ γ )) = 0 . (F23)  \nFinally, after collecting the pieces, one ends up with the last evolution equation (67)  \n˙ p ij = N √ γ   γ ij ( Π Π kl ψ kl +Π kl Π k m ψ ml -2Π kl p kl ) +Π ij Π kl ψ kl +Π i l Π j k ψ kl   + √ γ 2 ( D k ( ψ ij D k N -2 ψ k ( i D j ) N ) + γ ij ( ND k D l ψ kl -ψ kl D k D l N ) ) + N 2 √ γ δf δγ ij + √ γ L N ( p ij √ γ ) + p ij D k N k . (F24)  \nAt this point, one needs to know the explicit form of the function f to proceed further.', 'Appendix G: Application to General Relativity': "In this section, we will use the generic results to evaluate Einstein's theory. Let us set the function f to  \nf ( ρ µνρσ ) = g µρ g νσ ρ µνρσ . (G1)  \nSince ρ µνρσ has the symmetries of the Riemann tensor by assumption, one can easily show that  \nf = g ik g jl ρ ijkl +4 g 0 j g ik ρ 0 ijk +2 ρ 0 i 0 j ( g 00 g ij -g 0 i g 0 j ) . (G2)  \nWe can insert the components of the inverse spacetime metric to arrive at  \nf ( ρ µνρσ ) = γ ik γ jl ρ ijkl + 2 N 2 γ ij ( 2 N k ρ 0 ikj -ρ 0 i 0 j -N k N l ρ ikjl ) . (G3)  \nUsing the hypersurface projected tensor fields that we introduced before, we arrive at  \nf = γ ik γ jl ρ ijkl -2 γ ij Ω ij . (G4)  \nNow we can construct the constraint on the auxiliary field ψ ij , (53). It directly yields  \nψ ij = γ ij . (G5)  \nTherefore, the first set of evolution equations (57, 58) are related. Since ψ ij = γ ij , we can rewrite (58) as  \n˙ γ ij = -D i N j -D j N i + 2 N √ γ ( Π γ ij +Π ij -p ij -Π kl ∂f ∂ρ ikjl ) , (G6)  \nwhere we can use  \n∂f ∂ρ ikjl = γ pq γ nm ∂ρ pnqm ∂ρ ikjl . (G7)  \nWe have to preserve the symmetries on both sides of the equation. So, we have to express it in a more correct form as  \n∂f ∂ρ ikjl = 1 4 γ pq γ nm ∂ ∂ρ ikjl ( ρ pnqm + ρ qmpn -ρ npqm -ρ pnmq ) , (G8)  \nwhich yields  \nThen we have  \n∂f ∂ρ ikjl = 1 2 ( γ ij γ kl -γ il γ kj ) . (G9)  \n˙ γ ij = -D i N j -D j N i + 2 N √ γ ( 1 2 Π γ ij + 3 2 Π ij -p ij ) , (G10)  \nUsing the basic relation ˙ γ ij = -γ ik γ jl ˙ γ kl , one can directly rewrite the last expression as  \n˙ γ ij = D i N j + D j N i + 2 N √ γ ( -1 2 Π γ ij + p ij -3 2 Π ij ) . (G11)  \nFor consistency with the time evolution of the spatial metric  \n˙ γ ij = -2 N √ γ Π ij + D i N j + D j N i , (G12)  \none needs  \nΠ ij = 2 ( p ij -p n γ ij ) . (G13)", '1. Constraint equations': "One can reexpress the Hamiltonian density of General Relativity, using Ψ ij = γ ij and (G13) as  \nH = ( -p ij + 2 p n γ ij ) ˙ γ ij -L . (G14)  \nWe then introduce the new momentum  \nπ ij := -p ij + 2 p n γ ij , (G15)  \nwhich yields the trace  \nπ = n -2 n p. (G16)  \nThen, the Hamiltonian density becomes  \nH = π ij ˙ γ ij -L . (G17)  \nOne has the inverse relations  \nΠ ij = -2 π ij + 2 n -2 γ ij π, (G18)  \nand Π = 2 π/ ( n -2). Additionally, we can express  \np ij = -π ij + 2 n -2 γ ij π, p = n n -2 π. (G19)  \nFinally, the Hamiltonian and the momentum constraint equations (51, 52) reduce to  \nΦ 0 ( γ, π ) = 2 √ γ ( π 2 ij -π 2 n -2 ) -√ γ 2 Σ R, (G20)  \nand  \nΦ i ( γ, π ) = -2 D k π k i . (G21)  \n- [2] B. Tekin, A Tribute to S. Deser: Conserved Quantities in Generic Gravity Theories, [arXiv:2307.12758 [gr-qc]].\n- [3] N. Deruelle, M. Sasaki, Y. Sendouda, D. Yamauchi, Hamiltonian formulation of f(Riemann) theories of gravity, Prog. Theor. Phys. 123 :169-185, 2010.\n- [4] A. E. Fischer and J. E. Marsden, Linearization stability of the Einstein equations, Bull. Amer. Math. Soc., 79 , 997-1003 (1973).\n- [5] S. Dain, A New Geometric Invariant on Initial Data for the Einstein Equations, Phys. Rev. Lett. 93 , 23, 231101 (2004).\n- [6] J. A. V. Kroon and J. L. Williams, Dain's invariant on non-time symmetric initial data sets, Class. Quantum Grav. 34 , 12, 125013, (2017).\n- [7] E. Altas and B. Tekin, Nonstationary energy in general relativity, Phys. Rev. D 101 , no.2, 024035 (2020).\n- [8] R. Sansom and J. A. V. Kroon, Dain's invariant for black hole initial data, Class. Quant. Grav. 40 , no.11, 115002 (2023).\n- [9] E. Altas and B. Tekin, Nonstationary energy of perfect fluid sources in general relativity, Phys. Rev. D 109 , no.4, 4 (2024).\n- [10] B. Tekin, Particle Content of Quadratic and f ( R µνσρ ) Theories in ( A ) dS , Phys. Rev. D 93 , no.10, 101502 (2016).\n- [11] C. Senturk, T. C. Sisman and B. Tekin, Energy and Angular Momentum in Generic F(Riemann) Theories, Phys. Rev. D 86 , 124030 (2012).\n- [12] R. Arnowitt, S. Deser and C. W. Misner, Dynamical Structure and Definition of Energy in General Relativity, Phys. Rev. , 116 , 1322-1330, (1959).\n- [13] V. Moncrief, Spacetime symmetries, and linearization stability of the Einstein equations. I, J. Math. Phys. 16 , 493-498 (1975).\n- [14] R. Beig & P. T. Chruściel, Killing initial data, Class. Quantum Grav. 14 , A83 (1997).\n- [15] R. Bartnik, Phase space for the Einstein equations, Communications in Analysis and Geometry 13 , 845 (2005).\n- [16] A. Ashtekar, R. Maartens, Remembering Sergio Dain, Gen. Relativ. Gravit. 48, 81 (2016)."}