xPPN: An implementation of the parametrized post-Newtonian formalism using xAct for Mathematica

We present a package for the computer algebra system Mathematica, which implements the parametrized post-Newtonian (PPN) formalism. This package, named xPPN, is built upon the widely used tensor algebra package suite xAct, and in particular the package xTensor therein. The main feature of xPPN is to provide functions to perform a proper $3+1$ decomposition of tensors, as well as a perturbative expansion in so-called velocity orders, which are central tasks in the PPN formalism. Further, xPPN implements various rules for quantities appearing in the PPN formalism, which aid in perturbatively solving the field equations of any metric theory of gravity. Besides Riemannian geometry, also teleparallel and symmetric teleparallel geometry are implemented.


I. INTRODUCTION
The parametrized post-Newtonian (PPN) formalism [1][2][3][4][5][6][7] is an indispensable tool for testing the viability of gravity theories. It is used to characterize any given theory of gravity by a set of ten parameters, which are closely related to the phenomenological properties of the theory under consideration. This allows to compare the parameters which are obtained theoretical through the PPN formalism with their values measured in experiments within the solar system and related physical settings.
In order for the PPN formalism to be applicable, the gravity theory under consideration must satisfy a number of assumptions. The most important is the existence of a metric governing the motion of test bodies, which can be described by a perturbation around a flat Minkowski background. Another assumption concerns the source of gravity, which is chosen to be a fluid satisfying the covariant conservation of energy-momentum, encoded in the Euler equations of fluid dynamics. Further, it is assumed that the field equations are of second derivative order, or can be brought into this form, so that their solution take a known standard form in terms of particular Poisson-like integrals over the source matter distribution.
In order to determine the post-Newtonian limit of a given gravity theory, one must perform a 3 + 1 split of its field equations (which are, in general, tensor equations) into space and time components, and then perform a perturbative expansion around a fixed vacuum solution. Depending on the structure of the field equations, both tasks may be challenging to perform by hand, and so the use of computer algebra comes as a useful tool. Due to the common assumptions on which the PPN formalism is based, and the similar steps to be applied to different theories of gravity, it appears natural to implement these common tasks into a general computational tool, which can then be used to calculate the post-Newtonian limit for any given theory. A number of functions to achieve this has already been implemented in Maple [8]. The aim of xPPN [9] is to provide another implementation, based on a more extensive framework for performing tensor calculations.
A very powerful tensor algebra software is the xTensor package, which is part of the xAct suite of Mathematica packages [10]. It comes with numerous functions to define and manipulate tensorial expressions, and includes concepts such as induced metrics on hypersurfaces orthogonal to a vector field, or component calculations in xCoba, which can be used to achieve a 3+1 split of tensorial expressions. The former is employed, for example, for cosmological perturbation theory in the xPand package [11], in conjunction with the xPert package [12] for general tensor perturbations. It thus appears natural to follow a similar approach to implement also the PPN formalism in xAct. However, the latter comes with a few peculiarities which obstruct this strategy. The main difficulty lies in the different treatment of space and time in the PPN formalism, such as assigning different perturbation orders to space and time components of tensors. Another, albeit smaller issue is the different convention for counting the perturbation orders in xPert and the PPN formalism.
Even though it is possible to implement the PPN formalism also using the existing functionality mentioned above, it appears simpler to overcome the aforementioned difficulties by using a different approach both to the 3 + 1 split and the perturbative expansion, without using the induced metric framework or the xCoba and xPert packages. This is the approach followed by xPPN. The key idea is to complement every tensor defined on the spacetime manifold by a number of tensors on a purely spatial manifold, which additionally depend on an external time parameter, and which represent the separated time and space components of the original spacetime tensor. Also derivatives are split into spatial derivatives and derivatives with respect to the time parameter. This approach allows an essentially different treatment of perturbations for space and time components. The aim of this article is to present this approach, as well as the implementation of the PPN formalism which is based on this framework. This article is structured as follows. We start with a brief review of the PPN formalism and its extensions implemented by xPPN in section II. The general concepts on which this implementation is based are explained in section III. The main functionality of xPPN is laid out in section IV, where we display the geometric objects defined by xPPN, and section V, where we explain the functions provided for the common operations on tensor expressions. A complete usage example is given in section VI, which shows how to calculate the PPN parameters of a simple scalar-tensor theory of gravity. A summary and outlook towards possible future extensions are given in section VII. For further reference, we provide notes on the implementation and internal operation of xPPN in appendix A.
In order to make code listings more readable, the following syntax highlighting is used. Mathematica keywords are typeset in blue, xAct keywords are typeset in green and xPPN keywords are typeset in red. We use lowercase letters for coordinate indices and uppercase letters for Lorentz indices, where in both cases Greek indices run from 0 to 3 and belong to spacetime, while Latin indices run from 1 to 3 and belong to space only.

II. PARAMETRIZED POST-NEWTONIAN FORMALISM
The aim of the xPPN package is to provide an implementation of the parametrized post-Newtonian (PPN) formalism and several of its geometric extensions. Here we briefly summarize these theoretical foundations. We first discuss the standard PPN formalism and the perturbative expansion of the metric in section II A. An extended formulation using a tetrad as the fundamental field instead of the metric is shown in section II B. We then display its application to teleparallel gravity in section II C and symmetric teleparallel gravity in section II D.

