Continuity equations for general matter: applications in numerical relativity

In a curved spacetime, the flux of matter energy–momentum through a surface is in general not balanced by an equal increase in the energy and momentum of the matter contained within the enclosed volume. In the Newtonian picture this mismatch is thought of as resulting from the gravitational force acting on the fluid, but in a general relativistic description no such force exists, and instead the change is a result of the transfer of energy and momentum between the matter and the spacetime curvature. This latter description captures effects that cannot easily be described by Newtonian forces—a good example is a de Sitter expansion in FLRW cosmology, in which the total matter energy in a comoving volume grows over time with the expansion, despite no flux crossing the surface; in a flat space view, energy is created 'out of nothing'.

Describing this transfer of energy between curvature and matter is inherently problematic because of the difficulty of defining an energy for the gravitational field at some specific event; since the field is locally undetectable, its energy must be zero at any point in local inertial coordinates. One can get around this by considering integrated quantities over a finite spatial volume, but in doing so one has to choose a slicing of the four dimensional spacetime, and thus the description becomes inherently gauge dependent. Setting aside the special case of spherical symmetry, the only thing that can be unambiguously defined is the total Arnowitt Deser Misner (ADM) mass and momenta of the spacetime at spatial infinity (assuming asymptotically flatness); the division of this total energy between the spacetime and the matter unavoidably depends on the coordinates chosen to describe the non asymptotically flat region.

While acknowledging these difficulties in interpretation, it is nevertheless useful to define and inspect such quantities during numerical relativity (NR) simulations, or strong gravity simulations more generally, such as those with fixed background metrics. In particular, given the numerical slicing of the spacetime into spatial hypersurfaces, it is possible to reconcile the instantaneous fluxes and energy accumulation of the matter in an arbitrary spatial volume by an application of Gauss's law with an additional 'source' term that quantifies the transfer of energy or momentum from the spacetime curvature. The expressions so obtained are essentially the continuity equation and Euler equations for the conservation of the energy–momentum of the matter, with the effect of curvature treated as a source term.

The calculation of this additional source of energy–momentum from curvature is useful in numerical simulations of strong gravity environments for a number of practical reasons. Primarily, it provides a very useful check on the consistency and correctness of the calculation of the flux and energy accumulation—where the spacetime is not well resolved the quantities will not match up, and their disagreement is a useful measure of the non-convergence of the simulation as a whole. This is particularly important in simulations performed in the decoupling limit where matter is evolved on a fixed or dynamical metric for which the matter backreaction is neglected, because here one cannot check the satisfaction of the Hamiltonian and momentum constraints. Secondly, the source term is a volume integral and in situations where the system reaches a steady state, it tends to converge faster to its final value, and be less prone to boundary issues than the surface fluxes, making it easier to measure in time limited simulations (we will give an example of this in the context of dynamical friction simulations). Finally, it provides an intuitive quantification (albeit a gauge dependent one) of how the curvature and matter are interacting at different events in the spacetime, which can aid one's physical understanding of the evolution in complex dynamical cases.

The purpose of this note is to provide the relevant expressions for these sources for general matter in a format suitable for implementation in NR codes, and discuss applications and potential issues.

The expressions presented here are essentially the same as the energy–momentum conservation expressions that form the basis of the formulation of conservative hydrodynamical variables in NR codes (see any of the standard NR texts, e.g. section 7.3 of Alcubierre [1], section 5.2 of Baumgarte and Shapiro [4] and section 4.3 of Shibata [13]), where satisfying the relevant balance laws is built into the formulation. However, the expressions here are for a general source, and so, for example, they can be applied to real scalar fields and other types of matter with no natural flux-conservative form. In such cases it is relatively uncommon to see the relevant balance law being checked, except in cases with fixed background spacetimes containing Killing vectors (see e.g. [3, 7]), where the sources are identically zero.

The fundamental ideas around conservation in a curved spacetime date back to the Einstein's gravitational pseudotensor [8], and Landau and Lifshitz's alternative pseudotensor [10]. These ideas are thoroughly explained in chapters 19 and 20 of gravitation by Misner Thorne and Wheeler (MTW) [11], and [6] provides an interesting recent review of their philosophical difficulties.

The note is structured as follows:

• In section 2 we start by reviewing the derivation of conserved quantities in a spacetime with a Killing vector, using the 3 + 1 ADM picture.
• In section 3 we cover the case where the vector direction defining the current is no longer a symmetry of the spacetime and give expressions for the additional source terms obtained for time-like (energy) and spatial (linear momentum) directions in terms of the ADM quantities.
• In section 4 we discuss gauge dependence and relate the source to observed quantities. We emphasise in particular that the continuity equations derived focusses on the stress–energy of the matter and not of the spacetime so additional ingredients are required to infer, for example, the ADM mass or momenta of the system as a whole.
• In the final section, section 5 we will discuss three simple illustrative examples that give some additional physical insight to the role of the source: dynamical friction on a black hole (BH), real scalar field solitons ('oscillatons'), and FLRW cosmology.

A more detailed discussion of the application of Gauss's Law to the 4D ADM spacetime is given in the appendix.

Our notations follow the typical conventions in NR, in particular, the units are geometric, G = c = 1 and we use the ADM decomposition [2] ¹ , in which the metric g_μν is expressed in adapted coordinates in terms of a spatial metric γ_ij , a lapse α and a shift vector βⁱ as

$$\begin{equation}\mathrm{d}{s}^{2}=-({\alpha }^{2}-{\beta }^{i}{\beta }_{i})\mathrm{d}{t}^{2}+2{\beta }_{i}\enspace \mathrm{d}{x}^{i}\enspace \mathrm{d}t+{\gamma }_{ij}\enspace \mathrm{d}{x}^{i}\enspace \mathrm{d}{x}^{j},\end{equation} \tag{ 1 }$$

and the extrinsic curvature K_ij is related to the Lie derivative along the normal direction to the spatial slice n_μ = (−α, 0, 0, 0) as

$$\begin{equation}{K}_{ij}=-\frac{1}{2}{\mathcal{L}}_{\hat{n}}{\gamma }_{ij}=\frac{1}{2\alpha }(-{\partial }_{t}{\gamma }_{ij}+{D}_{i}{\beta }_{j}+{D}_{j}{\beta }_{i}).\end{equation} \tag{ 2 }$$

We express the stress energy tensor T_μν in terms of the quantities measured by the normal observers, energy density ρ, momentum density S_i and stress energy density S_ij , which are those typically calculated in NR codes

