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

There exist different versions of the PPN formalism. Here we adhere to its form presented in Will (1993). Its starting point is the assumption that the propagation of light and massive test bodies is governed by a pseudo-Riemannian metric $g_{\alpha\beta}$ 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

$$\Theta^{\alpha\beta}=(\rho+\rho\Pi+p)u^{\alpha}u^{\beta}+pg^{\alpha\beta}$$

(1)

of a perfect fluid with four-velocity $u^{\alpha}$, rest energy density $\rho$, specific internal energy $\Pi$ and pressure $p$. Further, a fixed Cartesian coordinate system $(x^{\alpha})=(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}|\ll c\equiv 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 $\mathcal{O}(n)$ in the literature¹¹1Note in particular that $\mathcal{O}(1)\sim|v|$ in this notation is the first velocity order, and not unity, which would be $\mathcal{O}(0)\sim 1$.. For the metric $g_{\alpha\beta}$, this expansion may be written explicitly in the form

$$g_{\alpha\beta}=\sum_{n=0}^{\infty}\accentset{n}{g}_{\alpha\beta}=\accentset{0}{g}_{\alpha\beta}+\accentset{1}{g}_{\alpha\beta}+\accentset{2}{g}_{\alpha\beta}+\accentset{3}{g}_{\alpha\beta}+\accentset{4}{g}_{\alpha\beta}+\mathcal{O}(5)\,,$$

(2)