A. Standard formalism for Riemannian geometry
There exist different versions of the PPN formalism. Here we adhere to its form presented in [5]. Its starting point is the assumption that the propagation of light and massive test bodies is governed by a pseudo-Riemannian metric g αβ of Lorentzian signature. This metric is further considered in a weak field limit as an asympotically flat perturbation around a flat Minkowski background. The source of this perturbation is assumed to be of compact support and modeled by the energy-momentum tensor Θ αβ = (ρ + ρΠ + p)u α u β + pg αβ (1) of a perfect fluid with four-velocity u α , rest energy density ρ, specific internal energy Π and pressure p. Further, a fixed Cartesian coordinate system (x α ) = (t, x a ) is used, denoted as the universe rest frame. It is then assumed that the velocity v a = u a /u 0 of the source matter in this coordinate system is small compared to the speed of light, |v a | ≪ c ≡ 1. This velocity takes the role of the perturbation parameter. Physical quantities, such as the metric, are expanded in terms of the source velocity, and any term in this perturbative expansion which is proportional to |v a | n is said to be of n'th velocity order, which is commonly denoted O(n) in the literature 1 . For the metric g αβ , this expansion may be written explicitly in the form where we used overscripts to indicate velocity orders n g αβ ∼ O(n). Terms beyond the fourth velocity order are usually not considered in the standard PPN formalism implemented by xPPN. The zeroth velocity order is assumed to be the flat Minkowski background, For the metric perturbations, one finds that the first velocity order vanishes, 1 g αβ = 0, since the lowest velocity order of the matter source is the second order, as we will argue below. Further, the components 2 g 0a , 3 g 00 , 3 g ab , 4 g 0a change their sign under time reversal, and are prohibited by energy-momentum conservation. It follows that only the components 2 g 00 , 2 g ab , 3 g 0a , 4 g 00 , 4 g ab (4) are non-vanishing. For the first four terms, a particular standard gauge is assumed, in which the metric components take the form For the last term 4 g ab no such expansion is assumed in the standard PPN formalism, but it may be included in an extended formalism [13,14]. The terms on the right hand side are the so-called PPN potentials and PPN parameters; see sections IV F and IV G for their definition and [5] for a physical explanation. Here the PPN potentials describe the matter distribution, and their form is independent of the theory under consideration, while the PPN parameters are independent of the matter distribution, and their values are determined by the gravity theory. In order to determine the values of the PPN parameters for a given gravity theory, one follows a well-defined procedure to expand the gravitational field equations into velocity orders, which are then solved successively, up to the fourth order. The virtue of the PPN formalism is the fact that this form of the metric is sufficiently generic to solve the metric field equations of a large number of gravity theories, in which the source of gravity is the matter energy-momentum. In order to obtain this solution, one must also decompose the energy-momentum tensor (1) into space and time components, as well as into velocity orders. For the matter variables, one assigns the orders ρ ∼ Π ∼ O(2) and p ∼ O(4), based on their order of magnitude in the solar system. Lowering the indices of the energy-momentum tensor, its components then take the form Further, one assumes that the gravitational field is quasi-static, which means that it changes only due to the motion of the source matter, excluding, e.g., transient gravitational waves. Hence, time derivatives of all physical quantities are weighted with an additional velocity order, ∂ 0 ∼ O(1). In particular, this affects derivatives of the metric, which enter the field equations through the Levi-Civita connection and its curvature tensor. Here and in the following, we denote terms related to the Levi-Civita connection with an empty circle, in order to distinguish them from other connections to be introduced later, following the standard convention in the literature on teleparallel and related gravity theories, where this distinction becomes relevant; note that this circle is a distinct symbol from a zero denoting the zeroth velocity order.

B. Tetrad extension
There exist numerous gravity theories in which one of the fundamental fields is a tetrad (or coframe) field θ Γ α , which defines the metric via the relation where η Γ∆ = diag(−1, 1, 1, 1) is again the Minkowski metric, here with Lorentz indices. Its post-Newtonian expansion can be obtained as follows [15][16][17]. Its zeroth order is given by a diagonal background tetrad, defined in the previous section. For any higher order terms k θ Γ α in its perturbative expansion, it turns out to be more convenient to work with the expressions [18,19] which carry only spacetime indices, where is the inverse background tetrad. While it is entirely possible to use the perturbations k τ αβ as fundamental variables, it turns out to be more convenient to decompose them into metric perturbations and another, independent degree of freedom. This can be achieved by noting that the metric perturbations follow from the relation (8) to be given by Hence, the n'th order metric perturbation encodes n g αβ the symmetric part 2 n τ (αβ) of the tetrad perturbation, up to lower order terms. Isolating the highest orders from the sum on the right hand side, we find that this symmetric part is given by For the antisymmetric part, which constitute the aforementioned independent degree of freedom, one may define another, antisymmetric tensor field by In summary, the perturbative orders n θ Γ α of the tetrad for n ≥ 1 are then defined as in terms of n g αβ and n a αβ . Following a similar line of argument as for the metric, the only non-vanishing components of the antisymmetric tensor up to the fourth velocity order which must be considered are given by Note finally that the tetrad θ Γ α comes also with an inverse, the frame field e Γ α . It follows from its relation with the coframe field that its background is the diagonal element 0 e Γ α = ∆ Γ α , while higher perturbation orders are recursively defined as where the tetrad perturbations on the right hand side are further expanded using the rule (15).

C. Teleparallel geometry
The tetrad and its perturbative expansion discussed above play a fundamental role as the dynamical field variable in a particular class of gravity theories known as teleparallel gravity [20]. In their covariant formulation [21], another fundamental field variable is used, which is the spin connection • ω Γ ∆α . Together with the tetrad and its inverse, it defines another connection, whose coefficients are given by 2 2 Here we follow the convention used by xTensor, in which the first lower index α on the connection coefficients • Γ γ αβ is the "derivative" index.
The spin connection is further restricted to be flat and antisymmetric, so that also the affine connection has vanishing curvature and is metric-compatible, but it possesses non-vanishing torsion. Under these conditions, the spin connection turns out to be a gauge quantity, so that one can locally choose a Lorentz gauge in which it vanishes, known as the Weitzenböck gauge [22]. Choosing this gauge, • ω Γ ∆α ≡ 0, the connecting coefficients (18) are expressed in terms of the tetrad and its inverse alone [18]. This allows to derive their perturbative expansion, and the perturbative expansion of any derived tensor fields, in terms of the perturbations (15) and (17), and hence in terms of n g αβ and n a αβ .

D. Symmetric teleparallel geometry
Finally, the third class of gravity theories discussed here has become known as symmetric teleparallel gravity theories [23]. In these theories yet another connection is employed, whose coefficients we denote by × Γ γ αβ , and which is assumed to have vanishing curvature and torsion, It follows from these properties that there exists a local coordinate system (x ′α ), called the coincident gauge, such that its connection coefficients × Γ ′γ αβ in this coordinate system vanish [24]. Hence, in the standard PPN coordinate system (x α ), the connection coefficients are given by The post-Newtonian expansion of these connection coefficients is derived from the assumption that the two coordinate systems (x α ) and (x ′α ) are related by an infinitesimal coordinate transformation, induced by a vector field ξ α , so that up to the quadratic order one can write From this assumption follows that the connection coefficients take the form also up to the quadratic order in the vector field ξ α . Finally, the vector field is expanded in velocity orders. It turns out that in the PPN formalism the only non-vanishing components are given by [25] 2 ξ a , From these components the perturbative expansion of the connection coefficients × Γ γ αβ is obtained.

III. MATHEMATICAL FOUNDATION
Two fundamental mathematical concepts form the basis of the PPN formalism detailed in the previous section, and so their implementation is an important part of any computational tool to address this formalism. Here we briefly discuss these concepts and their formulation we choose to implement in xPPN. This will serve as an explanation of the mathematical background for the following sections. We discuss the split of tensors on spacetime into their space and time parts in section III A and the perturbative expansion into velocity order in section III B.