$$\begin{equation}\rho ={n}^{\mu }{n}^{\nu }{T}_{\mu \nu }\quad {S}_{i}=-{n}^{\mu }{T}_{\mu i}\quad {S}_{ij}={T}_{ij}.\end{equation} \tag{ 3 }$$

From these the full stress energy tensor can be reconstructed as

$$\begin{equation}{T}^{\mu \nu }=\rho {n}^{\mu }{n}^{\nu }+{S}^{\mu }{n}^{\nu }+{S}^{\nu }{n}^{\mu }+{S}^{\mu \nu },\end{equation} \tag{ 4 }$$

noting that in the adapted coordinate basis, S^t = S^μt = 0 (while in general S_t ≠ 0 and S_μt ≠ 0).

We define a current J^μ by projecting the stress energy tensor T_μν in the direction ζ^μ

$$\begin{equation}{J}^{\mu }={\zeta }^{\nu }{T}_{\nu }^{\mu }.\end{equation} \tag{ 5 }$$

If ζ^ν is a Killing vector, ${\mathcal{L}}_{ \overrightarrow {\zeta }}{g}_{\mu \nu }=0$, then ∇_μ J^μ = 0 and the current is conserved throughout the spacetime in the usual 'flat space' sense that the energy accumulated in some volume is balanced by the flux through its surface. Explicitly we have

$$\begin{equation}\int {\mathrm{d}}^{4}x\enspace \sqrt{-g}{\nabla }_{\mu }{J}^{\mu }=0,\end{equation} \tag{ 6 }$$

from which we can convert the 4D volume integral to a 3D surface integral according to Gauss Law. The setup as applied to the ADM slicing is discussed in more detail in the appendix, but in summary one finds that on any 3D spatial hypersurface Σ of constant time coordinate t

$$\begin{equation}{\partial }_{t}{\int }_{{\Sigma}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }Q=-{\int }_{\partial {\Sigma}}{\mathrm{d}}^{2}x\enspace \sqrt{\sigma }F.\end{equation} \tag{ 7 }$$

Here we have identified a charge Q

$$\begin{equation}Q=-{n}_{\mu }{J}^{\mu }=-{n}_{\mu }{\zeta }^{\nu }{T}_{\nu }^{\mu }\end{equation} \tag{ 8 }$$

and a flux F out of the spatial boundary ∂Σ

$$\begin{equation}F=\alpha {N}_{i}{J}^{i}=\alpha {N}_{i}{T}_{\nu }^{i}{\zeta }^{\nu }.\end{equation} \tag{ 9 }$$

For an isolated system with zero net flux through the surface, the total amount of the charge Q is conserved. Here γ is the determinant of the spatial metric, Nⁱ is the outward normal direction to the 2D bounding surface ∂Σ on which the induced metric has determinant σ. Nⁱ is normalised such that γ_ij Nⁱ N^j = 1, and n^μ is the normal to Σ as defined above for which g_μν n^μ n^ν = −1 ² .

Let us first consider energy conservation, and choose to define the energy-current as the flux in the time-like direction defined by ζ^μ = t^μ = (1, 0, 0, 0), which is a Killing vector when ∂_t g_μν = 0, for example in the Kerr metric. Note that we could instead use the time-like direction that defines the surfaces on which we calculate the total conserved charge at some time t, that is n^μ —that of the normal observers. In numerical simulations, n^μ and t^μ are usually distinct for singularity avoidance (a non-constant lapse) and to mitigate the resultant slice-stretching (a non-zero shift). In a gauge with a non zero shift (e.g. the puncture gauge that is commonly used for BHs spacetimes) the normal direction does not correspond to a Killing direction once a stationary state is reached, whereas t^μ does. However, in a dynamical spacetime this is an arbitrary choice, and we will see later that the use of either ζ^μ = t^μ = (1, 0, 0, 0) or ζ_μ = n_μ /α = (1, 0, 0, 0) to define the current is related to using either the Einstein or Landau Lifshitz pseudotensors when defining the conserved gravitational stress energy tensor. For now, taking the case ζ^μ = t^μ = (1, 0, 0, 0), we have

$$\begin{equation}Q={Q}_{t}=-{n}_{\mu }{J}^{\mu }=-{n}_{\mu }{\zeta }^{\nu }{T}_{\nu }^{\mu }=\alpha {T}_{0}^{0}=-\alpha \rho +{\beta }_{k}{S}^{k}\end{equation} \tag{ 10 }$$

and

$$\begin{equation}F={F}_{t}=\alpha {N}_{i}{J}^{i}=\alpha {N}_{i}{T}_{0}^{i}={N}_{i}\left({\beta }^{i}(\alpha \rho -{\beta }^{j}{S}_{j})+\alpha ({\beta }^{k}{S}_{k}^{i}-\alpha {S}^{i})\right) .\end{equation} \tag{ 11 }$$

Using these expressions, we should be able to equate the time derivative of the volume integral of Q_t with (minus) the surface integral of the flux F_t .

One can also consider the case of a spatial Killing vector ${\zeta }^{\mu }={\delta }_{i}^{\mu }$, representing the conservation of momentum (where for example i = x for linear momentum in the x direction). Now we have momentum charge

$$\begin{equation}Q={Q}_{i}=\alpha {T}_{i}^{0}={S}_{i},\end{equation} \tag{ 12 }$$

and the flux is the i-stress in the direction normal to the surface ∂Σ (aka the i-momentum flux through the surface)

$$\begin{equation}F={F}_{i}=\alpha {N}_{j}{T}_{i}^{j}={N}_{j}(\alpha {S}_{i}^{j}-{\beta }^{j}{S}_{i}).\end{equation} \tag{ 13 }$$

Such expressions are commonly used as diagnostics in simulations of matter on a fixed Kerr or Schwarzschild BH background, e.g. in superradiance calculations [7], or accretion onto BHs [3]. Note that such calculations are always approximations—in a spacetime with a truly time independent metric, the matter must have also reached a stationary state, and so evolving matter on a fixed background is a fundamentally inconsistent thing to do (but often a very good approximation for small matter densities, as discussed in section 5).

Now let us consider the case where the current is defined using a direction ${~{\zeta }}^{\mu }$ that is not a Killing vector.

Using the same procedure as in the previous section, we start with the 4D integral of the divergence, which is no longer equal to zero but instead