where we used overscripts to indicate velocity orders $\accentset{n}{g}_{\alpha\beta}\sim\mathcal{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,

$$\accentset{0}{g}_{\alpha\beta}=\eta_{\alpha\beta}=\mathrm{diag}(-1,1,1,1)\,.$$

(3)

For the metric perturbations, one finds that the first velocity order vanishes, $\accentset{1}{g}_{\alpha\beta}=0$, since the lowest velocity order of the matter source is the second order, as we will argue below. Further, the components $\accentset{2}{g}_{0a},\accentset{3}{g}_{00},\accentset{3}{g}_{ab},\accentset{4}{g}_{0a}$ change their sign under time reversal, and are prohibited by energy-momentum conservation. It follows that only the components

$$\accentset{2}{g}_{00}\,,\quad\accentset{2}{g}_{ab}\,,\quad\accentset{3}{g}_{0a}\,,\quad\accentset{4}{g}_{00}\,,\quad\accentset{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

	$\displaystyle\accentset{2}{g}_{00}$	$\displaystyle=2U\,,$		(5a)
	$\displaystyle\accentset{2}{g}_{ab}$	$\displaystyle=2\gamma U\delta_{ab}\,,$		(5b)
	$\displaystyle\accentset{3}{g}_{0a}$	$\displaystyle=-\frac{1}{2}(3+4\gamma+\alpha_{1}-\alpha_{2}+\zeta_{1}-2\xi)V_{a}-\frac{1}{2}(1+\alpha_{2}-\zeta_{1}+2\xi)W_{a}\,,$		(5c)
	$\displaystyle\accentset{4}{g}_{00}$	$\displaystyle=-2\beta U^{2}+(2+2\gamma+\alpha_{3}+\zeta_{1}-2\xi)\Phi_{1}+2(1+3\gamma-2\beta+\zeta_{2}+\xi)\Phi_{2}$
		$\displaystyle\phantom{=}+2(1+\zeta_{3})\Phi_{3}+2(3\gamma+3\zeta_{4}-2\xi)\Phi_{4}-2\xi\Phi_{W}-(\zeta_{1}-2\xi)\mathcal{A}\,.$		(5d)

For the last term $\accentset{4}{g}_{ab}$ no such expansion is assumed in the standard PPN formalism, but it may be included in an extended formalism Richter and Matzner (1982a, b). The terms on the right hand side are the so-called PPN potentials and PPN parameters; see sections IV.6 and IV.7 for their definition and Will (1993) 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 $\rho\sim\Pi\sim\mathcal{O}(2)$ and $p\sim\mathcal{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

	$\displaystyle\Theta_{00}$	$\displaystyle=\rho\left(1+\Pi+|v|^{2}-\accentset{2}{g}_{00}\right)+\mathcal{O}(6)\,,$		(6a)
	$\displaystyle\Theta_{0a}$	$\displaystyle=-\rho v_{a}+\mathcal{O}(5)\,,$		(6b)
	$\displaystyle\Theta_{ab}$	$\displaystyle=\rho v_{a}v_{b}+p\delta_{ab}+\mathcal{O}(6)\,.$		(6c)

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, $\partial_{0}\sim\mathcal{O}(1)$. In particular, this affects derivatives of the metric, which enter the field equations through the Levi-Civita connection

$$\accentset{\circ}{\Gamma}^{\gamma}{}_{\alpha\beta}=\frac{1}{2}g^{\gamma\delta}(\partial_{\alpha}g_{\delta\beta}+\partial_{\beta}g_{\alpha\delta}-\partial_{\delta}g_{\alpha\beta})$$

(7)

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.

There exist numerous gravity theories in which one of the fundamental fields is a tetrad (or coframe) field $\theta^{\Gamma}{}_{\alpha}$, which defines the metric via the relation

$$g_{\alpha\beta}=\eta_{\Gamma\Delta}\theta^{\Gamma}{}_{\alpha}\theta^{\Delta}{}_{\beta}\,,$$

(8)

where $\eta_{\Gamma\Delta}=\mathrm{diag}(-1,1,1,1)$ is again the Minkowski metric, here with Lorentz indices. Its post-Newtonian expansion can be obtained as follows Nitsch and Hehl (1980); Smalley (1980); Hayward (1981). Its zeroth order is given by a diagonal background tetrad,

$$\accentset{0}{\theta}^{\Gamma}{}_{\alpha}=\Delta^{\Gamma}{}_{\alpha}=\mathrm{diag}(1,1,1,1)\,.$$

(9)

defined in the previous section. For any higher order terms $\accentset{k}{\theta}^{\Gamma}{}_{\alpha}$ in its perturbative expansion, it turns out to be more convenient to work with the expressions Ualikhanova and Hohmann (2019); Hohmann (2021)

$$\accentset{k}{\tau}_{\alpha\beta}=\eta_{\alpha\gamma}\Delta_{\Gamma}{}^{\gamma}\accentset{k}{\theta}^{\Gamma}{}_{\beta}\,,$$

(10)

which carry only spacetime indices, where

$$\Delta_{\Gamma}{}^{\alpha}=\mathrm{diag}(1,1,1,1)$$

(11)

is the inverse background tetrad. While it is entirely possible to use the perturbations $\accentset{k}{\tau}_{\alpha\beta}$ 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

$$\accentset{n}{g}_{\alpha\beta}=\eta_{\Gamma\Delta}\sum_{k=0}^{n}\accentset{k}{\theta}^{\Gamma}{}_{\alpha}\accentset{n-k}{\theta}^{\Delta}{}_{\beta}=\eta^{\gamma\delta}\sum_{k=0}^{n}\accentset{k}{\tau}_{\gamma\alpha}\accentset{n-k}{\tau}_{\delta\beta}\,.$$

(12)

Hence, the $n$’th order metric perturbation encodes $\accentset{n}{g}_{\alpha\beta}$ the symmetric part $2\accentset{n}{\tau}_{(\alpha\beta)}$ 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

$$\accentset{n}{\tau}_{\alpha\beta}+\accentset{n}{\tau}_{\beta\alpha}=\accentset{n}{g}_{\alpha\beta}-\eta^{\gamma\delta}\sum_{k=1}^{n-1}\accentset{k}{\tau}_{\gamma\alpha}\accentset{n-k}{\tau}_{\delta\beta}\,.$$

(13)

For the antisymmetric part, which constitute the aforementioned independent degree of freedom, one may define another, antisymmetric tensor field by

$$\accentset{n}{a}_{\alpha\beta}=2\accentset{n}{\tau}_{[\alpha\beta]}=\accentset{n}{\tau}_{\alpha\beta}-\accentset{n}{\tau}_{\beta\alpha}\,.$$

(14)

In summary, the perturbative orders $\accentset{n}{\theta}^{\Gamma}{}_{\alpha}$ of the tetrad for $n\geq 1$ are then defined as

$$\accentset{n}{\theta}^{\Gamma}{}_{\alpha}=\frac{1}{2}\eta^{\beta\gamma}\Delta^{\Gamma}{}_{\beta}\left(\accentset{n}{g}_{\gamma\alpha}+\accentset{n}{a}_{\gamma\alpha}-\eta_{\Delta\Theta}\sum_{k=1}^{n-1}\accentset{k}{\theta}^{\Delta}{}_{\gamma}\accentset{n-k}{\theta}^{\Theta}{}_{\alpha}\right)$$

(15)

in terms of $\accentset{n}{g}_{\alpha\beta}$ and $\accentset{n}{a}_{\alpha\beta}$. 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

$$\accentset{2}{a}_{ab}\,,\quad\accentset{3}{a}_{0a}\,,\quad\accentset{4}{a}_{ab}\,.$$

(16)

Note finally that the tetrad $\theta^{\Gamma}{}_{\alpha}$ comes also with an inverse, the frame field $e_{\Gamma}{}^{\alpha}$. It follows from its relation with the coframe field that its background is the diagonal element $\accentset{0}{e}_{\Gamma}{}^{\alpha}=\Delta_{\Gamma}{}^{\alpha}$, while higher perturbation orders are recursively defined as

$$\accentset{n}{e}_{\Gamma}{}^{\alpha}=-\Delta_{\Gamma}{}^{\beta}\sum_{k=0}^{n-1}\accentset{k}{e}_{\Delta}{}^{\alpha}\accentset{n-k}{\theta}^{\Delta}{}_{\beta}\,,$$

(17)

where the tetrad perturbations on the right hand side are further expanded using the rule (15).

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 Aldrovandi and Pereira (2013). In their covariant formulation Krssak et al. (2019), another fundamental field variable is used, which is the spin connection $\accentset{\bullet}{\omega}^{\Gamma}{}_{\Delta\alpha}$. Together with the tetrad and its inverse, it defines another connection, whose coefficients are given by²²2Here we follow the convention used by xTensor, in which the first lower index $\alpha$ on the connection coefficients $\accentset{\bullet}{\Gamma}^{\gamma}{}_{\alpha\beta}$ is the “derivative” index.

$$\accentset{\bullet}{\Gamma}^{\gamma}{}_{\alpha\beta}=e_{\Gamma}{}^{\gamma}\left(\partial_{\alpha}\theta^{\Gamma}{}_{\beta}+\accentset{\bullet}{\omega}^{\Gamma}{}_{\Delta\alpha}\theta^{\Delta}{}_{\beta}\right)\,.$$

(18)

The spin connection is further restricted to be flat and antisymmetric, so that also the affine connection has vanishing curvature and is metric-compatible,

$$\partial_{\alpha}\accentset{\bullet}{\Gamma}^{\gamma}{}_{\beta\delta}-\partial_{\beta}\accentset{\bullet}{\Gamma}^{\gamma}{}_{\alpha\delta}+\accentset{\bullet}{\Gamma}^{\gamma}{}_{\alpha\lambda}\accentset{\bullet}{\Gamma}^{\lambda}{}_{\beta\delta}-\accentset{\bullet}{\Gamma}^{\gamma}{}_{\beta\lambda}\accentset{\bullet}{\Gamma}^{\lambda}{}_{\alpha\delta}=0\,,\quad\accentset{\bullet}{\nabla}_{\gamma}g_{\alpha\beta}=0\,,$$

(19)

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 Golovnev et al. (2017). Choosing this gauge, $\accentset{\bullet}{\omega}^{\Gamma}{}_{\Delta\alpha}\equiv 0$, the connecting coefficients (18) are expressed in terms of the tetrad and its inverse alone Ualikhanova and Hohmann (2019). 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 $\accentset{n}{g}_{\alpha\beta}$ and $\accentset{n}{a}_{\alpha\beta}$.

Finally, the third class of gravity theories discussed here has become known as symmetric teleparallel gravity theories Nester and Yo (1999). In these theories yet another connection is employed, whose coefficients we denote by $\accentset{\times}{\Gamma}^{\gamma}{}_{\alpha\beta}$, and which is assumed to have vanishing curvature and torsion,

$$\partial_{\alpha}\accentset{\times}{\Gamma}^{\gamma}{}_{\beta\delta}-\partial_{\beta}\accentset{\times}{\Gamma}^{\gamma}{}_{\alpha\delta}+\accentset{\times}{\Gamma}^{\gamma}{}_{\alpha\lambda}\accentset{\times}{\Gamma}^{\lambda}{}_{\beta\delta}-\accentset{\times}{\Gamma}^{\gamma}{}_{\beta\lambda}\accentset{\times}{\Gamma}^{\lambda}{}_{\alpha\delta}=0\,,\quad\accentset{\times}{\Gamma}^{\gamma}{}_{\alpha\beta}-\accentset{\times}{\Gamma}^{\gamma}{}_{\beta\alpha}=0\,.$$

(20)

It follows from these properties that there exists a local coordinate system $(x^{\prime\alpha})$, called the coincident gauge, such that its connection coefficients $\accentset{\times}{\Gamma}^{\prime\gamma}{}_{\alpha\beta}$ in this coordinate system vanish Jiménez et al. (2018). Hence, in the standard PPN coordinate system $(x^{\alpha})$, the connection coefficients are given by

$$\accentset{\times}{\Gamma}^{\gamma}{}_{\alpha\beta}=\frac{\partial x^{\gamma}}{\partial x^{\prime\delta}}\frac{\partial x^{\prime\delta}}{\partial x^{\alpha}\partial x^{\beta}}\,.$$

(21)

The post-Newtonian expansion of these connection coefficients is derived from the assumption that the two coordinate systems $(x^{\alpha})$ and $(x^{\prime\alpha})$ are related by an infinitesimal coordinate transformation, induced by a vector field $\xi^{\alpha}$, so that up to the quadratic order one can write

$$x^{\prime\alpha}=x^{\alpha}+\xi^{\alpha}+\frac{1}{2}\xi^{\beta}\partial_{\beta}\xi^{\alpha}+\ldots$$

(22)

From this assumption follows that the connection coefficients take the form

$$\accentset{\times}{\Gamma}^{\gamma}{}_{\alpha\beta}=\partial_{\alpha}\partial_{\beta}\xi^{\gamma}+\frac{1}{2}\left(\xi^{\delta}\partial_{\alpha}\partial_{\beta}\partial_{\delta}\xi^{\gamma}+2\partial_{(\alpha}\xi^{\delta}\partial_{\beta)}\partial_{\delta}\xi^{\gamma}-\partial_{\alpha}\partial_{\beta}\xi^{\delta}\partial_{\delta}\xi^{\gamma}\right)+\ldots\,,$$

(23)

also up to the quadratic order in the vector field $\xi^{\alpha}$. 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 Flathmann and Hohmann (2021)

$$\accentset{2}{\xi}^{a}\,,\quad\accentset{3}{\xi}^{0}\,,\quad\accentset{4}{\xi}^{a}\,.$$

(24)

From these components the perturbative expansion of the connection coefficients $\accentset{\times}{\Gamma}^{\gamma}{}_{\alpha\beta}$ is obtained.

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^{\alpha})=(t,x^{a})$, as a product manifold $M_{4}\cong T_{1}\times S_{3}$, where $t$ is the time coordinate on the manifold $T_{1}\cong\mathbb{R}$, while $(x^{a})$ are the spatial coordinates on $S_{3}$. It follows that for every $t\in T_{1}$, one can find a map $i_{t}:S_{3}\to M_{4},(x^{a})\mapsto(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^{\alpha}$ on $M_{4}$ splits in the form

$$X^{\alpha}\partial_{\alpha}=X^{0}\partial_{0}+X^{a}\partial_{a}$$

(25)

into a temporal part $X^{0}\partial_{0}$, which is parallel to the coordinate lines generated by the time coordinate $t$, and a spatial part $X^{a}\partial_{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}\partial_{a})$, both now being tensor fields on $S_{3}$, hence functions of the spatial coordinates $(x^{a})$. The dependence of the components $X^{\alpha}$ 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 $\partial_{0}$ essentially differently from spatial derivatives $\partial_{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_{\alpha\beta}$, the $3+1$ split can be written as

$$A_{\alpha\beta}\mathrm{d}x^{\alpha}\otimes\mathrm{d}x^{\beta}=A_{00}\mathrm{d}t\otimes\mathrm{d}t+A_{0a}\mathrm{d}t\otimes\mathrm{d}x^{a}+A_{a0}\mathrm{d}x^{a}\otimes\mathrm{d}t+A_{ab}\mathrm{d}x^{a}\otimes\mathrm{d}x^{b}\,,$$

(26)

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.1 we explain how to select certain components from the $3+1$ decomposition of a tensor field, while section V.3 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.

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

$$A=\sum_{n=0}^{\infty}\accentset{n}{A}\,,$$

(27)

where each term is of the corresponding velocity order $\accentset{n}{A}\sim\mathcal{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 Brizuela et al. (2009), where the perturbative expansion takes the form

$$A=\sum_{n=0}^{\infty}\frac{\varepsilon^{n}}{n!}\Delta^{n}[A]\,,$$

(28)

and thus has the explicit form of a Taylor expansion in the perturbation parameter $\varepsilon$. We therefore list a few formulas which follow from the perturbative expansion (27). Most notably, for a tensor $A=A_{1}\cdots A_{N}$ given as a product of $N$ tensors $A_{1},\ldots,A_{N}$, one has the $n$’th order term given by

$$\accentset{n}{A}=\sum_{k_{1}+\ldots+k_{N}=n}\prod_{i=1}^{N}\accentset{k_{i}}{A_{i}}\,,$$

(29)

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 $\accentset{0}{A}$ is used,

$$f(A)=\sum_{k=0}^{\infty}\frac{\left(A-\accentset{0}{A}\right)^{k}}{k!}f^{(k)}\left(\accentset{0}{A}\right)\,,$$

(30)

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

$$f(A_{1},\ldots,A_{N})=\sum_{k_{1}=0}^{\infty}\cdots\sum_{k_{N}=0}^{\infty}f^{(k_{1},\ldots,k_{N})}\left(\accentset{0}{A}_{1},\ldots,\accentset{0}{A}_{N}\right)\prod_{i=1}^{N}\frac{\left(A_{i}-\accentset{0}{A}_{i}\right)^{k_{i}}}{k_{i}!}\,.$$

(31)

Finally, time derivatives are weighted with an additional velocity order, and so for $A^{\prime}=\partial_{0}A$ one has

$$\accentset{n}{A}^{\prime}=\partial_{0}\accentset{n-1}{A}\,,$$

(32)

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^{\prime}$.

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.4. A few notes on the implementation of this function are given in appendix A.3.

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

1. 1.

   MfTime is the one-dimensional time manifold $T_{1}$.

2. 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. 3.

   MfSpacetime is the four-dimensional spacetime manifold $M_{4}\cong T_{1}\times 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 $\mathbb{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 $\mathbb{L}$.

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:

1. 1.

   LI[0] 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.

2. 2.

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

3. 3.

   \[ScriptCapitalT] (printed as $\mathpzc{T}$) is the generic index on the Lorentz bundle $\mathbb{L}T_{1}$. It is used by xTensor in the same way as $\mathpzc{t}$ on $\mathbb{T}T_{1}$.

4. 4.

   Lowercase Latin letters $a,\ldots,z$ (character codes 97–122) are used to denote indices on $\mathbb{T}S_{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.

5. 5.

   Uppercase Latin letters $A,\ldots,Z$ (character codes 65–90) are used to denote indices on $\mathbb{L}S_{3}$. As for the lowercase indices listed in the previous item, xPPN defines symbols L3A, …, L3Z for these uppercase indices, and omits the prefix “L3” when printing the indices in output form.

6. 6.

   Lowercase Greek letters $\alpha,\ldots,\omega$ (character codes 945–969, 977, 981, 982, 1008, 1009, 1013) are used to denote indices on $\mathbb{T}M_{4}$. The corresponding symbols defined by xPPN are denoted T4\textalpha, …, T4\textomega.

7. 7.

   Uppercase Greek letters $\mathrm{A},\ldots,\Omega$ (character codes 913–929, 931–937) are used to denote indices on $\mathbb{L}M_{4}$. The corresponding symbols defined by xPPN are denoted L4\textAlpha, …, L4\textOmega.

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 $\mathbb{T}M_{4}$, use:

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:

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 1. 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^{\alpha}$ one has $A_{\alpha}A^{\alpha}=A^{\alpha}A^{\beta}\eta_{\alpha\beta}$, which can easily be verified by using the xTensor command SeparateMetric:

In order to raise and lower indices with the physical metric $g_{\alpha\beta}$, the metric (as defined in the following section) must be written out explicitly.

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.4.1. A similar role may be taken by the tetrad, as shown in section IV.4.2. Further, three covariant derivatives are defined, which represent the three connections used in different formulations of general relativity Jiménez et al. (2019): the Levi-Civita connection in section IV.4.3, the teleparallel connection in section IV.4.4 and the symmetric teleparallel connection in section IV.4.5.

In the PPN formalism it is assumed that the source of gravity is the energy-momentum tensor $\Theta_{\alpha\beta}$ of a perfect fluid (1), which defined by xPPN is a tensor on $M_{4}$ denoted by EnergyMomentum[-T4\textalpha, -T4\textbeta]. Note that it is defined with lower indices. In order to raise the indices, the physical metric $g_{\alpha\beta}$ must be used explicitly; implicitly raising the indices by writing EnergyMomentum[T4\textalpha, T4\textbeta] would use the background metric $\eta_{\alpha\beta}$, and hence not give the correct result. Also note that we use the symbol $\Theta$ 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

$$\bar{\Theta}_{\alpha\beta}=\Theta_{\alpha\beta}-\frac{1}{2}g_{\alpha\beta}g^{\gamma\delta}\Theta_{\gamma\delta}\,,$$

(37)

denoted by TREnergyMomentum[-T4\textalpha, -T4\textbeta], 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:

1. 1.

   Density[] is the rest mass density $\rho\sim\mathcal{O}(2)$.

2. 2.

   Pressure[] is the pressure $p\sim\mathcal{O}(4)$.

3. 3.

   InternalEnergy[] is the specific internal energy $\Pi\sim\mathcal{O}(2)$.

4. 4.

   Velocity[T3a] is the velocity $v^{a}\sim\mathcal{O}(1)$.

When the energy-momentum tensor $\Theta_{\alpha\beta}$ is expanded in velocity orders as shown in section V.4, the relations

$$\accentset{2}{\Theta}_{00}=\rho\,,\quad\accentset{4}{\Theta}_{00}=\rho\left(\Pi+|v|^{2}-\accentset{2}{g}_{00}\right)\,,\quad\accentset{3}{\Theta}_{0a}=-\rho v_{a}\,,\quad\accentset{4}{\Theta}_{ab}=\rho v_{a}v_{b}+p\delta_{ab}\,,$$

(38)

which follow directly from the expansion (6), are applied, while all other components vanish.

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 $\mathbb{T}S_{3}$. In particular, the following PPN potentials are defined:

1. 1.

   PotentialChi[] is the post-Newtonian superpotential

   $$\chi(t,\vec{x})=-\int\mathrm{d}^{3}x^{\prime}\rho(t,\vec{x}^{\prime})|\vec{x}-\vec{x}^{\prime}|\,.$$

   (39)

2. 2.

   PotentialU[] is the Newtonian potential

   $$U(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\frac{\rho(t,\vec{x}^{\prime})}{|\vec{x}-\vec{x}^{\prime}|}\,.$$

   (40)

3. 3.

   PotentialUU[-T3a, -T3b] is the anisotropic potential

   $$U_{ab}(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\frac{\rho(t,\vec{x}^{\prime})}{|\vec{x}-\vec{x}^{\prime}|^{3}}(x_{a}-x_{a}^{\prime})(x_{b}-x_{b}^{\prime})=\chi_{,ab}-\frac{1}{2}\triangle\chi\delta_{ab}\,.$$

   (41)

4. 4.

   PotentialV[-T3a] is the isotropic vector potential

   $$V_{a}(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\frac{\rho(t,\vec{x}^{\prime})v_{a}(t,\vec{x}^{\prime})}{|\vec{x}-\vec{x}^{\prime}|}\,.$$

   (42)

5. 5.

   PotentialW[-T3a] is the anisotropic vector potential

   $$W_{a}(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\frac{\rho(t,\vec{x}^{\prime})v_{b}(t,\vec{x}^{\prime})(x_{a}-x_{a}^{\prime})(x_{b}-x_{b}^{\prime})}{|\vec{x}-\vec{x}^{\prime}|^{3}}\,.$$

   (43)

6. 6.

   PotentialPhi1[] is the kinetic energy potential

   $$\Phi_{1}(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\frac{\rho(t,\vec{x}^{\prime})v(t,\vec{x}^{\prime})^{2}}{|\vec{x}-\vec{x}^{\prime}|}\,.$$

   (44)

7. 7.

   PotentialPhi2[] is the gravitational self-energy potential

   $$\Phi_{2}(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\frac{\rho(t,\vec{x}^{\prime})U(t,\vec{x}^{\prime})}{|\vec{x}-\vec{x}^{\prime}|}\,.$$

   (45)

8. 8.

   PotentialPhi3[] is the internal energy potential

   $$\Phi_{3}(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\frac{\rho(t,\vec{x}^{\prime})\Pi(t,\vec{x}^{\prime})}{|\vec{x}-\vec{x}^{\prime}|}\,.$$

   (46)

9. 9.

   PotentialPhi4[] is the pressure potential

   $$\Phi_{4}(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\frac{p(t,\vec{x}^{\prime})}{|\vec{x}-\vec{x}^{\prime}|}\,.$$

   (47)

10. 10.

    PotentialA[] is the anisotropic kinetic potential

    $$\mathcal{A}(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\frac{\rho(t,\vec{x}^{\prime})\left[v_{a}(t,\vec{x}^{\prime})(x_{a}-x_{a}^{\prime})\right]^{2}}{|\vec{x}-\vec{x}^{\prime}|^{3}}\,.$$

    (48)

11. 11.

    PotentialB[] is the potential related to change in velocity

    $$\mathcal{B}(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\frac{\rho(t,\vec{x}^{\prime})}{|\vec{x}-\vec{x}^{\prime}|}(x_{a}-x_{a}^{\prime})\frac{\mathrm{d}v_{a}(t,\vec{x}^{\prime})}{\mathrm{d}t}\,.$$

    (49)

12. 12.

    PotentialPhiW[] is the Whitehead potential

    $$\Phi_{W}(t,\vec{x})=\int\mathrm{d}^{3}x^{\prime}\int\mathrm{d}^{3}x^{\prime\prime}\rho(t,\vec{x}^{\prime})\rho(t,\vec{x}^{\prime\prime})\frac{x_{a}-x_{a}^{\prime}}{|\vec{x}-\vec{x}^{\prime}|^{3}}\left(\frac{x_{a}^{\prime}-x_{a}^{\prime\prime}}{|\vec{x}-\vec{x}^{\prime\prime}|}-\frac{x_{a}-x_{a}^{\prime\prime}}{|\vec{x}^{\prime}-\vec{x}^{\prime\prime}|}\right)\,.$$

    (50)

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.6.

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 2. 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 Zaglauer (1990); Helbig (1991). The latter might be implemented by allowing for PPN parameters which depend on TimePar instead of being declared constant Sanghai and Clifton (2017).

As explained in section III.1, 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:

1. 1.

   PPN[$h$][$i$], 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);

2. 2.

   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^{\alpha}$:

Note that the indices do not have to be in the natural position in which they have been for the tensor head $h$. The function PPN yields the same result as if the indices had been specified in their natural position, and then raised or lowed with the background metric BkgMetricS3 on $S_{3}$:

Tensor symmetries are taken into account, and indices are sorted into canonical order if possible:

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. 1.

   ApplyPPNRules[$X$] recursively applies the defined PPN rules to all tensors which appear in the expression $X$ and its subexpressions.

2. 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:

By default, xPPNdoes not associate any rules to the terms in the post-Newtonian expansion of user-defined tensors. Therefore, xPPN provides a number of functions to define and undefine such rules. This can be done with the following functions:

1. 1.

   OrderSet[PPN[$h$, $n$][$i$], $X$] defines or replaces a rule, such that ApplyPPNRules substitutes the $n$’th order term PPN[$h$, $n$][$i$] by $X$, where $X$ can be any tensorial expression which has the same free indices as defined by $i$.

2. 2.

   OrderUnset[PPN[$h$, $n$][$i$]] removes the PPN rule associated with the term PPN[$h$, $n$][$i$] in the post-Newtonian expansion.

3. 3.

   OrderClear[$h$] removes any PPN rules associated to the tensor head $h$.

Their application can be illustrated as follows.