A. 3 + 1 split of tensors
As explained in the previous section, one of the crucial steps in applying the PPN formalism is the 3+1 decomposition of all tensorial quantities. A proper, geometric interpretation of this split can be given by regarding the spacetime manifold M 4 , equipped with coordinates (x α ) = (t, x a ), as a product manifold M 4 ∼ = T 1 × S 3 , where t is the time coordinate on the manifold T 1 ∼ = R, while (x a ) are the spatial coordinates on S 3 . It follows that for every t ∈ T 1 , one can find a map i t : S 3 → M 4 , (x a ) → (t, x a ), whose image in M 4 is called the spatial slice, or constant time hypersurface at time t. The collection of all spatial slices then forms a foliation of M 4 .
While the construction appears abstract, it gives a clear interpretation of the 3 + 1 split of tensor fields as follows. A vector field with components X α on M 4 splits in the form into a temporal part X 0 ∂ 0 , which is parallel to the coordinate lines generated by the time coordinate t, and a spatial part X a ∂ a , which is tangent to the spatial hypersurfaces. Under purely spatial coordinate transformations, X 0 behaves as a scalar, while X a are the components of a vector field. For any fixed value of t, we may thus take the pullback along i t to obtain a scalar function i * t X 0 , as well as a vector field i * t (X a ∂ a ), both now being tensor fields on S 3 , hence functions of the spatial coordinates (x a ). The dependence of the components X α on the time coordinate t is still present as a dependence on the choice of the spatial hypersurface, and t is regarded as a parameter which corresponds to this choice, instead of a coordinate.
The interpretation of time t as a parameter instead of a coordinate is, of course, purely mathematical, but it allows for a number of simplification when it comes to the implementation of the 3 + 1 split in computer algebra. Most importantly for the PPN formalism, it allows to treat time derivatives ∂ 0 essentially differently from spatial derivatives ∂ a . This is necessary, because in the PPN formalism time derivatives are weighted with an additional velocity order, which is not the case for spatial derivatives. Another advantage becomes apparent when considering tensors with symmetries. For example, considering an antisymmetric tensor field A αβ , the 3 + 1 split can be written as where A 00 = 0 and A 0a = −A a0 . Hence, the independent components can be described by a single covector field A 0a and antisymmetric tensor field A ab on S 3 . Thus, the use of tensor symmetries in conjunction with the 3 + 1 split allows to immediately discard vanishing components, which helps to simplify the calculation. The outlined interpretation of the 3 + 1 split and the decomposition of tensor fields on spacetime M 4 into tensor fields on space S 3 , with an additional parameter dependence, and corresponding counting of velocity orders, are a central ingredient of xPPN. In section V A we explain how to select certain components from the 3 + 1 decomposition of a tensor field, while section V C shows how to decompose arbitrary tensorial expressions. The implementation of these operations is outlined in appendix A 1, which shows the internal representation of decomposed tensor fields, while in section A 2 we schematically explain the decomposition algorithm and its treatment of tensor symmetries.

B. Perturbative expansion in velocity orders
Another important ingredient of the PPN formalism is the perturbative expansion of tensor fields into velocity orders. Similarly to the expansion (2) for the metric, also any other tensor field A is expanded as where each term is of the corresponding velocity order n A ∼ O(n), and where we omit tensor indices for brevity. Note that this is different from the convention used, e.g., in the xPert package [12], where the perturbative expansion takes the form and thus has the explicit form of a Taylor expansion in the perturbation parameter ε. We therefore list a few formulas which follow from the perturbative expansion (27). Most notably, for a tensor A = A 1 · · · A N given as a product of N tensors A 1 , . . . , A N , one has the n'th order term given by where the orders k i of the factors A i run over non-negative integers. Another important relation is the expansion of expressions F = f (A), where f is a scalar, or more general a tensorial function. In this case a Taylor expansion of the function f around the background 0 A is used, where f (k) denotes the k'th derivative, and then the product formula (29) is applied. For a function of N arguments, this formula naturally generalizes to Finally, time derivatives are weighted with an additional velocity order, and so for A ′ = ∂ 0 A one has and analogously for higher derivative orders, where the left hand side is to be understood as the n'th order term in the expansion of A ′ . In xPPN, the perturbative expansion of tensor fields into velocity orders, which takes into account the rules outlined above, is implemented in the function VelocityOrder, which is explained in detail in section V D. A few notes on the implementation of this function are given in appendix A 3.

IV. GEOMETRIC OBJECTS DEFINED BY XPPN
The PPN formalism relies on the existence of a number of objects, such as the spacetime manifold, the physical metric and its background value, as well as a particular set of potentials and constant parameters, which are necessary for its application to any gravity theory. For convenience, these objects are already defined by xPPN when the package is loaded, so that the user must define only those additional geometric objects which are specific to the gravity theory under consideration. Here we give an overview of these pre-defined objects. The manifolds corresponding to the split of spacetime into space and time, as well as a number of bundles over these manifolds, are given in section IV A. To define tensors on these manifolds and write tensorial expressions, xAct relies on indices; section IV B lists the indices defined by xPPN for the given bundles. The geometric objects constituting the background geometry are listed in section IV C, while section IV D gives a detailed account of the dynamical geometry and its perturbative expansion. The energy-momentum tensor and its expansion in fluid variables is shown in section IV E. Finally, the objects constituting the standard PPN metric are the PPN potentials detailed in section IV F and the PPN parameters listed in section IV G.

A. Manifolds and bundles
Tensors in xTensor are defined with respect to a manifold. In order to implement the 3 + 1 decomposition discussed in section III A, xPPN defines three manifolds:

2.
MfSpace is the three-dimensional space manifold S 3 . All tensors for which the 3 + 1 split into time and space components has been carried out are defined on this manifold, with an additional dependence on the time parameter TimePar.

3.
MfSpacetime is the four-dimensional spacetime manifold M 4 ∼ = T 1 × S 3 . It is defined as a product manifold, so that sums over tensor indices of M 4 may automatically be decomposed into sums over T 1 and S 3 . All geometric objects which appear in the theory under consideration should be defined on this manifold.
Each of these manifolds is canonically equipped with a tangent bundle; xTensor defines these automatically, and they are named TangentMfTime, TangentMfSpace and TangentMfSpacetime, respectively, and denoted by prefixing the manifold symbol with T. In addition, xPPN defines another vector bundle over each of these three manifolds, whose rank is the same as the dimension of the manifold, and which is used to define quantities carrying Lorentz indices. These bundles are named LorentzMfTime, LorentzMfSpace and LorentzMfSpacetime, respectively, and denoted by prefixing the manifold symbol with L.

B. Indices
Indices play an important role in xTensor, since they are used to denote tensor expressions and establish the relation between tensor slots and vector bundles. Often a large number of indices is necessary for longer tensor expressions, and each of them must be defined as a symbol (possibly with an appended number). xPPN addresses this problem by pre-defining a large number of indices using symbols which are unlikely to collide with other notation. In particular, the following indices are defined: is used by xPPN to denote the time component of the 3 + 1 decomposed forms of tensors, and is printed as 0. From the viewpoint of xTensor, this is a label index, which is not associated with any vector bundle and has no implied meaning. It is used because 3 + 1 decompositions are defined on S 3 , and so only spatial indices can be associated with the tangent bundle. Also it may appear an arbitrary number of times and in any position in any tensor expression.

\[ScriptT]
(printed as t ) is the generic index on the tangent bundle TT 1 . xTensor will use this to automatically generate as many indices t 1 , t 2 , . . . as needed as an intermediate step when performing a 3 + 1 decomposition of indices, before replacing them with 0.