$$\begin{equation}\int {\mathrm{d}}^{4}\enspace x\sqrt{-g}{\nabla }_{\mu }{J}^{\mu }=\int {\mathrm{d}}^{4}x\enspace \sqrt{-g}{T}_{\nu }^{\mu }{\nabla }_{\mu }{~{\zeta }}^{\nu },\end{equation} \tag{ 14 }$$

from which we obtain

$$\begin{equation}{\partial }_{t}\left({\int }_{{\Sigma}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }Q\right)=-{\int }_{\partial {\Sigma}}{\mathrm{d}}^{2}x\enspace \sqrt{\sigma }F+{\int }_{{\Sigma}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }\mathcal{S}.\end{equation} \tag{ 15 }$$

Here in addition to the charge and the flux, which remain the same as in the previous section, we now have an additional 'source' $\mathcal{S}$

$$\begin{equation}\mathcal{S}=\alpha {T}_{\nu }^{\mu }{\nabla }_{\mu }{~{\zeta }}^{\nu }.\end{equation} \tag{ 16 }$$

Although the source is a scalar quantity, it relates to a particular set of observers for whom Q is the measured charge (in the case in which Q is the energy, for example, ${~{\zeta }}^{\mu }$ is the direction they call time). In most numerical situations one will be interested in quantities measured by coordinate observers, such that ${~{\zeta }}^{\mu }$ is one of the coordinate basis vectors (with components like (1,0,0,0) for time), then ${\partial }_{\mu }{~{\zeta }}^{\nu }=0$ and we can write

$$\begin{equation}\mathcal{S}=\alpha {T}_{\nu }^{\mu }{}^{(4)}{\Gamma}_{\mu \rho }^{\nu }{~{\zeta }}^{\rho }.\end{equation} \tag{ 17 }$$

This expression is calculated in a specific frame, and is manifestly coordinate dependent as one can see from the presence of the Christoffel symbol. While equation (17) covers the cases below, equation (16) is more general and would be required for more complicated scenarios (such as the treatment of angular momentum in Cartesian coordinates).

Physically, we can see from equation (16) that the source describes how the curvature of the spacetime guides the flow of energy–momentum in the fluid from other directions into (or away from) energy–momentum in the ${~{\zeta }}^{\mu }$ direction, thus breaking the conservation of the charge that we would have in flat space (or one possessing the appropriate Killing vector). This interpretation will be made more explicit in section 4.

In the case where we choose ${~{\zeta }}^{\mu }={\delta }_{t}^{\mu }$ we are considering the flux and accumulation of energy in the spacetime, as measured by the coordinate observers. Here the source is ³

$$\begin{equation}{\mathcal{S}}_{t}=\alpha {T}_{\nu }^{\mu }\enspace {}^{(4)}{\Gamma}_{\mu t}^{\nu }=-\rho {\partial }_{t}\alpha +{S}_{i}{\partial }_{t}{\beta }^{i}+\frac{\alpha }{2}{S}^{ij}{\partial }_{t}{\gamma }_{ij}.\end{equation} \tag{ 18 }$$

In the case of a spatial ${~{\zeta }}^{\mu }={\delta }_{i}^{\mu }$ we are considering the flux and accumulation of i-momentum in the spacetime. Here we obtain

$$\begin{equation}{\mathcal{S}}_{i}=\alpha {T}_{\nu }^{\mu }{}^{\enspace (4)}{\Gamma}_{\mu i}^{\nu }=-\rho {\partial }_{i}\alpha +{S}_{j}{\partial }_{i}{\beta }^{j}+\alpha {S}_{j}^{k}{}^{(3)}{\Gamma}_{ik}^{j}.\end{equation} \tag{ 19 }$$

In both cases one can see that the sources are zero in a spacetime with the appropriate Killing vector.

We should therefore be able to reconcile the flux across some surface with the rate of change of the charge in the enclosed volume, plus the volume integral of the source term, for any type of matter. This is an instantaneous expression obeyed for any volume at any time t and so it will be true in any spacetime, not just ones with a static background or those where the outer surface is taken to be in a region that is asymptotically flat (although in such general cases the components may be harder to interpret physically, as discussed further in section 4). Therefore one can integrate the quantities over time during a simulation to confirm that the above expression remains satisfied. This is usually a very good check of whether the different scales of interest are appropriately resolved.

When there is a singularity in the spacetime, we can no longer perform a volume integral over the whole spatial slice, and we therefore need to replace a volume containing the singularity Σ_i with the flux at some inner surface ∂Σ_i enclosing it, i.e.

$$\begin{align}{\partial }_{t}\left({\int }_{{\Sigma}-{{\Sigma}}_{i}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }Q\right)& =-\left({\int }_{\partial {\Sigma}}{\mathrm{d}}^{2}x\enspace \sqrt{\sigma }F-{\int }_{\partial {{\Sigma}}_{i}}{\mathrm{d}}^{2}x\enspace \sqrt{\sigma }F\right)\\ & \quad +{\int }_{{\Sigma}-{{\Sigma}}_{i}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }\enspace \mathcal{S}.\end{align} \tag{ 20 }$$

Σ_i is often chosen to be at or near the horizon, which allows us to relate the measured flux to energy or momentum accretion onto the BH, but this quantity is highly coordinate dependent, and in practice a slightly larger radius can be used (and may be required for a Schwarzschild type gauge where the horizon is a coordinate singularity).

Having worked out the effect of curvature on the matter content, we now describe how this is related to the total energy and momentum of the spacetime.

Einstein's pseudotensor for the gravitational field energy–momentum ${t}_{\nu }^{\mu }$ is defined as

$$\begin{align}{{t}_{\mu }}^{\nu }& =\frac{1}{16\pi \sqrt{-g}}\left({\partial }_{\mu }({g}^{\alpha \beta }\sqrt{-g})({{\Gamma}}_{\alpha \beta }^{\nu }-{\delta }_{\beta }^{\nu }{{\Gamma}}_{\alpha \sigma }^{\sigma })\right.\\ & \quad \left.-{\delta }_{\mu }^{\nu }{g}^{\alpha \beta }({{\Gamma}}_{\alpha \beta }^{\sigma }{{\Gamma}}_{\sigma \rho }^{\rho }-{{\Gamma}}_{\alpha \sigma }^{\rho }{{\Gamma}}_{\beta \rho }^{\sigma })\sqrt{-g}\right) ,\end{align} \tag{ 21 }$$

which can be shown to satisfy the conservation law

$$\begin{equation}{\partial }_{\mu }(\sqrt{-g}({T}_{\nu }^{\mu }+{t}_{\nu }^{\mu }))=0.\end{equation} \tag{ 22 }$$

The pseudotensor ${t}_{\nu }^{\mu }$ is a second order quantity in the first derivatives of g_μν , and so while in local inertial coordinates it will reduce to zero, its integral over a volume is directly related to the energy momentum of the spacetime curvature. One can show that in an asymptotically flat space the ADM mass and momenta P_μ agree at infinity to the volume integral of the relevant charge, i.e.,

$$\begin{equation}{P}_{\mu }=\underset{r\to \infty }{\mathrm{lim}}\int {\mathrm{d}}^{3}x\enspace \sqrt{-g}({T}_{\mu }^{0}+{t}_{\mu }^{0}),\end{equation} \tag{ 23 }$$

assuming the usual conditions regarding asymptotically Minkowski coordinates for the ADM quantities are met (see Gourgoulhon section 8.2.1 [9]).

We have calculated the matter half of this conservation law, that is, (assuming that ${~{\zeta }}^{\mu }$ is the ν coordinate direction) our source is simply

$$\begin{equation}\sqrt{\gamma }{\mathcal{S}}_{\nu }=\sqrt{-g}{T}_{\rho }^{\mu }{}^{\enspace (4)}{\Gamma}_{\mu \nu }^{\rho }={\partial }_{\mu }(\sqrt{-g}{T}_{\rho }^{\mu }{\delta }_{\nu }^{\rho })={\partial }_{\mu }(\sqrt{-g}{T}_{\nu }^{\mu }).\end{equation} \tag{ 24 }$$

and therefore

$$\begin{equation}\sqrt{\gamma }{\mathcal{S}}_{\nu }=-{\partial }_{\mu }(\sqrt{-g}{t}_{\nu }^{\mu }).\end{equation} \tag{ 25 }$$

Thus what we have called the source ${\mathcal{S}}_{\nu }$ does indeed tell us the transfer of energy and momentum between the source and the curvature within the volume, in the framework of the Einstein pseudotensor.

Note that to infer from this the total energy or momentum of the spacetime we need further information—specifically, we need to know the initial volume integral of the curvature charge and any curvature fluxes (i.e., gravitational waves) that exit the volume under consideration. If this is required, one could add the relevant expressions from equation (21) into the simulation, or, if the outer boundary of the volume lies in an asymptotically flat region, recover them from the gravitational wave flux and deficit in the total ADM mass at this surface.

An alternative is provided by the Landau Lifshitz pseudotensor ${t}_{LL}^{\mu \nu }$

$$\begin{align}{t}_{LL}^{\mu \nu }& =\frac{1}{16\pi }\left((2{{\Gamma}}_{\alpha \beta }^{\sigma }{{\Gamma}}_{\sigma \rho }^{\rho }-{{\Gamma}}_{\alpha \rho }^{\sigma }{{\Gamma}}_{\beta \sigma }^{\rho }-{{\Gamma}}_{\alpha \sigma }^{\sigma }{{\Gamma}}_{\beta \rho }^{\rho })({g}^{\mu \alpha }{g}^{\nu \beta }-{g}^{\mu \nu }{g}^{\alpha \beta })\right.\\ & \quad +{g}^{\mu \alpha }{g}^{\beta \sigma }({{\Gamma}}_{\alpha \rho }^{\nu }{{\Gamma}}_{\beta \sigma }^{\rho }+{{\Gamma}}_{\beta \sigma }^{\nu }{{\Gamma}}_{\alpha \rho }^{\rho }-{{\Gamma}}_{\sigma \rho }^{\nu }{{\Gamma}}_{\alpha \beta }^{\rho }-{{\Gamma}}_{\alpha \beta }^{\nu }{{\Gamma}}_{\sigma \rho }^{\rho })\\ & \quad +{g}^{\mu \alpha }{g}^{\beta \sigma }({{\Gamma}}_{\alpha \rho }^{\nu }{{\Gamma}}_{\beta \sigma }^{\rho }+{{\Gamma}}_{\beta \sigma }^{\nu }{{\Gamma}}_{\alpha \rho }^{\rho }-{{\Gamma}}_{\sigma \rho }^{\nu }{{\Gamma}}_{\alpha \beta }^{\rho }-{{\Gamma}}_{\alpha \beta }^{\nu }{{\Gamma}}_{\sigma \rho }^{\rho })\\ & \quad +{g}^{\nu \alpha }{g}^{\beta \sigma }({{\Gamma}}_{\alpha \rho }^{\mu }{{\Gamma}}_{\beta \sigma }^{\rho }+{{\Gamma}}_{\beta \sigma }^{\mu }{{\Gamma}}_{\alpha \rho }^{\rho }-{{\Gamma}}_{\sigma \rho }^{\mu }{{\Gamma}}_{\alpha \beta }^{\rho }-{{\Gamma}}_{\alpha \beta }^{\mu }{{\Gamma}}_{\sigma \rho }^{\rho })\\ & \quad +{g}^{\nu \alpha }{g}^{\beta \sigma }({{\Gamma}}_{\alpha \rho }^{\mu }{{\Gamma}}_{\beta \sigma }^{\rho }+{{\Gamma}}_{\beta \sigma }^{\mu }{{\Gamma}}_{\alpha \rho }^{\rho }-{{\Gamma}}_{\sigma \rho }^{\mu }{{\Gamma}}_{\alpha \beta }^{\rho }-{{\Gamma}}_{\alpha \beta }^{\mu }{{\Gamma}}_{\sigma \rho }^{\rho })\\ & \left.\quad \left.+{g}^{\alpha \beta }{g}^{\sigma \rho }({{\Gamma}}_{\alpha \sigma }^{\mu }{{\Gamma}}_{\beta \rho }^{\nu }-{{\Gamma}}_{\alpha \beta }^{\mu }{{\Gamma}}_{\sigma \rho }^{\nu })\right)+{g}^{\alpha \beta }{g}^{\sigma \rho }({{\Gamma}}_{\alpha \sigma }^{\mu }{{\Gamma}}_{\beta \rho }^{\nu }-{{\Gamma}}_{\alpha \beta }^{\mu }{{\Gamma}}_{\sigma \rho }^{\nu })\right) .\end{align} \tag{ 26 }$$

(See appendix F of [13] or chapter 20 of MTW [11] for a more thorough description.) This allows one to define a conservation law from the relation

$$\begin{equation}{\partial }_{\mu }(-g({T}^{\mu \nu }+{t}_{LL}^{\mu \nu }))=0.\end{equation} \tag{ 27 }$$