\[ScriptCapitalT]
(printed as T ) is the generic index on the Lorentz bundle LT 1 . It is used by xTensor in the same way as t on TT 1 .

Lowercase
Latin letters a, . . . , z (character codes 97-122) are used to denote indices on TS 3 . xPPN automatically defines the symbols T3a, . . . , T3z to denote these indices, so that there is no need for the user to manually declare any indices, while at the same time avoiding possible name conflicts and leaving single letters free for other use. The prefix "T3" is omitted when the symbol is printed in output form, so that only the index letter itself appears in printed output. To obtain a list of pre-defined indices for any of these bundles, the xTensor command IndicesOfVBundle may be used. For example, to list all indices defined on TM 4 , use:

Uppercase Latin letters
The Mathematica command Names, together with a string pattern, may also be used to list these indices. This will give a similar result as the aforementioned command: C. Background geometry xPPN automatically defines a number of geometric objects corresponding to the background geometry, around which the perturbative PPN expansion is performed. On each of the three manifolds T 1 , S 3 , M 4 a metric, a tetrad and an inverse tetrad are defined. They are denoted as given in table I. Each of these geometric objects is defined to be constant with respect to the partial derivative PD on its respective manifold, so that PD applied to any of these objects vanishes. Further, contractions of a tetrad and its inverse are carried out automatically, and yield the corresponding delta tensor. NB! Note that each of the background metrics is defined as the first metric on its respective manifold. Hence, it is also used by xTensor in order to raise and lower indices, such that for a vector defined as A α one has A α A α = A α A β η αβ , which can easily be verified by using the xTensor command SeparateMetric: In order to raise and lower indices with the physical metric g αβ , the metric (as defined in the following section) must be written out explicitly.

D. Dynamical geometry
The background geometry detailed above is defined independently of any gravity theory and matter configuration. This is contrasted with the dynamical geometry which we discuss next. These are the dynamical fields which are used to mediate the gravitational interaction, and which are related to the matter source by the gravitational field equations. The most important of these geometric objects is the metric, which we discuss in section IV D 1. A similar role may be taken by the tetrad, as shown in section IV D 2. Further, three covariant derivatives are defined, which represent the three connections used in different formulations of general relativity [26]: the Levi-Civita connection in section IV D 3, the teleparallel connection in section IV D 4 and the symmetric teleparallel connection in section IV D 5.

Metric
As detailed in section II A, the central object in the PPN formalism is the metric g αβ , which is denoted by Met[-T4α, -T4β] in xPPN. Its inverse g αβ is written as InvMet[T4α, T4β]. Also here it must be emphasized that this is different from Met[T4α, T4β], which would be interpreted differently: A number of rules are defined for the perturbative expansion of the metric, and automatically applied when this expansion is performed. The zeroth order reduces to the background Minkowski metric (3), and so its space and time components are given by The remaining components, up to the fourth velocity order, are set to vanish, except for the components (4), for which no further rules are defined. These must be solved for in order to determine the PPN parameters. See section V B for more information how to apply these rules using the function ApplyPPNRules.

Tetrad
In order to implement the tetrad extension to the PPN formalism detailed in section II B, xPPN defines the tetrad θ Γ α as the object Tet[L4Γ, -T4α], together with an appropriate set of rules to evaluate its perturbative expansion using the formula (15). The antisymmetric tensor field a αβ occurring in this expansion is denoted by Asym[-T4α, -T4β] in xPPN. Further, next to the tetrad, also its inverse e Γ α is defined in xPPN, which is entered as InvTet[-L4Γ, T4α]. A number of automatically applied rules are defined for these inverse tetrads, so that they are contracted with tetrads if possible, and derivatives are evaluated according to Further, perturbations of the tetrad are evaluated following the formula (17), together with expanding the tetrad perturbation using the expansion (15).

Levi-Civita connection
Together with the metric Met[-T4α, -T4β], xPPN defines its unique metric-compatible, torsion-free Levi-Civita connection, whose coefficients are defined from the metric by the well-known formula (7). The corresponding covariant derivative is entered as CD[-T4α], and it is written in postfix notation using a semicolon ";", as well as using the symbol • ∇ in prefix notation, to distinguish it from other covariant derivatives introduced later. xTensor automatically defines a number of tensors for any covariant derivative. Most notably are the Christoffel "tensors" ChristoffelCD[T4γ, -T4α, -T4β], which are defined as the difference between the connection coefficients • Γ γ αβ and the (vanishing) coefficients of a fiducial connection associated to the partial derivatives ∂ α with respect to the coordinates. Further, xTensor defines a number of curvature tensors, such as the Riemann, Ricci, Einstein and Weyl tensors. xPPN defines the corresponding perturbative expansions in terms of the metric perturbations for all tensor fields derived from the covariant derivative CD[-T4α].

Teleparallel connection
In order to work with teleparallel gravity theories, as discussed in section II C, xPPN defines also the teleparallel connection (18), along with its torsion tensor and their perturbative expansion. The corresponding covariant derivative is entered as FD[-T4α], and is denoted • ∇ in prefix notation and a bar "|" in postfix notation. Working in the Weitzenböck gauge implies that the connection coefficients ChristoffelFD[T4γ, -T4α, -T4β] take the simple form and their perturbative expansion is defined accordingly.

Symmetric teleparallel connection
Finally, to accommodate the PPN formalism for symmetric teleparallel theories of gravity shown in section II D, xPPN provides yet another pre-defined connection, which is characterized by vanishing torsion and curvature, but which is not compatible with the metric g µν . Its covariant derivative is entered as ND[-T4α] and denoted by the symbol × ∇ in prefix notation, and a hash "#" in postfix notation. The perturbative expansion (23) of its connection coefficients is expressed in terms of a vector field ξ α , which is defined under the name Xi[T4α] in xPPN. Its perturbative expansion is implemented such that the only non-vanishing components are the three components (24). Also for this connection, xAct automatically defines torsion and curvature tensors, but both are set to vanish identically. The only tensorial quantity describing this connection and its relation to the Levi-Civita connection is the nonmetricity, which is given by In xPPN, the nonmetricity is called NonMet[-T4α, -T4β, -T4γ], and its perturbative expansion is obtained from the definition above.

E. Energy-momentum variables
In the PPN formalism it is assumed that the source of gravity is the energy-momentum tensor Θ αβ of a perfect fluid (1), which defined by xPPN is a tensor on M 4 denoted by EnergyMomentum[-T4α, -T4β]. Note that it is defined with lower indices. In order to raise the indices, the physical metric g αβ must be used explicitly; implicitly raising the indices by writing EnergyMomentum[T4α, T4β] would use the background metric η αβ , and hence not give the correct result. Also note that we use the symbol Θ instead of the more conventional T in order to avoid confusion with the torsion tensor.
For convenience, xPPN also defines the trace-reversed energy-momentum tensor denoted by TREnergyMomentum[-T4α, -T4β], and which is likewise defined as a tensor on M 4 . In order to describe the 3 + 1 split of the energy-momentum tensor and its decomposition into velocity orders, xPPN further defines the following tensors on S 3 depending on TimePar, which represent the variables describing the perfect fluid: When the energy-momentum tensor Θ αβ is expanded in velocity orders as shown in section V D, the relations which follow directly from the expansion (6), are applied, while all other components vanish.