It has the advantage over the Einstein pseudotensor that it is index symmetric, making it more appropriate for defining a global angular momentum of the spacetime. Note that the conservation law contains the determinant of the metric g and not the square root, which means that the resulting volume integrals are integrals of scalar densities of weight 1 rather than scalars. However, in cases where the asymptotic behaviour of the spacetime is flat with determinant 1, the volume integral of the charges again correspond at infinity to the ADM charges ⁴ . A second point to note is that the energy current must now be defined by ${J}^{\mu }=(-g){T}^{\mu \nu }{\delta }_{\nu }^{t}$, that is ${~{\zeta }}_{\mu }=(1,0,0,0)$ rather than ${~{\zeta }}^{\mu }=(1,0,0,0)$. This implies that the conservation of charges in the normal way is associated with the case where the normal direction n^μ is a Killing symmetry rather than the time-like direction as in the Einstein pseudotensor description. As discussed above, in most numerical simulations, the stationary state obtained corresponds to the metric being constant in time, not in the normal direction n^μ . Nevertheless one can construct the equivalent expressions in this formalism and use them in the same way, for example in the case of energy conservation one obtains that

$$\begin{equation}{\partial }_{t}\left({\int }_{{\Sigma}}{\mathrm{d}}^{3}x\enspace \mathcal{Q}\right)=-{\int }_{\partial {\Sigma}}{\mathrm{d}}^{2}x\enspace \mathcal{F}+{\int }_{{\Sigma}}{\mathrm{d}}^{3}x\mathcal{S},\end{equation} \tag{ 28 }$$

with the charge

$$\begin{equation}\mathcal{Q}=(-g){T}^{00}=\gamma \rho \end{equation} \tag{ 29 }$$

and flux (using the flat space normal direction s_i )

$$\begin{equation}\mathcal{F}=(-g){T}^{i0}{s}_{i}=\gamma (\alpha {S}^{i}-{\beta }^{i}\rho ){s}_{i}.\end{equation} \tag{ 30 }$$

The source is now

$$\begin{align}\mathcal{S}& =(-g)\left({T}^{\mu \nu }{}^{(4)}{\Gamma}_{\mu \nu }^{t}-{T}^{\mu t}{}^{(4)}{\Gamma}_{\mu \nu }^{\nu }\right)\\ & =\gamma \left[\rho \left({\partial }_{i}{\beta }^{i}-\alpha K-\frac{{\beta }^{i}}{\alpha }{\beta }^{j}{K}_{ij}\right)+\alpha {K}^{ij}{S}_{ij}+{S}^{i}(\alpha {{\Gamma}}_{ij}^{j}-{\partial }_{i}\alpha )\right].\end{align} \tag{ 31 }$$

Note the additional Christoffel symbol resulting from the scalar density weight. The behaviour of the different formalisms is compared in figure 2 for the case of an oscillaton (as discussed in section 5).

To summarise—assuming that we have an isolated system in an asymptotically flat region, we have proposed splitting the contributions to the ADM mass into a matter and a curvature component, and noted that by tracking the matter contribution we can infer the changes to the curvature contribution. However, as set out in the introduction, such a split is gauge dependent—the total ADM charges are the only unambiguous physical observables. In a more general non asymptotically flat dynamical spacetime, the problem is even worse, since the total itself is dependent on both the gauge and the formalism chosen (i.e., which curvature pseudotensor we define). So is it in fact useful to calculate such a split?

We would argue yes—it is still useful to track the volume integrals of the energy and source terms in order to gain an intuition about the transfer of energy between the spacetime and the matter, although one should take care not to make quantitative conclusions without considering the gauge dependence. In section 5 we will give several simple examples that illustrate how this approach can be informative. If nothing else, the reconciliation of the matter charges is a useful check on the accuracy of the simulation, and often helps to identify issues such as boundary reflections, underresolved horizons and general code bugs.

In this section we discuss several specific simple examples that give an intuition for the continuity equation in general spacetimes, illustrate some of the points made above, and highlight its potential applications in numerical simulations of strong gravity environments.

5.1. Dynamical friction: the limit of small backreaction

Dynamical friction is the gravitational drag force on an object due to the build up of an overdensity of matter behind it. A classic example is a BH moving through a fluid and building up a gravitational 'wake' of particles behind it. In the Newtonian description of such a system, the dynamical friction force Fⁱ is simply the integral of the gravitational force of each fluid element acting on the object, that is

$$\begin{equation}{F}_{i}=\int \frac{GM\rho }{{r}^{2}}{s}_{i}\enspace \mathrm{d}V,\end{equation} \tag{ 32 }$$

with ${s}_{i}=\frac{{x}^{i}}{r}$ defining the projection of the force in the i-direction. In the fully relativistic case the concept of a gravitational force is less straightforward to define—gravity is not a force but a curvature of spacetime—so one must instead define the force as the rate of change of the ADM momentum of the spacetime F_i = ∂_t P_i . This is assumed to be measured by some observer in an asymptotically flat region where their proper time τ = t. (Strictly ADM charges are defined at infinity, where the fluxes should be zero, but the relations below hold approximately at a finite radius. To be clear we will write these approximate ADM charges as ${~{P}}_{\mu }$.) Assuming the fluid has reached a steady state, we can get the force directly by using the surface integral equivalent to equation (32)

$$\begin{equation}{F}_{i}=\int {T}_{i}^{j}\enspace \mathrm{d}{S}_{j}={\partial }_{t}{~{P}}_{i},\end{equation} \tag{ 33 }$$

where the last equality assumes that the surface integral is evaluated in the asymptotically flat region, such that the result will hold also in GR. One issue with this is that what we really want is the force on the BH due to the fluid, that is given that

$$\begin{equation}{~{P}}_{i}={P}_{i}^{\text{curvature}}+{P}_{i}^{\text{matter}},\end{equation} \tag{ 34 }$$

we want the part associated with ${\partial }_{t}{P}_{i}^{\text{curvature}}$ only. If the matter has reached a steady state ${\partial }_{t}{P}_{i}^{\text{matter}}=0$ then equation (33) gives the correct result. However, in a dynamical simulation often the matter takes time to reach a steady state around the BH, and it is useful to track the quantity specifically associated with the curvature in the period leading up to the stationary state being reached. One usually finds that where a steady state is eventually reached, this term settles to a constant value well before the surface integral of the flux does so, which can reduce the time to solution for states that take a long time to settle.

Consider a case where the stress energy of the fluid is small T_μν ∼ O(), then to zeroth order in the spacetime is described by a vacuum solution ${g}_{\mu \nu }^{(0)}$. At first order in the propagation of the stress energy tensor is described by its evolution on the vacuum background

$$\begin{equation}{\nabla }_{\rho }({g}_{\mu \nu }^{(0)}){T}^{\rho \sigma }=0.\end{equation} \tag{ 35 }$$

There will be a resulting first order correction to the metric δg_μν ∼ O(), which will have the effect of changing the energy–momentum of the spacetime curvature ${P}_{\mu }^{\text{curvature}}$—it is this that results in the gravitational 'force' on the BH. This correction will affect the propagation of the stress energy tensor, but only at higher orders in , so we can neglect this backreaction effect, along with the gravitational wave flux ${t}_{\mu \nu }^{\text{GW}}\sim O({{\epsilon}}^{2})$. Taking care to include the effect of the flux of the matter across the BH horizon ∂Σ_BH (which we take to contribute to ${P}_{i}^{\text{curvature}}$ for the BH), we find that to first order in

$$\begin{equation}{F}_{i}\equiv {\partial }_{t}{P}_{i}^{\text{curvature}}=-{\int }_{{\Sigma}-{{\Sigma}}_{\text{BH}}}\mathrm{d}V\enspace {\mathcal{S}}_{i}-{\int }_{\partial {{\Sigma}}_{\text{BH}}}\mathrm{d}S\enspace \alpha {N}_{j}{T}_{i}^{j}.\end{equation} \tag{ 36 }$$

A few comments on this result:

• The force is composed of two contributions, the first being the dynamical friction effect of the cloud of matter behind the BH, and the second relating to accretion of the fluid's momentum onto the BH.
• The split between these components is highly gauge dependent—in a Schwarzschild type metric clearly any flux measured at the horizon will be exactly zero, unlike in a horizon penetrating gauge.
• As discussed above, the source ${\mathcal{S}}_{i}$ quantifies the transfer of momentum between matter and curvature, which is the relativistic analogy of the Newtonian gravitational force being equal and opposite between two bodies.
• One can easily show that equation (19) for ${\mathcal{S}}_{i}$ reduces to equation (33) in the Newtonian limit.

The agreement in the quantities extracted and their evolution over time is illustrated in figure 1, coming from simulations carried out in [14].

Zoom In Zoom Out Reset image size

5.2. Oscillatons: quasi stable dynamical spacetime

Quasi-stable localised configurations of massive, real scalar fields, called oscillatons [12] balance the gravitational tendency of the field to collapse under gravity with its tendency to disperse due to gradient pressure. These are long lived localised solitons—a roughly Gaussian-looking profile in the field that oscillates in time. Unlike the complex scalar case, the real scalar cannot support a time invariant stress–energy tensor and so both the metric and the field undergo oscillations in time. On each oscillation energy is transferred between the field and the spacetime curvature, and back again.

Due to the spherical symmetry of these solutions, no gravitational radiation is emitted, but scalar radiation is emitted on extremely long timescales, resulting in their eventual decay. On the timescales of an NR simulation, such a flux is negligible and so the initial ADM mass of the oscillaton is roughly conserved. The transfer of energy between the field and the metric can be tracked by considering equation (18) for the matter energy source, and the result is plotted in figure 2 for the formalism based on both the Einstein and the Landau Lifshitz pseudotensors.

Zoom In Zoom Out Reset image size

One can see for both cases that there is an initial transient in the source, which is unphysical, and relates solely to the solution settling into its steady state in the dynamical puncture gauge used in the simulations, which is different from the areal polar gauge in which the initial data are derived. Eventually the source and mass of the matter settle into a steady oscillating state. Naively we might say that the total energy in the matter has decreased—but this cannot have happened since, as seen in the plots, there has been no flux out of the surface of the enclosing volume. In fact, since we have a new gauge and thus a new slicing of the 4D spacetime, what has changed is simply the split of what we denote to be matter energy and curvature energy. This illustrates clearly the gauge dependence of the split and the danger of interpreting such results too literally. However, once a quasi stable state is reached in the chosen gauge, the split is reasonably well defined and does tell us something physical about the oscillations of the energy between curvature and matter.

5.3. FLRW spacetime expansion: no asymptotically flat space

Another interesting illustration of the ideas above comes from their application to a spatially flat FLRW spacetime, in which there is no asymptotically flat space and no time-like Killing vector.

One can show that due to the spatial flatness and isotropy, the total energy charge $E=-g({T}_{0}^{0}+{t}_{0}^{0})$ in any spatial volume is zero, and any volume traced out by the comoving observers does not have a flux though its boundary. However, while the total energy is fixed (at zero), the energy content in the matter and the spacetime curvature can change over time, leading to the well known dependence of the behaviour of cosmological spacetimes on matter content.

Working in the comoving gauge with time defined by the proper time of the comoving observers,

$$\begin{equation}\mathrm{d}{s}^{2}=-\mathrm{d}{t}^{2}+{a}^{2}(t){\delta }_{ij}\enspace \mathrm{d}{x}^{i}\enspace \mathrm{d}{x}^{j},\end{equation} \tag{ 37 }$$

our continuity equation for the matter energy in a unit coordinate spatial volume (with the source equation (18)) tells us that

$$\begin{equation}{\partial }_{t}({a}^{3}\rho )={a}^{3}\left(\frac{1}{2}{S}^{ij}{\partial }_{t}{\gamma }_{ij}\right)=-3{a}^{2}\dot {a}P.\end{equation} \tag{ 38 }$$

Here P and ρ are the pressure and energy density of the perfect fluid matter. This is the usual continuity equation of cosmological expansion, but written in a way that highlights the fact that what breaks 'normal' energy conservation is the presence of pressure in the fluid—dust-like fluids with P = 0 maintain the same total mass-energy even when the volume expands (the particles just get more spread out) ⁵ . A positive pressure will tend to transfer energy from the matter to the spacetime curvature (e.g. in the case of radiation), whereas a negative pressure (e.g. a cosmological constant) will do the reverse.

We note that the requirement that E = 0 on each spatial slice translates to the first Friedmann equation

$$\begin{equation}\int {\mathrm{d}}^{3}x\sqrt{-g}({T}_{0}^{0}+{t}_{0}^{0})=0\enspace \Longrightarrow\enspace {H}^{2}=\frac{8\pi \rho }{3},\end{equation} \tag{ 39 }$$

which highlights the additional information required for full knowledge of the total energy–momentum of the spacetime (the requirement that ∂_t E = 0 gives the second Friedmann equation, but the same information can be obtained from a combination of equations (39) and (38)).