F. Post-Newtonian potentials
The post-Newtonian potentials are a central ingredient to the PPN formalism. In xPPN they are defined as tensors on the space manifold S 3 with an additional time dependence. Most of them are scalars, but there are also vector and tensor potentials, which therefore carry indices in the tangent bundle TS 3 . In particular, the following PPN potentials are defined: 1. PotentialChi[] is the post-Newtonian superpotential 2.
PotentialPhi2[] is the gravitational self-energy potential 8.
PotentialPhi3[] is the internal energy potential 9.
PotentialPhi4[] is the pressure potential 10.
PotentialA[] is the anisotropic kinetic potential 11. PotentialB[] is the potential related to change in velocity 12.
PotentialPhiW[] is the Whitehead potential Note that these integrals are not implemented directly in xPPN, which makes no explicit use of the coordinates apart from the assumption that partial derivatives of the tensors defining the background geometry vanish. However, they enter into the formalism indirectly, since they define the relations between the potentials and the source terms which are explained in detail in section V F.

G. Post-Newtonian parameters
The PPN parameters are represented in xPPN by objects which are declared as constants within xTensor, so that any derivatives acting on them vanish. A full list is given in table II. Note that assuming the PPN parameters to be constant, which is one of the standard assumptions of the PPN formalism, means that it cannot be directly applied to, e.g., theories with massive fields mediating the gravitational interaction, since in such theories the values of the PPN parameters depend on the distance between the source and the test mass; also time variation of PPN parameters depending on a cosmological background is not considered. The former may be included in a future version of xPPN by introducing additional Yukawa-type potentials, which depend on a mass determining the interaction scale [27,28]. The latter might be implemented by allowing for PPN parameters which depend on TimePar instead of being declared constant [29].

V. UTILITY FUNCTIONS
We now present a number of utility functions defined by xPPN, which can be used to manipulate terms which typically appear in the post-Newtonian expansion, perform necessary computational steps and solve the gravitational field equations in terms of the post-Newtonian potentials and parameters listed in the previous sections. Which of these functions need to be used highly depends on the gravity theory under consideration, and we show this for a few of them in an explicit example in section VI; here we give an overview of the functions provided and how they are applied. In section V A we show how to refer to specific parts of tensors in their space-time decomposition and perturbative expansion. The definition and application of substitution rules for such terms is explained in section V B. The decomposition of general tensorial expressions into space and time components is shown in section V C, and further into velocity orders in section V D. The transformation of terms using the Euler equations derived from energy-momentum conservation is displayed in section V E, and applied to post-Newtonian potentials in section V F. Finally, utility functions for sorting derivatives prior to applying these rules are shown in section V G.

A. Selecting space-time components and velocity orders of tensors
As explained in section III A, xPPN defines for each tensor declared on the spacetime manifold M 4 a number of tensors on the space manifold S 3 , which represent the 3 + 1 decomposition of the initial tensor and its velocity orders. In order to access these components, xPPN defines the utility function PPN to easily obtain them from a tensor head, a sequence of indices and optionally a velocity order. This function can be used in two ways, taking the following arguments: where h is a tensor head and i is a sequence of indices, such that each index either belongs to S 3 or is LI[0] (possible with a minus sign);

PPN[h, n][i]
, where n is a non-negative integer and h and i are as above.
The application of this function is illustrated by the following example of a vector A α : In

B. Definition and application of replacement rules
For pre-defined tensors on the spacetime manifold M 4 , such as the metric or the energy-momentum tensor, xPPN defines a number of substitution rules for the terms in their post-Newtonian expansion. In order to apply these rules to any given tensorial expression, xPPN provides the function ApplyPPNRules, which can be invoked in two different ways: 1. ApplyPPNRules[X] recursively applies the defined PPN rules to all tensors which appear in the expression X and its subexpressions.
2. ApplyPPNRules[X, h] applies the defined PPN rules only to those tensors within X with head h.
For example, the zeroth order of the physical metric is given by the Minkowski background metric: xPPN defines two functions handling the 3 + 1 decomposition of arbitrary tensorial expressions on spacetime into their space and time components, as explained in section III A: SpaceTimeSplit and SpaceTimeSplits. While the former calculates one specific component of the 3 + 1 decomposition, the latter yields an array with all components. In both functions the decomposition is applied to both free and dummy indices.
The function SpaceTimeSplit takes two arguments: a tensor on M 4 (possibly containing free indices) and a list of replacement rules, which assigns to each free index in a bundle over M 4 either an index over the corresponding bundle over S 3 or the label index LI[0] (possibly dressed with a minus sign to indicate a lower index). For example, for an expression of the form A α β with free indices T4α, -T4β, the second argument may be any of the following: Of course, other names than a, b may be used for the indices on S 3 , but their position must remain the same; their role is to specify how the free indices in the resulting expression will be named.
Observe that also the dummy index γ has been split into a sum over dummy indices in the 3 + 1 decomposition. The function SpaceTimeSplits is similar, but in its second argument only indices of S 3 (and hence no LI[0]) are allowed as the right hand sides of the rules. The function then returns an array in which the free indices of the original expression are replaced either by the specified spatial indices or the time index 0, as shown in the following example: The order of the rules in the replacement list corresponds to the order of the dimensions of the resulting array. In the example above, the first dimension corresponds to α, while the second dimension corresponds to β. Hence, selects the component where α is replaced by 0 (the first element of (0, a)), while β is replaced by b (the second element of (0, b)). Similarly, yields the transposed array. Finally, we remark that both SpaceTimeSplit and SpaceTimeSplits also operate on partial derivatives PD[-T4α] with respect to spacetime coordinates. These are converted as follows:

In[ ]:= % == {ParamD[TimePar][PPN[A][]], PD[-T3a][PPN[A][]]} Out []= True
Observe that the time derivative is not represented as a partial derivative, but as a parameter derivative with respect to the time parameter TimePar. Finally, we remark that the automatic split is implemented only for partial derivatives PD[-T4α]; any other derivatives most be converted with ChangeCovD or LieDToCovD, as appropriate.

D. Decomposition into velocity orders
As explained in section III B, an important part of the PPN formalism is the series expansion of expressions in velocity orders. In order to select a single term n X ∼ O(n) from an expression X, xPPN provides the function VelocityOrder. The term n X is then expressed as VelocityOrder[X, n]. Here n must be a non-negative integer, while X can be any tensorial expression on S 3 , which may in addition depend on TimePar. It makes use of a number of standard relations for the velocity order in the PPN formalism: orders are distributed over products, and time derivatives are weighted with an additional velocity order. This is illustrated in the following example: Observe that in the first case the rule 2 Θ 00 → ρ is applied, but not in the second case. For large expressions the default setting UsePPNRules →True is significantly faster, since the PPN rules imply the vanishing of many terms in the post-Newtonian expansion, which are thus immediately discarded when the rules are applied.