In this note we have provided an overview of the source terms that break the conservation of energy and momentum for general matter evolution in curved spacetimes. We have provided the relevant expressions for these sources in the ADM picture, and discussed applications and potential issues for several examples in the context of strong gravity simulations. We hope that this provides a useful reference on the topic, in particular for new students of NR.

KC thanks Jean Alexandre and her GRChombo collaborators (https://grchombo.org/#people), in particular Jamie Bamber, Robin Croft, Thomas Helfer, Eugene Lim, Miren Radia and Dina Traykova, for useful discussions, and acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No. 693024). The simulations presented in this paper used DiRAC resources under the Projects ACSP218 and ACTP238. This work was performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (https://dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital Grants ST/P002307/1 and ST/R002452/1 and STFC operations Grant ST/R00689X/1. In addition used the DiRAC at Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (https://dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital Grants ST/P002293/1 and ST/R002371/1, Durham University and STFC operations Grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

No new data were created or analysed in this study.

The simplest way of applying Gauss's Law to a 4D spacetime is to start with a conservation law manipulated into a form in which it is expressed with partial derivatives only

$$\begin{equation}{\partial }_{\mu }{~{J}}^{\mu }=0.\end{equation} \tag{ 40 }$$

For example, this can be achieved by defining ${~{J}}^{\mu }=\sqrt{-g}{J}^{\mu }=\sqrt{-g}{T}_{\nu }^{\mu }{\zeta }^{\nu }$ with ζ^ν a Killing vector of the spacetime. Then one can write

$$\begin{equation}{\partial }_{t}{~{J}}^{0}=-{\partial }_{i}{~{J}}^{i}\end{equation} \tag{ 41 }$$

and integrate both sides over a 3D spatial coordinate volume Σ bounded by the 2D surface ∂Σ—in which case one can directly apply Gauss theorem in 3D to find that

$$\begin{equation}{\partial }_{t}{\int }_{{\Sigma}}{~{J}}^{0}\enspace \mathrm{d}V=-{\int }_{\partial {\Sigma}}{~{J}}^{i}\enspace \mathrm{d}{S}_{i}^{f},\end{equation} \tag{ 42 }$$

where $\mathrm{d}{S}_{i}^{f}$ is the coordinate surface integral element (i.e., it does not contain factors of the metric, unlike the proper physical surface element dS_i ), and we assume that the spatial volume Σ is a constant coordinate volume with respect to time. For example in the case of the energy current for the coordinate observers, where ${\zeta }^{\nu }={\delta }_{0}^{\nu }$, we obtain

$$\begin{equation}{\partial }_{t}{\int }_{V}\alpha {T}_{0}^{0}\sqrt{\gamma }\enspace \mathrm{d}V=-{\int }_{\partial V}\alpha {T}_{0}^{i}\sqrt{\sigma }\enspace \mathrm{d}{S}_{i},\end{equation} \tag{ 43 }$$

where we have expanded out the$\sqrt{-g}$ factors in terms of the lapse and spatial metric to make the volume elements for the integrals physically meaningful. From this we can identify directly the conserved charge as $\alpha {T}_{0}^{0}$ and the flux as $\alpha {T}_{0}^{i}$.

A couple of implementation notes:

• Do not forget that J^μ is a four vector with non zero time-like components, and not a spatial three-vector, so one cannot raise and lower the indices on Jⁱ with the spatial metric. In particular, this means that ${J}^{i}\enspace \mathrm{d}{S}_{i}={T}_{0}^{i}\mathrm{d}{S}_{i}\ne {T}_{i0}\enspace \mathrm{d}{S}^{i}={J}_{i}\enspace \mathrm{d}{S}^{i}$.
• In practice it is often useful to use the form with the flat surface integral element as this avoids the need to appropriately normalise dS_i and calculate explicitly the induced metric on the surface σ. One can simply use the fact that $\sqrt{\sigma }{N}_{i}=\sqrt{\gamma }{s}_{i}$ where s_i is a unit vector in the coordinate basis, i.e. δ^ij s_i s_j = 1.

Despite the elegance of this method, it is quite an instructive exercise to apply Gauss theorem directly to the 4D spacetime volume, particularly in the context of the ADM decomposition where we may have a non trivial shift and lapse, and convince ourselves that all the relevant metric factors work out and the four vectors reduce appropriately 3D spatial vectors. Below we derive explicitly the continuity equation using Gauss Law on the 4D volume in the ADM decomposition. The set up is illustrated in figure 3. This is adapted from section 3.5 of [4], but we will make some of the choices regarding volumes more explicit.

Zoom In Zoom Out Reset image size

Our starting point is the 4D volume integral written in covariant form

$$\begin{equation}\int {\mathrm{d}}^{4}x\enspace \sqrt{-g}{\nabla }_{\mu }{J}^{\mu }=0.\end{equation} \tag{ 44 }$$

We assume that we have a conserved current J^μ such that ∇_μ J^μ = 0. Note that ${J}^{\mu }={T}_{\nu }^{\mu }{\zeta }^{\nu }$ is now properly a vector, not a tensor density like ${~{J}}^{\mu }$ was.

In the ADM decomposition there are two natural 4D volumes we could consider, volume 1 is the volume between the slice Σ₁ and Σ₁(t + dt), defined by the congruence of the normal observers from one time slice to the next. Volume 2 is between Σ₁ and Σ₂, which have a constant coordinate volume in the spatial hyperslices, and are thus defined by the congruence of the time like observers. We will consider both in turn and show that the results agree—the considerations in each case are useful reminders of the physical properties of the slicing.

For volume 1:

• The top and bottom spatial hyperslices, Σ₁ and Σ₁(t + dt), are two 3D slices of constant t. They have normals ±n^μ and the related metric determinant is γ.
• The spatial hyperslices are each bounded by a 2D surface ∂Σ on which the induced metric has the determinant σ.
• The sides of the volume trace the normal lines between the two slices. The normal to these sides lie wholly in the spatial surface, n^μ N_μ = 0, thus N^μ = (0, Nⁱ ), with the direction of Nⁱ defined via the 3D covector s_i = ∂_i r (so s_i = (1, 0, 0) in spherical coordinates or xⁱ /r in Cartesian ones) and normalised such that g_μν N^μ N^ν = γ_ij Nⁱ N^j = 1. Thus explicitly,
  $$\begin{equation}{N}^{i}=\frac{{\gamma }^{ij}{s}_{j}}{\sqrt{{\gamma }^{ij}{s}_{i}{s}_{j}}}.\end{equation} \tag{ 45 }$$
  Note that we find from lowering the indices on N^μ that N₀ = βⁱ N_i .
• The separation between the two spatial hypersurfaces in the normal n^μ direction has proper length ds = αdt.

For volume 2:

• The top and bottom slices, Σ₁ and Σ₂, are two (3D) slices of constant t. They have normals ±n^μ and the metric determinant is γ.
• The spatial hyperslices are each bounded by a 2D surface ∂Σ', which is a surface of constant r and metric determinant σ'. (These coincide with ∂Σ and σ on Σ₁ but not in the rest of the volume.)
• The sides of the volume trace the time vector t^μ = (1, 0, 0, 0)—giving a surface of constant radial coordinate r.
• The normal to the sides of the volume is the normal to the surface of constant r, thus it has a direction defined via the 4D covector s_μ = ∂_μ r (so s_μ = (0, 1, 0, 0) in spherical coords) and normalised such that ${g}_{\mu \nu }{~{N}}^{\mu }{~{N}}^{\nu }=1$. One can show that this leads to the relation
  $$\begin{equation}{~{N}}_{i}=\frac{\alpha {N}_{i}}{\sqrt{{\alpha }^{2}+{({\beta }^{i}{N}_{i})}^{2}}}.\end{equation} \tag{ 46 }$$
• The separation between the two spatial hypersurfaces in the time-like (t^μ = (1, 0, 0, 0)) direction has proper length $\mathrm{d}s={({\alpha }^{2}+{({\beta }^{i}{N}_{i})}^{2})}^{1/2}\mathrm{d}t$.

Finally, note that in a 4D Lorentzian spacetime a minus sign enters for space-like boundaries but not for time-like ones in Gauss's Law.

Starting with volume 1, we convert the 4D volume integral to a sum of 3D surface integrals

$$\begin{align}\int {\mathrm{d}}^{3}x({\mathcal{N}}^{\mu }{J}_{\mu })& =-{\int }_{{{\Sigma}}_{1}(t+\mathrm{d}t)}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }({n}_{\mu }{J}^{\mu })+{\int }_{{{\Sigma}}_{1}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }({n}_{\mu }{J}^{\mu })\\ & \quad +{\int }_{t}^{t+\mathrm{d}t}\mathrm{d}~{t}{\int }_{\partial {\Sigma}}{\mathrm{d}}^{2}x\enspace \alpha \sqrt{\sigma }({N}_{\mu }{J}^{\mu })=0.\end{align} \tag{ 47 }$$