E. Euler equations
An important assumption of the PPN formalism is that the energy-momentum tensor Θ αβ introduced in section IV E satisfies the covariant conservation equation

TimePiToEuler[X] applies the replacement
The functions can be successively applied in order to transform terms involving higher order time derivatives.

F. Transformation of PPN potentials
The post-Newtonian potentials defined in section IV F and their derivatives satisfy a number of relations among each other and with the energy-momentum variables defined in section IV E, some of which follow directly from the definition of the potentials, while others can be derived from the Euler equations discussed in the previous section. xPPN defines a number of utility functions in order to implement these relations and use them to turn PPN potentials and their derivatives into either other PPN potentials or terms constructed from the energy-momentum tensor. The most important of these is PotentialToSource[X], which performs the following replacements on all subexpressions of X: The remaining functions act on specific PPN potentials or their derivatives. Their effects are listed in table III.

G. Sorting of derivatives
In xAct, derivatives are, in general, not sorted automatically. This poses two difficulties, which may lead to problems in simplifications: 1. Mathematica does not recognize terms which are mathematically equal if they contain derivatives in different order. For example, expressions of the form ∂ α ∂ β X and ∂ β ∂ α X are equal, since partial derivatives commute, but since they are represented differently, Mathematica does not recognize this.

Function Transformation
PotentialChiToU △χ → −2U PotentialChiToPhiAB 2. In order to apply the transformations listed in section V F, pattern matching is applied to find divergences, time derivatives and Laplace operators acting on tensors. However, due to the way how derivatives are represented in xTensor, Mathematica recognizes such terms only if they are not interspersed with other derivatives. For example, ∂ α A α is found in ∂ β ∂ α A α , but not in ∂ α ∂ β A α , even though these terms are mathematically equal.
The first of these problems can be solved by defining a canonical order for derivatives and keeping them always sorted in canonical order. Such a canonical ordering is implemented as part of the xTensor function ToCanonical, so that the following terms are recognized by Mathematica as equal and canceled:

In[ ]:= PD[-T4α][PD[-T4β][A[]]] -PD[-T4β][PD[-T4α][A[]]]
However, this does not solve the second problem, since a different order of derivatives is required, depending on the type of pattern to be matched; for example, matching a time derivative would require the order ∂ α ∂ 0 A α , while matching a divergence would require the order ∂ 0 ∂ α A α . Hence, keeping all indices in a fixed, canonical order would not allow for such patterns to be found. Hence, xPPN implements a different method by defining different functions, which allow sorting derivatives on demand in any of the necessary orders. The following functions are defined: 1. SortPDs[X] sorts all derivatives in canonical order by applying the following rules: (a) Spatial derivatives are applied before time derivatives, i.e., moved to the right.  The effect of these functions on an expression with multiple derivatives acting on a tensor is shown by the following code:

In[ ]:= expr = PD[T3b][ParamD[TimePar][PD[-T3a][PD[-T3c][PD[-T3b][PPN[A][T3c]]]]]]
The implementation of these functions is inspired by the field theory package xTras [30], which defines a similar set of functions, which take into account also non-commuting derivatives by adding suitable correction terms. Also note that the replacement functions listed in section V F make use of these sorting functions automatically in order to establish the necessary order of indices before pattern matching is performed.

VI. EXAMPLE: SCALAR-TENSOR GRAVITY
In order to demonstrate the use of the xPPN package, we provide a practical example in form of a complete, commented session, which calculates the PPN parameters of a scalar-tensor class of gravity theories, thereby showing the use of some of the functions detailed in the previous sections. The precise steps which are needed to determine the PPN parameters highly depend on the specific theory under consideration, and must be chosen accordingly. We briefly present the action and field equations of this class in section VI A. A few preliminary commands, for loading the package and setup, are given in section VI B. In section VI C, the necessary tensors and constants for the calculation are defined. These are used in the definition of the field equations in section VI D. Their post-Newtonian expansion is derived in section VI E. In section VI F, the perturbative solution of the resulting equations is obtained. Finally, in section VI G, the PPN metric and parameters are calculated.

A. Action and field equations
In the following we discuss a class of scalar-tensor theories of gravity, whose action is given by [31] in Brans-Dicke like parametrization in the Jordan conformal frame. Here S m denotes the matter part of the action, where we collectively denoted by χ the set of matter fields. The gravitational part contains a free function ω of the scalar field ψ. Each theory of this class is defined by a particular choice of this free function ω. By variation of this action with respect to the metric and the scalar field as well as subtraction of a suitable multiple of the trace of the metric field equation one obtains the field equations where ∇ ν is the d'Alembert operator on the physical spacetime M 4 .

B. Package loading and preliminaries
In order to load and use xPPN, as the first prerequisite a working installation of xAct [10] is needed. The files provided by xAct then reside in a directory named xAct in the Mathematica package search path. The xPPN package comes as a directory named xPPN, which must be placed inside the xAct directory. If both packages are installed correctly, xPPN can be loaded with the command This should load xPPN and its dependencies from the xAct package suite. Note that loading may take some time, since xPPN calculates the perturbative expansion of the pre-defined tensor fields upon package loading. To suppress $ symbols in the index notation, it is useful to set the following xAct printing option: In[ ]:= $PrePrint = ScreenDollarIndices; Finally, we define two utility functions which help to create rules from equations; this functionality, which is simply a shorthand notation for the versatile function MakeRule from the xTensor package, are not specific to xPPN, and are defined here only to shorten the notation in the later course of this example: We continue by defining the xAct objects which we will be using for the calculation and their notation. We start with the scalar field ψ, which is a tensor without any indices.
In[ ]:= DefScalarFunction[omega, PrintAs → "ω"] We then come to the field equations (56), which we write in the form E αβ = 0 and E = 0, by moving the energymomentum tensor to the left hand side. The two tensors representing these field equations are then defined as follows. Finally, in this example we will solve the field equations by making an ansatz for the solution in terms of PPN potentials and unknown, constant coefficients, which we will then determine. For this purpose, we define a function which creates such constant coefficients on demand as follows. In the next step, we enter the field equations (56). As mentioned before, we must pay attention to explicitly use the inverse metric g αβ for terms such as the kinetic term ∂ ρ ψ∂ ρ ψ of the scalar field or the trace Θ ρ ρ of the energy-momentum tensor. Taking these into account, we can enter the metric field equation (56a) as follows. The most crucial and resource intensive part of the calculation is the derivation of the post-Newtonian expansion of the field equations. For this purpose, we must first define the expansion of the scalar field. Here we impose the relations In xPPN, they are implemented as follows: In With these definitions in place, we can continue with the calculation. We will do so in two steps. First, we perform a 3 + 1 decomposition of the field equations. For this purpose, we convert the Levi-Civita covariant derivatives to partial derivatives and Christoffel symbols. Also we perform a number of simplifications. Again we start with the metric field equations.

F. Perturbative solution
We can now use the post-Newtonian expansion of the field equations and solve them order by order. For this purpose we make use of prior knowledge that the field equations of the theory at hand can be solved by a particular ansatz for the metric and scalar field, so that at every step of the calculation we are left with solving for a number of constant coefficients; this ansatz and solution method must be adapted if the package is applied to other theories. We show the zeroth velocity order (the background vacuum solution) in section VI F 1, the second velocity order in section VI F 2, the third velocity order in section VI F 3 and the fourth velocity order in section VI F 4.

Zeroth velocity order
First, we verify that the background vacuum equations are solved identically by the background geometry. This can be checked as follows. These equations take the form The easiest way to solve these equations is using an ansatz for the metric and scalar field perturbations in terms of the post-Newtonian potentials, with arbitrary, constant coefficients. This ansatz can be defined as follows: In  To obtain the equations to be solved, one inserts this ansatz into the field equations. Since the field equations contain derivatives of the metric and scalar field perturbations, this will yield derivatives of the post-Newtonian potentials U and U ab . These must be matched with the terms arising from the energy-momentum tensor Θ αβ . For this purpose, one applies the functions listed in section V F. For the potentials at hand it is most convenient to first convert them to derivatives of the superpotential χ, before transforming them to terms involving the matter density ρ. Inspecting the obtained equations shows that they are given by These equations are solved if and only if the coefficients of each of the terms ρ, ρδ ab and △χ ,ab vanishes individually. This yields four equations for the four coefficients a 1 , . . . , a 4 . However, one finds that these are not linearly independent, due to the gauge freedom arising from the diffeomorphism invariance of the theory. Hence, they must be supplemented with an additional, gauge fixing equation. The common choice in the PPN formalism is to eliminate the term U ab from the component 2 g ab , thus setting a 3 = 0. We can thus determine the component equations: This yields the solution , a 2 = κ 2 ω(Ψ) + 1 2πΨ(2ω(Ψ) + 3) , a 3 = 0 , a 4 = κ 2 4π(2ω(Ψ) + 3) .
We may check that this indeed solves the equations: In

Third velocity order
We then come to the third velocity order. Also here we proceed in full analogy to the second velocity order shown above. First, we isolate the field equation This equation takes the form We then choose an ansatz for the third-order metric perturbation To obtain the equation to be solved, we insert this ansatz into the field equation, together with the second-order solution obtained in the previous step. Also here we obtain derivatives acting on the post-Newtonian potentials, which we can simplify by using their interrelations, and finally convert them to terms matching the energy-momentum source. Here it is most useful to first eliminate the potential W a . The resulting terms can then be converted to source terms and terms involving derivatives acting on U only: Inspecting this equation shows the following form: This equation only determines the sum a 5 + a 6 of the coefficients we introduced. Here we leave their difference as an undetermined parameter, which we will solve for when we come to the fourth velocity order, which will fix the post-Newtonian gauge. We thus determine the solution: In The solution is given by We check that this solves the equations: Together with the ansatz for the third-order metric component Finally, we also check that this solves the third-order field equation we started with: In As we shall see below, this equation can be solved by the following ansatz: In In[ ]:= ans4ru = mkrg[ans4def]; We then insert this ansatz into the fourth-order field equation, alongside the previously determined solutions at lower velocity order. In order to convert the post-Newtonian potentials to a unified form, from which we can read off the independent equations for the constant coefficients in the ansatz, it is sufficient to replace the divergence terms V a,a and W a,a by time derivatives of U , before replacing all potentials by source terms. This is done as follows: This yields the solution We check that it is indeed a solution: This takes the form To solve this equation, we take its positive root: This yields the solution With this solution, we can now determine the PPN parameters. These are the equations to be solved. To determine the PPN parameters, we isolate the constant coefficients in front of the post-Newtonian potentials This is of course the well-known post-Newtonian limit of scalar-tensor gravity with a massless scalar field [31].

VII. SUMMARY AND OUTLOOK
We have presented the Mathematica package xPPN, which is based on the tensor algebra suite xAct and which implements the PPN formalism in its formulation presented in [5], with extensions to adapt it to the geometries employed in the three formulations of general relativity [26] and modifications thereof. Besides discussing its underlying mathematical concepts and giving a detailed account on its usage, we provided an example application to a simple scalar-tensor theory of gravity. We believe that this package will find wide application to assess the viability of gravity theories by comparing their post-Newtonian limit to observations in the solar system.
Despite the generality of xPPN, being applicable to gravity theories based on different geometric foundations, various extensions and modifications are possible and may be implemented in future versions. For example, one may consider the following types of extensions: 1. Alternative formulations of the PPN formalism: A newer approach towards the PPN formalism employs a different density variable [7,Sec. 4], which has the advantage that it simplifies the gauge transformation of the PPN potentials compared to the original formulation. Another approach to simplify the issue of gauge transformations is to resort to a formulation which makes use of gauge-invariant perturbation theory, thereby resolving the necessity of gauge considerations and simplifying the application of the PPN procedure [32]. Such modifications may be included by changing the expansion of the metric in terms of PPN potentials.
2. Additional post-Newtonian potentials and parameters: Various theories of gravity exhibit a post-Newtonian limit in which the metric perturbation cannot be expressed only in terms of the post-Newtonian potentials used in the standard formalisms, and so more general potentials and corresponding parameters have been introduced. The presence of massive fields leads to the appearance of Yukawa-type potentials [27,28], while higher-order derivatives can be accommodated by further integrals in the post-Newtonian potentials [33,34]. Another class of terms arises from theories which include parity-violating contributions [35][36][37], or by violation of diffeomorphism invariance [38][39][40]. Further, one may introduce also fourth-order tensor potentials to expand the term 4 g ab , which is neglected in the standard PPN formalism, but may be used to describe higher-order corrections to light propagation [13,14]. Such additional potentials may easily be included by defining the corresponding tensor fields and the equations which relate them to the corresponding source terms, in full analogy to the standard PPN potentials.
3. Higher than fourth order in the post-Newtonian expansion: While in the standard PPN formalism implemented in the present version of xPPN the metric (and possible other fields present in the theory) are expanded only up to the fourth velocity order, one may also consider higher order terms, and address problems such as the post-Newtonian dynamics of systems involved in the generation of gravitational waves [7,41]. Care must be taken since the Euler equations governing the dynamics of the source matter attain dissipative terms at higher orders, and also the gravitational field exhibits dissipative behavior due to the loss of energy by emitted gravitational waves. Further, such as extension requires the definition of higher-order post-Newtonian potentials, along with corresponding post-Newtonian parameters, whose relation to observations and associated observables must be defined, so that such an extension would be a major task.
4. Generalized post-Newtonian expansion: The standard PPN formalism assumes the validity of a post-Newtonian expansion around a flat Minkowski background. This assumption may be generalized, by allowing for a timedependent Friedmann-Lemaître-Robertson-Walker background [29], or by including non-perturbative effects arising from the Vainshtein screening mechanism [42].
5. More general geometric frameworks: Another assumption of the standard PPN formalism states that the motion of test masses is governed by a single metric, and so the observable effects on test masses are fully covered by a post-Newtonian expansion of this single metric. This assumption may be relaxed if there are several dynamical metrics present in the considered gravity theory, which lead to more complex test matter dynamics and generalized sources of gravity [43][44][45]. Also more general connections may be considered, such as a general teleparallel connection [46], as well as Riemann-Cartan [16,33,34,47,48] or metric-affine geometry [49]. In theories of this class additional matter currents may appear from a coupling of spin to gravity, which gives rise to additional source terms, which must be accommodated for by further PPN potentials, together with their respective equations of motion and conservation laws generalizing the Euler equations [50][51][52][53][54][55][56][57].