We want to consider the limit of this volume when dt is small, but we need to take care about the effect of the changing volume of the slices. In particular:

$$\begin{equation}{\partial }_{t}{\int }_{{{\Sigma}}_{1}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }({n}_{\mu }{J}^{\mu })=\underset{\mathrm{d}t\to 0}{\mathrm{lim}}\enspace \frac{1}{\mathrm{d}t}\left({\int }_{{{\Sigma}}_{2}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }({n}_{\mu }{J}^{\mu })-{\int }_{{{\Sigma}}_{1}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }\enspace ({n}_{\mu }{J}^{\mu })\right)\end{equation} \tag{ 48 }$$

$$\begin{align}& =\underset{\mathrm{d}t\to 0}{\mathrm{lim}}\frac{1}{\mathrm{d}t}\left({\int }_{{{\Sigma}}_{1}(t+\mathrm{d}t)}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }({n}_{\mu }{J}^{\mu })\right.\\ & \quad \left.-{\int }_{{{\Sigma}}_{1}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }({n}_{\mu }{J}^{\mu })\right)+{\int }_{\partial {{\Sigma}}_{1}}{\mathrm{d}}^{2}x\enspace \sqrt{\sigma }({\beta }^{i}{N}_{i})({n}_{\mu }{J}^{\mu }).\end{align} \tag{ 49 }$$

This means that, dividing equation (47) by dt and taking the limit as dt → 0 we obtain

$$\begin{equation}0={\partial }_{t}{\int }_{{{\Sigma}}_{1}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }({n}_{\mu }{J}^{\mu })-{\int }_{\partial {{\Sigma}}_{1}}{\mathrm{d}}^{2}x\enspace \sqrt{\sigma }(\alpha {N}_{\mu }{J}^{\mu }+{\beta }^{i}{N}_{i}{n}_{\mu }{J}^{\mu }).\end{equation} \tag{ 50 }$$

Thus, remembering that N₀ = βⁱ N_i , we have

$$\begin{equation}{\partial }_{t}{\int }_{{{\Sigma}}_{1}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }({n}_{\mu }{J}^{\mu })={\int }_{\partial {{\Sigma}}_{1}}{\mathrm{d}}^{2}x\enspace \sqrt{\sigma }(\alpha {N}_{i}{J}^{i}),\end{equation} \tag{ 51 }$$

from which we can recover equation (43).

Now for fun, let us check that we get the same result using volume 2 instead ⁶ . Now we have:

$$\begin{align}& {\int }_{{{\Sigma}}_{2}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }({n}_{\mu }{J}^{\mu })-{\int }_{{{\Sigma}}_{1}}{\mathrm{d}}^{3}x\sqrt{\gamma }({n}_{\mu }{J}^{\mu })-{\int }_{t}^{t+\mathrm{d}t}\mathrm{d}~{t}\\ & \qquad {\times}{\int }_{\partial {{\Sigma}}^{\prime }}{\mathrm{d}}^{2}x({({\alpha }^{2}+{({\beta }^{i}{N}_{i})}^{2})}^{1/2})\sqrt{{\sigma }^{\prime }}({~{N}}_{\mu }{J}^{\mu })=0.\end{align} \tag{ 52 }$$

Since now the coordinate volume of the spatial slices does not change, the first part is less fiddly. So we get directly that on the initial slice (where σ = σ')

$$\begin{equation}\left.0={\partial }_{t}{\int }_{{{\Sigma}}_{1}}{\mathrm{d}}^{3}x\enspace \sqrt{\gamma }({n}_{\mu }{J}^{\mu })-{\int }_{\partial {{\Sigma}}_{1}}{\mathrm{d}}^{2}x\enspace \sqrt{\sigma }({({\alpha }^{2}+{({\beta }^{i}{N}_{i})}^{2})}^{1/2}){~{N}}_{i}{J}^{i}\right) .\end{equation} \tag{ 53 }$$

Using the relation between the two normals in equation (46) we see that we recover the same result as in equation (51).