ACKNOWLEDGMENTS
The author thanks Yuri Obukhov and Dirk Puetzfeld for helpful comments on the manuscript and pointing out related work. He gratefully acknowledges the full support by the Estonian Research Council through the Personal Research Funding project PRG356, as well as the European Regional Development Fund through the Center of Excellence TK133 "The Dark Side of the Universe".

Appendix A: Implementation notes
Although xTensor offers a possibility to split tensor indices for product manifolds and to project tensors to submanifolds, while xPert provides an implementation of a perturbative expansion of tensor fields, these approaches are not the most well-adapted to the mathematical concepts discussed in section III. While it may still be possible to use these packages for the task at hand, xPPN provides its own implementations of these concepts, in order to match as closely as possible with the PPN formalism, and to simplify the different perturbative treatment of space and time components of tensors and in particular their derivatives. This appendix contains a few notes how the functions detailed in the main part of this article are implemented. The internal representation of split tensor fields is shown in section A 1. In section A 2, we give an overview of the algorithm used to obtain this decomposition for arbitrary tensor fields. Finally, in section A 3 we give a few notes on the implementation of the perturbative expansion in velocity orders. The aim of this technical appendix is to give an insight to the internal data structures and algorithms used by 4. SymmetryGroupOfTensor returns the remaining symmetry group of the 3 + 1 decomposed tensor. 5. PrintAs yields the same symbol PrintAs[h] which is used for printing h if no velocity order is given, and otherwise adds the velocity order n as an overscript n h.
Further, PPNTensor automatically applies rules which are obtained from index symmetries. These are calculated whenever a tensor field is defined, as explained below.

3 + 1 spacetime split algorithm
In order to retain only independent components in the 3 + 1 split of a tensor, xPPN examines tensor symmetries whenever a new tensor on the spacetime manifold is defined. It does so by extending the xTensor function DefTensor, using the extensibility framework implemented in xAct, to perform the following steps: 1. Determine which index slots of the tensor are associated to the spacetime manifold. If a tensor carries additional, internal indices, then they will be treated separately and ignored in the following steps.
2. Calculate the symmetry group acting on the spacetime indices only, i.e., factor out any possible symmetries acting on internal indices and also ignore them in the following steps.
3. Create a list of all possible ways to split the spacetime indices into space and time components.
4. Consider the action of the symmetry group on the splits from the previous step and determine the orbits of this action. For example, if a tensor carries a symmetry (or antisymmetry) in two indices, then there is no difference between choosing the first index to be temporal and the second index spatial or vice versa, and so these two possible splits belong to the same orbit; these two possible ways to arrange the indices of the split tensor are equivalent.
5. For each orbit, perform the following steps: (a) Pick a representative, i.e., from all equivalent ways to arrange the indices, choose the one which comes first in canonical (lexicographic) order.
(b) Check whether any of the elements of the subgroup which permutes only the temporal indices comes with an antisymmetry, i.e., changes the sign of the tensor. If this is the case, this component must vanish identically, since permuting only temporal indices must leave the tensor invariant, and so it must be equal to its negative. In this case, define all index combinations in this orbit to be an alias for the zero tensor Zero. (c) If the chosen representative does not vanish, conclude with the following steps: i. Determine the remaining symmetry group acting on the spatial indices only. ii. Define a tensor on a purely spatial manifold, whose temporal and spatial indices match with the representative, which carries the symmetry determined in the previous step in the spacetime indices and any inherit symmetry in internal indices, and which in addition depends on a time parameter.
iii. Define all other arrangements of temporal and spatial indices which are equivalent to the representative as alias for the previously defined tensor or its negative, taking into account possible sign changes when indices are permuted.
To illustrate these steps, consider the case of defining an antisymmetric tensor A αβ = A [αβ] on the spacetime manifold, by invoking Then the following steps are carried out automatically by xPPN : 1. There are no internal indices, and so the spacetime related symmetry operations acting on the tensor indices are the identity and the swapping of the two indices, where the latter also changes the sign of the tensor.
2. Each of the two indices of A αβ splits into time and space. Hence, there are four possible decomposed tensor fields: A 00 , A 0a , A a0 , A ab .
3. The two components A 0a and A a0 are transformed into each other under the index symmetry of the original tensor fields; hence, they belong to the same orbit of the action of the symmetry group. The remaining components are the sole elements in their respective orbits.
4. The three orbits are examined separately: (a) For A 00 one finds that the swap operation acts on temporal indices only, but also changes the sign of the tensor. Hence, A 00 = −A 00 = 0 and the tensor is defined as an alias of Zero: (b) For A 0a and A a0 the former arrangement of indices is chosen as a representative, since the indices are in canonical (lexicographic) order. A purely spatial tensor A 0a is defined, which has one free index slot. A a0 is defined as an alias for −A 0a : (c) Another spatial tensor A ab is defined, which is antisymmetric in its two index slots. The tensor symmetry is again associated to the symbol A.
Note that all properties of these tensors, such as their symmetries and tensor slots, are associated with the symbol (tensor head) which is used in the original call to DefTensor to define the tensor on the spacetime manifold, and no further symbols are introduced.

Velocity order decomposition algorithm
In order to implement the rules for the perturbative expansion in velocity orders detailed in section III B, a few special cases have been defined for the function VelocityOrder shown in section V D. For the product rule (29), it is most convenient to rely on Mathematica's pattern matching and recursive function application functionality. Given a product A = A 1 · · · A N of N tensors, one may split off the first factor, and then recursively apply the rule This is repeated until only a single factor is left. The implementation then takes the following simple form. Next, we come to the relation (31) for the perturbative expansion of expressions which are given as functions of an arbitrary number of arguments. To evaluate this formula, it is useful to define the formal series for the tensor fields A 1 , . . . , A N . Observe that the n'th order series coefficient 1 n! d n dǫ n f (Ã 1 (ǫ), . . . ,Ã N (ǫ)) ǫ=0 = p ik f (p11+...+p1n,...,pN1+...+pNn) 0 A 1 , . . . , where the sum runs over p ik ∈ {0, . . . , q k } , is exactly the term of n'th velocity order in the expansion (31). Hence, it can be obtained as follows.