Implications of Symmetry and Pressure in Friedmann Cosmology. I. Formalism

Given Assumptions 1 and 2, we work exclusively on ${ \mathcal U }$. In the context of gravity, the action principle demands that

$$\begin{eqnarray}&&\delta \left\{{S}_{G}\left[{g}_{\mu \nu }\right]+{S}_{M}\left[{g}_{\mu \nu },{\boldsymbol{\Psi }}(\eta ,{\boldsymbol{x}})\right]\right\}{| }_{{\boldsymbol{\Psi }}}\equiv 0,\end{eqnarray} \tag{ 16 }$$

where S_G is a gravitational action, S_M is a matter action, and ${\boldsymbol{\Psi }}$ is held fixed for the variation δ. By the definition given in Weinberg (1972, Section 12.2), the symmetric, rank-(2,0), stress tensor ${T}^{\mu \nu }$ enters Einstein's equations through

$$\begin{eqnarray}&&\delta {S}_{M}{| }_{{\boldsymbol{\Psi }}}\equiv \displaystyle \frac{1}{2}{\int }_{{ \mathcal U }}{T}^{\mu \nu }\delta {g}_{\mu \nu }\sqrt{-g}\,{d}^{4}x.\end{eqnarray} \tag{ 17 }$$

Here, g is the metric determinant. Note that we have flipped indices, relative to Weinberg's definition, and thus incur a minus sign through the metric variation

$$\begin{eqnarray}&&\delta \left({{\delta }^{\alpha }}_{\beta }\right)=0=\delta \left({g}^{\alpha \mu }{g}_{\beta \mu }\right).\end{eqnarray} \tag{ 18 }$$

Weinberg's definition encodes the standard assumption that S_M comes from a Lagrange scalar density ${{\mathscr{L}}}_{M}$,

$$\begin{eqnarray}&&\delta {S}_{M}={\int }_{{ \mathcal U }}\delta \left({{\mathscr{L}}}_{M}\sqrt{-g}\right)\,{d}^{4}x\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&={\int }_{{ \mathcal U }}\left[\displaystyle \frac{1}{\sqrt{-g}}\displaystyle \frac{\delta }{\delta {g}_{\mu \nu }}\left({{\mathscr{L}}}_{M}\sqrt{-g}\right)\right]\delta {g}_{\mu \nu }\sqrt{-g}\,{d}^{4}x.\end{eqnarray} \tag{ 20 }$$

Comparison of Equation (17) with (20) reveals that, under the integral, one may identify

$$\begin{eqnarray}&&{T}^{\mu \nu }=\displaystyle \frac{2}{\sqrt{-g}}\displaystyle \frac{\delta }{\delta {g}_{\mu \nu }}\left\{{{\mathscr{L}}}_{M}\left[{\boldsymbol{\Psi }}(\eta ,{\boldsymbol{x}}),{g}_{\mu \nu }\right]\sqrt{-g}\right\}.\end{eqnarray} \tag{ 21 }$$

The stress tensor components are, therefore, a mixture of the metric and the non-gravitational contributions ${\boldsymbol{\Psi }}(\eta ,{\boldsymbol{x}})$. The non-gravitational contributions ${\boldsymbol{\Psi }}(\eta ,{\boldsymbol{x}})$ have no explicit dependence on ${\boldsymbol{g}}$, and thus are unaltered by the model builder's choice of metric model. Of course, the ${\boldsymbol{\Psi }}(\eta ,{\boldsymbol{x}})$ are implicitly constrained by the eventual equations of motion. Any premature constraint, before equations of motion consistent with the action principle are determined, risks introducing inconsistencies into the field equations.

Let R be the Ricci scalar, built from the Levi-Civita connection of any metric. The EH gravitational action is

$$\begin{eqnarray}&&{S}_{G}\equiv \displaystyle \frac{1}{16\pi G}{\int }_{{ \mathcal U }}R\sqrt{-g}\,{d}^{4}x.\end{eqnarray} \tag{ 22 }$$

We use the representation of the metric given in Equation (14). For the purposes of computation, we now parameterize ${ \mathcal U }$. While there is observational evidence that the universe began in a nearly singular hot and dense point, it is not clear how far back in time one can "safely" apply GR. We thus restrict consideration to

$$\begin{eqnarray}&&\eta \in \left[{\eta }_{i},{\eta }_{f}\right].\end{eqnarray} \tag{ 23 }$$

Here ${\eta }_{i}$ is some arbitrary initial time and ${\eta }_{f}\gt {\eta }_{i}$ can either be arbitrary, in the case of matter domination, or any time below the asymptotic value in the case of dark energy domination. Note that because ${ \mathcal U }$ is compact, an action defined over it will never diverge. For the moment, consider a closed 3-ball ${ \mathcal V }\subset {{\mathbb{R}}}^{3}$ so that

$$\begin{eqnarray}&&{ \mathcal U }\equiv { \mathcal V }\times \left[{\eta }_{i},{\eta }_{f}\right].\end{eqnarray} \tag{ 24 }$$

We will often use ${ \mathcal V }$ to represent both this compact submanifold and its volume, which is the integral of the Euclidean volume form over this submanifold.

We now work only through O(1). From Equation (14), to leading order, the approximation is

$$\begin{eqnarray}&&{g}_{\mu \nu }={a}^{2}(\eta )\left[{\eta }_{\mu \nu }+O(\epsilon )\right].\end{eqnarray} \tag{ 25 }$$

This gives, for the metric determinant and inverse metric variation,

$$\begin{eqnarray}&&\sqrt{-g}={a}^{4}+O(\epsilon ),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\delta {g}^{\mu \nu }=-2{\eta }^{\mu \nu }{a}^{-3}\delta a+O(\epsilon ).\end{eqnarray} \tag{ 27 }$$

Following Wald (2010, Equation (D.9), p. 446), we may exploit that the flat RW metric is a conformal rescaling of flat space to write

$$\begin{eqnarray}&&R=-6{\eta }^{\mu \nu }{a}^{-3}{\partial }_{\mu }{\partial }_{\nu }a+O(\epsilon ),\end{eqnarray} \tag{ 28 }$$

where the derivatives are simple partials. Note that we have used Assumption 2 to regard the Ricci scalar of the bare (not multiplied by the scale factor) metric as $O(\epsilon )$.

Substitution into the definition of the gravitational action S_G gives

$$\begin{eqnarray}&&{S}_{G}=-\displaystyle \frac{3}{8\pi G}{\int }_{{ \mathcal U }}{\eta }^{\mu \nu }a({\partial }_{\mu }{\partial }_{\nu }a)\,{d}^{4}x+O(\epsilon ).\end{eqnarray} \tag{ 29 }$$

Performing the variation gives

$$\begin{eqnarray}&&\delta {S}_{G}=-\displaystyle \frac{3}{8\pi G}{\int }_{{ \mathcal U }}\left[\delta a({\partial }^{\mu }{\partial }_{\mu }a)+a({\partial }^{\mu }{\partial }_{\mu }\delta a)\right]{d}^{4}x\end{eqnarray} \tag{ 30 }$$

plus terms $O(\epsilon )$. Integrating by parts twice gives

$$\begin{eqnarray}&&\delta {S}_{G}=-\displaystyle \frac{3}{4\pi G}{\int }_{{ \mathcal U }}\delta a{\partial }^{\mu }{\partial }_{\mu }a\,{d}^{4}x+{ \mathcal B }+O(\epsilon ),\end{eqnarray} \tag{ 31 }$$

where we have separated off a total divergence term

$$\begin{eqnarray}&&{ \mathcal B }=-\displaystyle \frac{3}{8\pi G}{\int }_{{ \mathcal U }}{\partial }_{\mu }\left\{{\eta }^{\mu \nu }\left[a{\partial }_{\nu }\delta a-\delta a{\partial }_{\nu }a\right]\right\}{d}^{4}x.\end{eqnarray} \tag{ 32 }$$

This term becomes a boundary term, which is usually discarded. We must be more careful in our RW setting, however, because we do not have spatial control of the variations $\delta a(\eta )$. Detailed consideration of the boundary term can be found in Appendix A, where we show that the term should still be discarded. Consistent with this result, we now set ${ \mathcal B }=0$.

We now arrive at the crucial step, which reveals the importance of working directly from the EH action. Only variations consistent with the RW model symmetries are permitted. In other words, the variations $\delta a$ of the scale factor $a(\eta )$ must be functions of conformal time alone: $\delta a(\eta )$. Since the integrand in Equation (31) depends only on η, we have that

$$\begin{eqnarray}&&\delta {S}_{G}=-\displaystyle \frac{3}{4\pi G}{\int }_{{\eta }_{i}}^{{\eta }_{f}}\delta a{\partial }^{\mu }{\partial }_{\mu }a\,d\eta {\int }_{{ \mathcal V }}\,{d}^{3}x+O(\epsilon ).\end{eqnarray} \tag{ 33 }$$

Here we have used Fubini's theorem to write the compactly supported integral as an iterated integral. Note that the spatial integration is just the comoving volume ${ \mathcal V }$, so

$$\begin{eqnarray}&&\delta {S}_{G}=-\displaystyle \frac{3}{4\pi G}{\int }_{{\eta }_{i}}^{{\eta }_{f}}\delta a{\partial }^{\mu }{\partial }_{\mu }a\,{ \mathcal V }d\eta +O(\epsilon ).\end{eqnarray} \tag{ 34 }$$

To compute the variation of the matter action, first note that

$$\begin{eqnarray}&&{T}^{\mu \nu }\delta {g}_{\mu \nu }=2a{\eta }_{\mu \nu }{T}^{\mu \nu }\delta a+O(\epsilon ),\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&=\,2{a}^{-1}\delta a\left[{g}_{\mu \nu }-O(\epsilon )\right]{T}^{\mu \nu }+O(\epsilon ),\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&=\,2{a}^{-1}{{T}^{\mu }}_{\mu }\delta a+O(\epsilon ),\end{eqnarray} \tag{ 37 }$$

where ${{T}^{\mu }}_{\mu }$ is the trace of the stress tensor. Substitution into Equation (17), restricted to ${ \mathcal U }$, gives

$$\begin{eqnarray}&&\delta {S}_{M}={\int }_{{\eta }_{i}}^{{\eta }_{f}}{\int }_{{ \mathcal V }}{a}^{3}{{T}^{\mu }}_{\mu }(\eta ,{\boldsymbol{x}})\delta a\,{d}^{3}x\,d\eta +O(\epsilon ).\end{eqnarray} \tag{ 38 }$$

Again, because a and $\delta a$ depend only on η, we have that

$$\begin{eqnarray}&&\delta {S}_{M}={\int }_{{\eta }_{i}}^{{\eta }_{f}}{a}^{3}\delta a\left[{\int }_{{ \mathcal V }}{{T}^{\mu }}_{\mu }(\eta ,{\boldsymbol{x}})\,{d}^{3}x\right]\,d\eta +O(\epsilon ).\end{eqnarray} \tag{ 39 }$$

The trace of the stress tensor cannot be moved through the integral. The trace continues to be position-dependent, even though the metric at zero order is position-independent, because the stress tensor depends on ${\boldsymbol{\Psi }}(\eta ,{\boldsymbol{x}})$. These fields are not explicitly constrained by our (or any) choice of metric model. The quantity in the square brackets, because of the integration over ${ \mathcal V }$, is formally independent of position. Note, however, that it could be sensitive to the choice of the region ${ \mathcal V }$. This will be elaborated upon in Section 2.3, where we discuss the physical significance of ${ \mathcal V }$.

We now combine the variation of the gravitational action Equation (34) with that of the matter action Equation (39). The action principle as stated in Equation (16) becomes

$$\begin{eqnarray}&&{\int }_{{\eta }_{i}}^{{\eta }_{f}}\delta a\left[-\displaystyle \frac{3}{4\pi G}{\partial }^{\mu }{\partial }_{\mu }a\,{ \mathcal V }+{a}^{3}{\int }_{{ \mathcal V }}{{T}^{\mu }}_{\mu }(\eta ,{\boldsymbol{x}})\,{d}^{3}x\right]d\eta \equiv 0.\end{eqnarray} \tag{ 40 }$$

Applying the Fundamental Theorem of Variational Calculus then gives the consistent equation of motion

$$\begin{eqnarray}&&\displaystyle \frac{3}{4\pi G}{\partial }^{\mu }{\partial }_{\mu }a\,{ \mathcal V }={a}^{3}{\int }_{{ \mathcal V }}{{T}^{\mu }}_{\mu }(\eta ,{\boldsymbol{x}})\,{d}^{3}x.\end{eqnarray} \tag{ 41 }$$

Dividing by the constants on the left naturally reveals

$$\begin{eqnarray}&&{\partial }^{\mu }{\partial }_{\mu }a=\displaystyle \frac{4\pi G}{3}{a}^{3}\displaystyle \frac{1}{{ \mathcal V }}{\int }_{{ \mathcal V }}{{T}^{\mu }}_{\mu }(\eta ,{\boldsymbol{x}})\,{d}^{3}x\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{4\pi G}{3}{a}^{3}{\left\langle {{T}^{\mu }}_{\mu }(\eta ,{\boldsymbol{x}})\right\rangle }_{{ \mathcal V }},\end{eqnarray} \tag{ 43 }$$

which is the spatial-slice average of the stress tensor's trace. This expression is manifestly coordinate-invariant:

• 1.  
  the trace of any tensor is a general coordinate scalar; and
• 2.  
  by Assumption 1, the notion of a spatial slice is geometrically well-defined (e.g., O'neill 1983, Section 12, p. 342).

The spatial integration in this average, and the volume to normalize it, come directly from the action principle. Expanding the derivatives and trace gives Friedmann's acceleration equation

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}a}{d{\eta }^{2}}=\displaystyle \frac{4\pi G}{3}{a}^{3}{\left\langle \rho (\eta ,{\boldsymbol{x}})-\displaystyle \sum _{i=1}^{3}{{ \mathcal P }}_{i}(\eta ,{\boldsymbol{x}})\right\rangle }_{{ \mathcal V }}\end{eqnarray} \tag{ 44 }$$

expressed in conformal time. Though the trace is invariant, $\rho (\eta ,{\boldsymbol{x}})$ and ${{ \mathcal P }}_{i}(\eta ,{\boldsymbol{x}})$ are the diagonal components of the generic stress tensor in the RW preferred coordinate system. This is the coordinate system in which the metric takes the form of Equation (14). One may introduce the definitions

$$\begin{eqnarray}&&\rho (\eta )\equiv \langle \rho (\eta ,{\boldsymbol{x}}){\rangle }_{{ \mathcal V }},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&3{ \mathcal P }(\eta )\equiv {\left\langle \displaystyle \sum _{i=1}^{3}{{ \mathcal P }}_{i}(\eta ,{\boldsymbol{x}})\right\rangle }_{{ \mathcal V }},\end{eqnarray} \tag{ 46 }$$

to arrive at the typical presentation of Friedmann's acceleration equation in conformal time. Note, however, that these defined $\rho (\eta )$ and ${ \mathcal P }(\eta )$ must be spatial-slice averages of the unapproximated quantities.

The spatial average contains the arbitrary closed 3-ball ${ \mathcal V }$. This enters the construction because we must guarantee that the action integral remains well-defined. This means that Equation (44) represents a continuum of Friedmann models.

We may understand the parameter ${ \mathcal V }$ as follows. Fix a radius b and pick some event $P\in { \mathcal U }$. Compute the spatial averages, at all η, of the stress trace over a round ball of radius b centered at the spatial location of P. The resulting source will determine how $a(\eta )$ evolves in conformal time. With the same b, repeat this procedure at a distinct event $Q\in { \mathcal U }$, which lies in the same spatial slice as P. If the Q spatial averages do not agree with the P spatial averages, the dynamics for $a(\eta )$ need not agree. This can happen because matter is not distributed uniformly on small scales. In other words, well-defined field equations for $a(\eta )$ require a volume large enough that the averaged density and pressures are position-independent.

Observations (e.g., Hogg et al. 2005; Scrimgeour et al. 2012; Nadathur 2013) suggest that, at the present epoch, there is a radius

$$\begin{eqnarray}&&{b}_{* }\sim 180\,\mathrm{Mpc},\end{eqnarray} \tag{ 47 }$$

beyond which the averaged quantities are position-independent. We will assume that ${ \mathcal V }$ contains a round ball of at least radius b_*. Note that the specific value of b_* is unimportant for the purposes of our argument. Strong evidence for hierarchical formation of structure tells us that, throughout time, these averaged quantities remain position-independent. In other words, any such ${ \mathcal V }$ suffices for all earlier epochs.

Friedmann's equation, as given in Equation (44) with ${ \mathcal V }$ interpreted according to Section 2.3, resolves all ambiguity in our understanding of how to use Einstein's equations to extract physics on the largest scales. Consistent with our motivation of cosmological backreaction in Section 1, the action includes an integration of degrees of freedom present in the metric over the spatial volume ${ \mathcal V }$. At zero order, however, the only such degree of freedom is the scale factor $a(\eta )$. It is already constrained by the RW model assumption (not the universe) to be position-independent. The required integration produces the spatial-slice volume ${ \mathcal V }$, as seen in Equation (34). In other words, given Assumptions 1 and 2, the action principle produces equations of motion without any cosmological backreaction ${X}^{\mu \nu }(\eta )$.

It is important to note that, given a different starting ansatz, one can arrive at a different cosmological model. This different cosmological model could possibly exhibit a backreaction: additive and position-independent contributions of geometric origin. The presence or absence of such terms is an artifact of the particular model and is logically distinct from the actual metric tensor on ${ \mathcal U }$. Our Assumption 1 reduces to the RW metric when truncated at N = 0, and to the typical covariant linear perturbation metric when truncated at N = 1. What we have shown is that, within these frameworks, there is zero cosmological backreaction ${X}^{\mu \nu }(\eta )$.

Inspection of Equations (44) and (46) reveals that the position-independent pressure ${ \mathcal P }(\eta )$ of any Friedmann model must include all contributions from pressures interior to compact objects on ${ \mathcal V }$. In other words, given Assumptions 1 and 2, using Birkhoff's theorem to construct a cosmological source to Friedmann's equations is inconsistent with the action principle. In general, use of Birkhoff's theorem in cosmological contexts is only appropriate under very restricted settings, which we derive and discuss in Section 3.1. We summarize these statements for clarity.

Our Friedmann's equation is an unambiguous and coordinate-invariant consequence of Assumptions 1, 2, and the EH action.

• 1.  
  There is no backreaction. Inhomogenieties in the spatial distribution of matter do not affect the source to Friedmann's equation.
• 2.  
  Pressures interior to all compact objects contribute to Friedmann's equation.

These consequences arise because the RW metric cannot distinguish spatial regions. Perturbatively consistent equations of motion thus require a zero-order source without explicit or implicit notions of interior.

It is instructive to show how careful application of symmetry considerations also produces the consistent conservation of the stress–energy statement. To derive the appropriate relation, note that a general coordinate scalar, such as the matter action, cannot change under coordinate transformations. It is shown in Weinberg (1972, Section 12.3) that this restriction implies

$$\begin{eqnarray}&&{\rm{\Delta }}{S}_{M}=0={\int }_{{ \mathcal U }}{\epsilon }^{\lambda }{{\rm{\nabla }}}_{\nu }{{T}^{\nu }}_{\lambda }(\eta ,{\boldsymbol{x}})\sqrt{-g}\,{d}^{4}x,\end{eqnarray} \tag{ 48 }$$

where ${\epsilon }^{\lambda }$ is an arbitrary (infinitesimal) vector field. Since ${\epsilon }^{\lambda }$ is arbitrary, the above equation immediately leads to the familiar statement of covariant conservation of stress–energy:

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\nu }{{T}^{\nu }}_{\lambda }(\eta ,{\boldsymbol{x}})=0.\end{eqnarray} \tag{ 49 }$$

Note that we deviate slightly from Weinberg's notation and use Δ to distinguish the change in S_M induced by a coordinate change from any hypothetical change induced by the variational differential δ. We continue to interpret ${ \mathcal V }$ consistent with Section 2.3 For clarity, we will continue to explicitly write the time and position dependence of the generic stress tensor.

At zero order, the position-independent vector fields ${\epsilon }^{\lambda }(\eta )$ are the only ones permissible by RW model symmetries. For consistency, we must restrict ourselves to these vector fields. Note that Equation (48) is true for arbitrary ${\epsilon }^{\lambda }$, so it is true for position-independent ${\epsilon }^{\lambda }(\eta )$. Thus, we may commute though the spatial integral as before:

$$\begin{eqnarray}\begin{array}{rcl}0 & = & {\int }_{{\eta }_{i}}^{{\eta }_{f}}{\epsilon }^{\lambda }(\eta ){a}^{4}{\int }_{{ \mathcal V }}\left[{\partial }_{\nu }{{T}^{\nu }}_{\lambda }(\eta ,{\boldsymbol{x}})+{{{\rm{\Gamma }}}^{\nu }}_{\nu \rho }{{T}^{\rho }}_{\lambda }(\eta ,{\boldsymbol{x}})\right.\\ & & -\,\left.{{{\rm{\Gamma }}}^{\rho }}_{\lambda \nu }{{T}^{\nu }}_{\rho }(\eta ,{\boldsymbol{x}})\right]\,{d}^{3}x\,d\eta +O(\epsilon ).\end{array}\end{eqnarray} \tag{ 50 }$$

This is permissible because the preferred RW frame is a function of η alone. We now expand the derivative term and apply Stokes' theorem on the spatial slice

$$\begin{eqnarray}&&\begin{array}{l}{\int }_{{\eta }_{i}}^{{\eta }_{f}}{\epsilon }^{\lambda }(\eta ){a}^{4}\left\{{\int }_{{ \mathcal V }}\left[{\partial }_{0}{{T}^{0}}_{\lambda }(\eta ,{\boldsymbol{x}})+{{{\rm{\Gamma }}}^{\nu }}_{\nu \rho }{{T}^{\rho }}_{\lambda }(\eta ,{\boldsymbol{x}})\right.\right.\\ \quad -\,\left.{{{\rm{\Gamma }}}^{\rho }}_{\lambda \nu }{{T}^{\nu }}_{\rho }(\eta ,{\boldsymbol{x}})\right]\,{d}^{3}x\\ \quad +\,\left.{\int }_{\partial { \mathcal V }}{{T}^{k}}_{\lambda }(\eta ,{\boldsymbol{x}})\,{d}^{2}{x}_{k}\right\}\,d\eta =0.\end{array}\end{eqnarray} \tag{ 51 }$$

Because we are interested in the conservation condition, set $\lambda =0$ so that the vector ${{T}^{k}}_{0}$ is an energy flow. By Einstein's equations ${{T}^{k}}_{0}\sim {{G}^{k}}_{0}$, but ${{G}^{k}}_{0}\sim O(\epsilon )$ by Assumption 2, so the boundary term must be dropped at this order

$$\begin{eqnarray}\begin{array}{rcl}0 & = & {\int }_{{\eta }_{i}}^{{\eta }_{f}}{a}^{4}{\epsilon }^{0}{\int }_{{ \mathcal V }}\left[{\partial }_{0}{{T}^{0}}_{0}(\eta ,{\boldsymbol{x}})+{{{\rm{\Gamma }}}^{\nu }}_{\nu \rho }{{T}^{\rho }}_{0}(\eta ,{\boldsymbol{x}})\right.\\ & & -\,\left.{{{\rm{\Gamma }}}^{\rho }}_{0\nu }{{T}^{\nu }}_{\rho }(\eta ,{\boldsymbol{x}})\right]\,{d}^{3}x\,d\eta .\end{array}\end{eqnarray} \tag{ 52 }$$

A standard calculation gives

$$\begin{eqnarray}&&{{{\rm{\Gamma }}}^{0}}_{0\nu }={H}_{\eta }{{\delta }^{0}}_{\nu }\qquad {{{\rm{\Gamma }}}^{k}}_{0\nu }={H}_{\eta }{{\delta }^{k}}_{\nu }\qquad {{{\rm{\Gamma }}}^{\nu }}_{\nu \rho }=4{H}_{\eta }{{\delta }^{0}}_{\rho },\end{eqnarray} \tag{ 53 }$$

where ${H}_{\eta }\equiv {a}^{-1}{da}/d\eta$. All of these connection components are dependent on η alone, so that dividing through by ${ \mathcal V }$ gives

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial \eta }{\left\langle {{T}^{0}}_{0}(\eta ,{\boldsymbol{x}})\right\rangle }_{{ \mathcal V }}+3{H}_{\eta }{\left\langle {{T}^{0}}_{0}(\eta ,{\boldsymbol{x}})\right\rangle }_{{ \mathcal V }}-{H}_{\eta }{\left\langle {{T}^{k}}_{k}(\eta ,{\boldsymbol{x}})\right\rangle }_{{ \mathcal V }}=0,\end{eqnarray} \tag{ 54 }$$

which is the appropriate continuum of conservation of stress–energy statements.

The most general RW metric takes the form

$$\begin{eqnarray}&&{{ds}}^{2}=-b{\left(\tau \right)}^{2}d{\tau }^{2}+a{\left(\tau \right)}^{2}{{dx}}^{i}{{dx}}^{j}{\gamma }_{{ij}},\end{eqnarray} \tag{ 55 }$$

where ${\gamma }_{{ij}}$ can be the metric for any 3-space of constant curvature and $b(\tau )$ is a gauge degree of freedom. In the previous computation, we have fixed the gauge by defining

$$\begin{eqnarray}&&b(\tau )\,d\tau \equiv a(\eta )\,d\eta .\end{eqnarray} \tag{ 56 }$$

If one keeps the $b(\tau )$ degree of freedom, and varies with respect to it, the typical Friedmann energy equation is obtained (see Suzuki et al. 1996, Equation (9)). In this way, one can reproduce Equation (54) entirely within the variational formalism. Non-flat spatial slices lead to the same conclusions with respect to computation of the source terms.

The formalism developed in this section has also been applied to the next order in . The advantage is that source terms to the field equations are generated unambiguously. This has been verified by Croker (2018, Section 3) at $O(\epsilon )$ for the scalar modes, where a lengthy calculation in longitudinal gauge reproduces the standard equations.

First, fix some time of interest ${\eta }_{0}$ and consider some later time ${\eta }_{0}+{\rm{\Delta }}\eta$. We may express Equation (14) at ${\eta }_{0}+{\rm{\Delta }}\eta$ as

$$\begin{eqnarray}&&{g}_{\mu \nu }={a}^{2}({\eta }_{0}+{\rm{\Delta }}\eta )\left[{\eta }_{\mu \nu }+...\right],\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}&&=\,\left\{{a}^{2}({\eta }_{0})+2{\rm{\Delta }}\eta \left(a\displaystyle \frac{{da}}{d\eta }\right){| }_{{\eta }_{0}}+...\right\}\left[{\eta }_{\mu \nu }+...\right],\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&=\,{a}^{2}({\eta }_{0})\left\{1+2{\rm{\Delta }}\eta \left(\displaystyle \frac{1}{a}\displaystyle \frac{{da}}{d\eta }\right){| }_{{\eta }_{0}}+...\right\}\left[{\eta }_{\mu \nu }+...\right].\end{eqnarray} \tag{ 59 }$$

If we define a reciprocal timescale

$$\begin{eqnarray}&&{H}_{\eta }^{0}\equiv \left(\displaystyle \frac{1}{a}\displaystyle \frac{{da}}{d\eta }\right){| }_{{\eta }_{0}},\end{eqnarray} \tag{ 60 }$$

then the series expansion of Equation (14) about ${\eta }_{0}$ gives

$$\begin{eqnarray}&&{g}_{\mu \nu }{| }_{u}={a}^{2}({\eta }_{0})\left[{\eta }_{\mu \nu }+\mathop{\mathrm{lim}}\limits_{N\to \infty }\displaystyle \sum _{n=1}^{N}{\epsilon }^{n}{h}_{\mu \nu }^{(n)}(\eta ,{\boldsymbol{x}})+O\left({\rm{\Delta }}\eta {H}_{\eta }^{0}\right)\right].\end{eqnarray} \tag{ 61 }$$

Now choose ${\rm{\Delta }}\eta$ such that the correction is small

$$\begin{eqnarray}&&{\rm{\Delta }}\eta \ll {\left({H}_{\eta }^{0}\right)}^{-1}.\end{eqnarray} \tag{ 62 }$$

For any fixed $\epsilon \lt 1$, the $O\left({\rm{\Delta }}\eta {H}_{\eta }^{0}\right)$ term allows the approximation to be cut off at finite N:

$$\begin{eqnarray}&&N\lt \displaystyle \frac{\mathrm{log}\left[{\rm{\Delta }}\eta {H}_{\eta }^{0}\right]}{\mathrm{log}\epsilon }.\end{eqnarray} \tag{ 63 }$$

For physical reasons, we may choose a spatial volume $v\subset { \mathcal V }$, whose light-crossing time is less than ${\rm{\Delta }}\eta$

$$\begin{eqnarray}&&{v}^{1/3}\lt {\rm{\Delta }}\eta .\end{eqnarray} \tag{ 64 }$$

Note that mathematically, however, the following arguments will not depend on this spatial restriction. We will restrict our attention to $\eta \in [{\eta }_{0},{\eta }_{0}+{\rm{\Delta }}\eta ]$ and define a spacetime domain $u\subset { \mathcal U }$

$$\begin{eqnarray}&&u\equiv [{\eta }_{0},{\eta }_{0}+{\rm{\Delta }}\eta ]\times v.\end{eqnarray} \tag{ 65 }$$

Note that u is a spacetime "hockey puck," which encloses a region of the spatial slice ${ \mathcal V }$ required for the spatial averages in Equation (44). Under a constant rescaling of the representation

$$\begin{eqnarray}&&\eta ^{\prime} \equiv \eta a({\eta }_{0})\qquad {\boldsymbol{x}}^{\prime} \equiv {\boldsymbol{x}}a({\eta }_{0}),\end{eqnarray} \tag{ 66 }$$

Equation (61) becomes

$$\begin{eqnarray}&&{g}_{\mu \nu }{| }_{u}={\eta }_{\mu \nu }+\mathop{\mathrm{lim}}\limits_{N\to \infty }\displaystyle \sum _{n=1}^{N}{\epsilon }^{n}{h}_{\mu \nu }^{(n)}\left(\eta ^{\prime} ,{\boldsymbol{x}}^{\prime} \right)+O\left({\rm{\Delta }}\eta {H}_{\eta }^{0}\right).\end{eqnarray} \tag{ 67 }$$

This is the same metric, so the dynamics are unchanged. In other words, it may be possible to use asymptotically flat models within a universe described by Equation (14). Conclusions drawn from such models, however, are only valid on the restricted spacetime u. In particular, conclusions drawn on u cannot be trivially extended off of u. We summarize for clarity. If the real world satisfies Assumptions 1 and 2, then for short intervals asymptotically flat models are permissible, but only for short intervals.

This is another reason why Birkhoff's theorem cannot be used to infer cosmological behavior. Given a sequence of times $\{{\eta }_{j}\}$ and suitable symmetry, one can use Birkhoff's theorem to construct a sequence of Schwarzschild exterior spacetimes characterized uniquely by a sequence of masses $\{{M}_{j}\}$. Without global knowledge of the actual metric, however, no relation between the M_j can be established.

We emphasize, however, that the converse relation holds. In other words, it is always correct to proceed globally from Equation (14) to locally, as represented by Equation (67). This follows because $u\subset { \mathcal U }$. This means Friedmann's equation remains valid at all events in u, where the metric takes the form given in Equation (67). This means that any source on u, which produces a metric of the form Equation (67), will influence the dynamics of $a(\eta )$ at ${\eta }_{0}$, in the manner given by Equation (44).

For all subsequent discussion, we fix some time ${\eta }_{0}$ and we fix some $\epsilon \lt 1$. Because we have the large-scale solution by virtue of Assumptions 1 and 2, we have a well-defined ${H}_{\eta }^{0}({\eta }_{0})$. Now choose ${\rm{\Delta }}\eta$ and $N\gt 0$ so that Equation (63) is satisfied. We now have a well-defined u given by Equation (65). Strong sources defined in u will the contribute to the dynamics of $a({\eta }_{0})$ through the spatial average over u at ${\eta }_{0}$, as required by Friedmann's equations given in Equation (44). We now study some specific sources of interest.

LIGO (e.g., Abbott et al. 2016) and the Event Horizon Telescope Collaboration et al. (2019) have established the existence of ultracompact objects consistent with the Kerr geometry. In this section, we will establish that the cosmological formalism developed continues to operate near to ultrarelativistic, spinning, sources. Since the Kerr spacetime is asymptotically flat, the considerations of Section 3.1 necessarily restrict use of the Kerr geometry to timescales $\ll 1/{H}_{\eta }^{0}$. In particular, we determine an appropriate submanifold $u^{\prime} \subset u$ defined in Equation (65), where the Kerr spacetime satisfies Assumptions 1 and 2 in the approximate sense of Equation (67). For a convenient reference on the Kerr spacetime, we refer the reader to Visser (2007).

In Kerr–Schild coordinates, the Kerr metric is

$$\begin{eqnarray}&&{{ds}}^{2}=\left({\eta }_{\mu \nu }+\displaystyle \frac{{R}_{s}r{\left({\boldsymbol{x}}\right)}^{3}}{r{\left({\boldsymbol{x}}\right)}^{4}+{A}^{2}{z}^{2}}{{\ell }}_{\mu }{{\ell }}_{\nu }\right){{dx}}^{\mu }{{dx}}^{\nu },\end{eqnarray} \tag{ 68 }$$

where R_s is the Schwarzschild radius for some fixed mass M, A encodes dimensionful information about the spin,⁶ and $r({\boldsymbol{x}})$ is a function of the position coordinates x, y, and z. We will henceforth drop the explicit indication of ${\boldsymbol{x}}$ dependence for r. The covector ${{\ell }}_{\mu }$ is defined to be null with respect to ${\eta }_{\mu \nu }$ and has the form

$$\begin{eqnarray}&&{{\ell }}_{\mu }\equiv \left(1,\displaystyle \frac{{rx}+{Ay}}{{r}^{2}+{A}^{2}},\displaystyle \frac{{ry}-{Ax}}{{r}^{2}+{A}^{2}},\displaystyle \frac{z}{r}\right).\end{eqnarray} \tag{ 69 }$$

Since ${\eta }^{\mu \nu }{{\ell }}_{\mu }{{\ell }}_{\nu }\equiv 0$, it follows that

$$\begin{eqnarray}&&1={\left(\displaystyle \frac{{rx}+{Ay}}{{r}^{2}+{A}^{2}}\right)}^{2}+{\left(\displaystyle \frac{{ry}-{Ax}}{{r}^{2}+{A}^{2}}\right)}^{2}+{\left(\displaystyle \frac{z}{r}\right)}^{2},\end{eqnarray} \tag{ 70 }$$

which implicitly defines r. Since each term on the right side is positive, we see that the components of ${{\ell }}_{\mu }$ satisfy

$$\begin{eqnarray}&&| {{\ell }}_{\mu }| \leqslant 1.\end{eqnarray} \tag{ 71 }$$

It follows immediately that

$$\begin{eqnarray}&&{{\ell }}_{k}{{\ell }}_{j}\leqslant {{\ell }}_{0}| {{\ell }}_{j}| \leqslant {{\ell }}_{0}{{\ell }}_{0}=1.\end{eqnarray} \tag{ 72 }$$

To avoid a naked singularity, we require that

$$\begin{eqnarray}&&A\lt M,\end{eqnarray} \tag{ 73 }$$

which maintains real-valued solutions for the horizon surfaces.

By Equation (72), a necessary and sufficient condition to write the metric in the form of Equation (67) is

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{s}{r}^{3}}{{r}^{4}+{A}^{2}{z}^{2}}\leqslant \epsilon .\end{eqnarray} \tag{ 74 }$$

This restriction is saturated when

$$\begin{eqnarray}&&\epsilon {r}^{4}-{R}_{s}{r}^{3}+\epsilon {A}^{2}{z}^{2}=0.\end{eqnarray} \tag{ 75 }$$

By Descartes' rule of signs, this equation implicitly defines two (or zero) non-intersecting surfaces ${r}_{\pm }({\boldsymbol{x}})$. As $\epsilon \to 1$, ${r}_{+}$ and ${r}_{-}$ become, respectively, the outer and inner ergosurfaces of Kerr. Define a submanifold $u^{\prime} \subset u$ such that

$$\begin{eqnarray}&&u^{\prime} \equiv [{\eta }_{0},{\eta }_{0}+{\rm{\Delta }}\eta ]\times v/\left\{{\boldsymbol{x}}:\epsilon \left({r}_{+}^{4}+{A}^{2}{z}^{2}\right)-{R}_{s}{r}_{+}^{3}\lt 0\right\}.\end{eqnarray} \tag{ 76 }$$

Then $u^{\prime}$ has an additional boundary, which becomes the outer ergosurface as $\epsilon \to 1$. We now establish that this additional boundary will not alter Friedmann's equations. As discussed in Appendix A, the dynamically relevant contribution to Friedmann's acceleration equation from this boundary is the second term in Equation (129). This term continues to vanish at zero order on $u^{\prime}$ because of the RW model symmetries. As shown in Equation (51), an additional term also appears in the conservation equation

$$\begin{eqnarray}&&\int {\epsilon }^{0}(\eta ){a}^{4}\left\{{\int }_{\partial u^{\prime} }{{T}^{k}}_{0}(\eta ,{\boldsymbol{x}})\,{d}^{2}{x}_{k}\right\}\,d\eta .\end{eqnarray} \tag{ 77 }$$

Since Kerr is a vacuum solution, however, ${{T}^{k}}_{0}=0$ on $\partial u^{\prime}$ and this term also vanishes.

Now that we have a suitable domain, we may write Equation (68) in the form of Equation (67):

$$\begin{eqnarray}&&{h}_{\mu \nu }^{(1)}\equiv \left(\displaystyle \frac{{R}_{s}}{r}\right){\left[1+\displaystyle \frac{{A}^{2}{z}^{2}}{{r}^{4}}\right]}^{-1}{{\ell }}_{\mu }{{\ell }}_{\nu },\end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray}&&{h}_{\mu \nu }^{(n\gt 1)}\equiv 0.\end{eqnarray} \tag{ 79 }$$

In summary, we have established a domain $u^{\prime}$ where the Kerr solution can be accommodated under Assumptions 1 and 2. This implies that the spatial averages appearing in Friedmann's equation remain valid in the ultrarelativistic vicinity of a Kerr BH. The domain $u^{\prime}$ can be taken to include all but the region enclosed by the outer ergosurface. Consequently, for the domain we have constructed, a Kerr BH contributes nothing to Friedmann's equations. Local observers in $u^{\prime}$ will still perceive a Kerr BH with mass M and spin A for ${r}_{+}\lt r\lt {v}^{1/3}$ and $\eta \in [{\eta }_{0},{\eta }_{0}+{\rm{\Delta }}\eta ]$. This is unsurprising, because the cutting procedure we have performed is mathematical. Visser (2007) emphasizes that the inner horizon and the enclosed inner ergoregion are extremely pathological and should not be regarded as physically relevant. Our result suggests that, additionally, the entire region below the outer ergosurface should be replaced with a distinct interior solution.

In this section, we show that the interior region of Schwarzschild's constant-density sphere satisfies Assumptions 1 and 2. The exterior region of Schwarzschild's constant-density sphere is an A = 0 Kerr solution (i.e., Schwarzschild's BH), which has already been treated. The constant-density solution can be chosen to dominate the actual energy densities for many spherically symmetric, static, sources. We will show that, for such objects with physical radius $R\gt 3{GM}$, all energy densities and pressures in u will influence $a(\eta )$ according to Equation (44). This will establish that localized pressures, interior to a large class of compact and relativistic objects, contribute to the global Friedmann average.

The Schwarzschild constant-density sphere on ${B}_{3}(0,R)$ may be expressed as

$$\begin{eqnarray}&&{{ds}}^{2}=-\exp (2{\rm{\Phi }})\,d{\eta }^{2}+\exp (2{\rm{\Lambda }})\,{{dr}}^{2}+{r}^{2}d{{\rm{\Omega }}}^{2}.\end{eqnarray} \tag{ 80 }$$

The functions Φ and Λ are defined in terms of the Schwarzschild radius R_s and physical radius R of the object

$$\begin{eqnarray}&&\exp \left(2{\rm{\Phi }}\right)\equiv \displaystyle \frac{1}{4}{\left[3\sqrt{1-\displaystyle \frac{{R}_{s}}{R}}-\sqrt{1-\displaystyle \frac{{R}_{s}}{R}{\left(\displaystyle \frac{r}{R}\right)}^{2}}\right]}^{2},\end{eqnarray} \tag{ 81 }$$

$$\begin{eqnarray}&&\exp \left(2{\rm{\Lambda }}\right)\equiv {\left[1-{\left(\displaystyle \frac{r}{R}\right)}^{2}\displaystyle \frac{{R}_{s}}{R}\right]}^{-1}.\end{eqnarray} \tag{ 82 }$$

We will show that there exists a natural such that the coefficients of a Taylor expansion of Equation (80) take the form of Equation (67).

Denote the Taylor expansion coefficients of $\sqrt{1-x}$ on $0\leqslant x\lt 1$ by q_j:

$$\begin{eqnarray}&&\sqrt{1-x}=1-\displaystyle \frac{x}{2}-\displaystyle \sum _{j=2}^{\infty }\displaystyle \frac{(2j-3)!}{{2}^{2(j-1)}j!(j-2)!}{x}^{j},\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&\equiv \,{q}_{0}-{q}_{1}x-\displaystyle \sum _{j=2}^{\infty }{q}_{j}{x}^{j}.\end{eqnarray} \tag{ 84 }$$

Furthermore, extend the notation to negative indices by

$$\begin{eqnarray}&&{q}_{k}\equiv 0\qquad k\lt 0.\end{eqnarray} \tag{ 85 }$$

Then it can be shown that the expansion of Equation (81) in R_s/R becomes

$$\begin{eqnarray}\begin{array}{rcl}\exp (2{\rm{\Phi }}) & = & 1-\displaystyle \frac{{R}_{s}}{R}\left[\displaystyle \frac{3}{2}-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{r}{R}\right)}^{2}\right]\\ & & -\,\displaystyle \frac{3}{2}\displaystyle \sum _{n=2}^{\infty }{\left(\displaystyle \frac{{R}_{s}}{R}\right)}^{n}\left[\displaystyle \sum _{j=0}^{\infty }{q}_{j}{q}_{n-j}{\left(\displaystyle \frac{r}{R}\right)}^{2j}\right].\end{array}\end{eqnarray} \tag{ 86 }$$

The terms in these infinite sums can be commuted and rearranged because $\sqrt{1-x}$ converges uniformly for $| x| \lt 1$. To guarantee that the coefficient of the linear term remains bounded by 1 on $r\in [0,R]$, we must have

$$\begin{eqnarray}&&\epsilon \geqslant \displaystyle \frac{3{R}_{s}}{2R}.\end{eqnarray} \tag{ 87 }$$

Since $\epsilon \lt 1$, this becomes a restriction on the object's physical radius

$$\begin{eqnarray}&&R\gt 3{GM}.\end{eqnarray} \tag{ 88 }$$

Note that this is the radius of the innermost stable orbit for a photon.⁷ Saturating this bound, we find

$$\begin{eqnarray}\begin{array}{rcl}\exp (2{\rm{\Phi }}) & = & 1-\epsilon \left[1-\displaystyle \frac{1}{3}{\left(\displaystyle \frac{r}{R}\right)}^{2}\right]\\ & & -\,\displaystyle \sum _{n=2}^{\infty }{\epsilon }^{n}\left[\displaystyle \sum _{j=0}^{\infty }{q}_{j}{q}_{n-j}{\left(\displaystyle \frac{2}{3}\right)}^{n-1}{\left(\displaystyle \frac{r}{R}\right)}^{2j}\right].\end{array}\end{eqnarray} \tag{ 89 }$$

To bound the remaining coefficients of ${\epsilon }^{n}$, note that for each coefficient in n,

$$\begin{eqnarray}&&\displaystyle \sum _{j=0}^{\infty }{q}_{j}{q}_{n-j}{\left(\displaystyle \frac{2}{3}\right)}^{n-1}{\left(\displaystyle \frac{r}{R}\right)}^{2j}\lt \displaystyle \sum _{j=0}^{\infty }{q}_{j}{\left(\displaystyle \frac{r}{R}\right)}^{2j}=\sqrt{1-{\left(\displaystyle \frac{r}{R}\right)}^{2}}.\end{eqnarray} \tag{ 90 }$$

The inequality follows because ${q}_{n-j}\lt 1$ for all j. Since $\sqrt{1-x}\leqslant 1$ for $x\in [0,1]$, each coefficient of ${\epsilon }^{n}$ for $n\geqslant 2$ is bounded above by 1 on $r\in [0,R]$. This was to be shown. The result for $\exp (2{\rm{\Lambda }})$ follows immediately from use of the binomial expansion, followed by substitution of Equation (87). In Appendix B, we establish that the Schwarzschild constant-density sphere also satisfies Assumption 2.

In summary, we have established that the Schwarzschild constant-density sphere is of the form Equation (67). In other words, all pressures interior to any astrophysical object bounded by the Schwarzschild constant-density sphere solution will contribute to the cosmological average in Equation (44).

In this section, we show that the interior of an isolated de-Sitter sphere satisfies Assumptions 1 and 2. The exterior region of such a sphere is an A = 0 Kerr solution (i.e., Schwarzschild's BH), which has already been treated. The isolated de-Sitter sphere is the simplest model of a GEneric Object of Dark Energy (GEODE). Such objects, like the solution of Dymnikova (1992) or the gravastar of Mazur & Mottola (2015), have been proposed as possible BH replacements. Related, but dynamic, GEODEs called "vacuum bubbles," have also been considered as inflationary relics by Berezin et al. (1987). We will demonstrate shortly that the strong negative pressure inside the de-Sitter sphere will influence $a(\eta )$ according to Equation (44). This will establish that a physically realistic GEODE could contribute to the cosmological ${ \mathcal P }=-\rho$.

Consider a de-Sitter patch in static coordinates

$$\begin{eqnarray}&&{{ds}}^{2}=-\left(1-\displaystyle \frac{{r}^{{\prime} 2}}{{R}_{s}^{2}}\right)d{\eta }^{{\prime} 2}+{\left(1-\displaystyle \frac{{r}^{{\prime} 2}}{{R}_{s}^{2}}\right)}^{-1}{{dr}}^{{\prime} 2}+{r}^{{\prime} 2}d{{\rm{\Omega }}}^{2}.\end{eqnarray} \tag{ 91 }$$

We have introduced primes on the $\eta ^{\prime}$ and $r^{\prime}$ coordinates for reasons that will become apparent. Denote these coordinates by ${\rm{\Xi }}^{\prime}$. This patch admits a timelike Killing vector field for $r^{\prime} \lt {R}_{s}$. This means that a sphere of Dark Energy with radius $R\lt {R}_{s}$ is static, regardless of the coordinate representation.

Again, expand in a series, then multiply and divide by :

$$\begin{eqnarray}\begin{array}{rcl}{g}_{\mu \nu } & = & {\eta }_{\mu \nu }(\eta ^{\prime} ,{\boldsymbol{x}}^{\prime} )+{\epsilon }^{2}{\left(\displaystyle \frac{r^{\prime} }{{R}_{s}\epsilon }\right)}^{2}{{\delta }^{0}}_{\mu }{{\delta }^{0}}_{\nu }\\ & & +\ \displaystyle \sum _{n=1}^{\infty }{\epsilon }^{2n}{\left(\displaystyle \frac{r^{\prime} }{{R}_{s}\epsilon }\right)}^{2n}{{\delta }^{1}}_{\mu }{{\delta }^{1}}_{\nu }.\end{array}\end{eqnarray} \tag{ 92 }$$

We see that $\left|{h}_{\mu \nu }^{(n)}\right|\leqslant 1$, so Equation (92) satisfies Equation (67), provided that

$$\begin{eqnarray}&&r^{\prime} \leqslant {R}_{s}\epsilon .\end{eqnarray} \tag{ 93 }$$

In other words, Equation (91) has a convergent representation on ${B}_{3}(0,{R}_{s}\epsilon )$. Contrast this situation with the Schwarzschild BH, which converges on $v/{B}_{3}(0,{R}_{s}{\epsilon }^{-1})$. Because Einstein's equations hold everywhere on ${ \mathcal U }$, we know there exist coordinates ζ, which cover $r^{\prime} \gt {R}_{s}\epsilon$. This means there exists $z\gt 0$ such that

$$\begin{eqnarray}&&{\mathscr{O}}^{\prime} \equiv \mathrm{dom}({\rm{\Xi }}^{\prime} )\cap \mathrm{dom}(\zeta ),\end{eqnarray} \tag{ 94 }$$

$$\begin{eqnarray}&&=\ \left({\eta }_{0}^{{\prime} },{\eta }_{0}^{{\prime} }+{\rm{\Delta }}\eta ^{\prime} \right)\times \left({R}_{s}\epsilon -z,{R}_{s}\epsilon \right)\times {S}^{2},\end{eqnarray} \tag{ 95 }$$

is the overlap between the ${\rm{\Xi }}^{\prime}$ and ζ charts. We will continue to use the ${\rm{\Xi }}^{\prime}$ chart on this region.

The relation between the Ξ coordinates of the exterior Schwarzschild spacetime and the ${\rm{\Xi }}^{\prime}$ coordinates of the de-Sitter sphere is determined by Mazur & Mottola (2015, Equations (5.1)–(5.3)). They find that Ξ and ζ are the same, and that ${\rm{\Xi }}^{\prime}$ satisfies

$$\begin{eqnarray}&&r^{\prime} =r,\end{eqnarray} \tag{ 96 }$$

$$\begin{eqnarray}&&d\eta ^{\prime} =2\,d\eta .\end{eqnarray} \tag{ 97 }$$

In other words, the spatial slices are unaltered, but time interior to the sphere runs twice as fast. This means that our static approximation is only valid for ${\rm{\Delta }}\eta /2$. Apart from a reduction in v to $v/8$, nothing changes with respect to the spatial-slice integration. We may integrate Equation (97) to find

$$\begin{eqnarray}&&\eta ^{\prime} (\eta )=2\eta +C,\end{eqnarray} \tag{ 98 }$$

where C is an integration constant. Note that we may choose $C\equiv -{\eta }_{0}$ so that

$$\begin{eqnarray}&&a\left[\eta ^{\prime} ({\eta }_{0})\right]=a({\eta }_{0}).\end{eqnarray} \tag{ 99 }$$

This shows that the unit redefinition used to write Equation (67) remains consistent inside the sphere.

We now cut out the $r\gt {R}_{s}\epsilon -z/2$ region, leaving a boundary at the cut radius. It again follows that Equation (129) continues to vanish at zero order. The required conservation condition in Equation (50) continues to hold because ${{T}^{\mu }}_{\nu }\,\propto {{\delta }^{\mu }}_{\nu }$, so the contribution vanishes identically. Consider the domain $u^{\prime\prime} \subset u$

$$\begin{eqnarray}\begin{array}{rcl}u^{\prime\prime} & \equiv & \left[{\eta }_{0},{\eta }_{0}+\displaystyle \frac{{\rm{\Delta }}\eta }{2}\right]\\ & & \times \ \displaystyle \frac{v}{8}/\left\{\left[{R}_{s}\epsilon -\displaystyle \frac{z}{2},\displaystyle \frac{{R}_{s}}{\epsilon }+\displaystyle \frac{z}{2}\right]\times {S}^{2}\right\}.\end{array}\end{eqnarray} \tag{ 100 }$$

We have established that Equation (67) is satisfied everywhere on $u^{\prime\prime}$ and that there are no additional boundary contributions to Equation (44). Finally, static de-Sitter space satisfies Assumption 2 because its stress tensor is constant. By adjusting z and in Equation (100), we conclude that nearly all of the de-Sitter region contributes ${ \mathcal P }=-\rho$ to Equation (44).

The cosmological evolution of local energies is a necessary consequence of Assumptions 1 and 2, and the EH action. The essential question is now observational: can this effect be measured? Because this section is observational, it is sometimes clearer to use redshift z as a time variable. We thus convert Equation (108) to redshift

$$\begin{eqnarray}&&{E}_{i}\propto {\left(1+z\right)}^{3w},\end{eqnarray} \tag{ 109 }$$

using that $a={(1+z)}^{-1}$. Where appropriate, we use Planck collaboration cosmological parameters from Aghanim et al. (2018): ${H}_{0}=67.4\pm 0.5\,\mathrm{km}\,{\mathrm{Mpc}}^{-1}\,{{\rm{s}}}^{-1}$, ${{\rm{\Omega }}}_{b}=(2.24\pm 0.01)\,\times {10}^{-2}{h}^{-2}\simeq 0.05$, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}\equiv 1-{{\rm{\Omega }}}_{m}=0.685\pm 0.007$. We now consider specific model objects relevant for astrophysical observation, using fixed w approximations. We proceed in order of increasingly relativistic w.

4.1.1. Stars

Consider a population of typical stars, at fixed (comoving) coordinate positions. Each typical star will contribute an extremely small positive pressure to the cosmological average. This follows because the pressure is everywhere positive within a star, so an integral over the stellar pressure cannot vanish. This is true even in simplified stellar models, where fluid packets are radially static. Each packet's contribution to the pressure is non-zero because gravitational momentum flux is not included in ${T}_{\mu \nu }$. The equation of state w for a typical star can be approximated from the ideal gas law,

$$\begin{eqnarray}&&{w}_{\mathrm{star}}\sim \displaystyle \frac{{kT}}{{m}_{{\rm{p}}}},\end{eqnarray} \tag{ 110 }$$

where k is Boltzmann's constant, T is the temperature of the star, and m_p is the proton mass. To develop an upper bound, set $T\sim {10}^{6}\,{\rm{K}}$, representative of core temperatures. Then

$$\begin{eqnarray}&&{w}_{\mathrm{star}}\sim {10}^{-7}.\end{eqnarray} \tag{ 111 }$$

Thus, the energy (as perceived by RW observers) of Sun-like stars cosmologically evolves as

$$\begin{eqnarray}&&{E}_{\mathrm{star}}\propto {\left(1+z\right)}^{3\times {10}^{-7}}.\end{eqnarray} \tag{ 112 }$$

How large is this effect? Suppose a star is produced at z_i = 2 and observed at z_f = 0. Define the fractional shift in energy as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}E}{E}\equiv \displaystyle \frac{E({z}_{f})-E({z}_{i})}{E({z}_{i})},\end{eqnarray} \tag{ 113 }$$

then

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{E}_{\mathrm{star}}}{{E}_{\mathrm{star}}}={\left(\displaystyle \frac{1}{1+2}\right)}^{3\times {10}^{-7}}-1\simeq -{10}^{-7}.\end{eqnarray} \tag{ 114 }$$

This change occurs over $10\,\mathrm{Gyr}$ from z_i = 2 to z_f = 0. In other words, it is dominated by other stellar processes and is thus unobservable. This also establishes that the effect is unobservable for any other material with $w\lt {w}_{\mathrm{star}}$.

What is the reciprocal effect on the zero-order expansion? Stars contribute $\sim 2{{\rm{\Omega }}}_{b}/5$ to Friedmann's equation. The cumulative adjustment to $\rho (a)$ from a conservative first light of z_i = 40 to z_f = 0 is then $\sim -{10}^{-8}$. The effect is unobservable.

4.1.2. Galaxy Clusters

Consider a population of clusters, at fixed (comoving) coordinate positions. The velocity dispersion of clusters is roughly bounded by $\sigma \sim {10}^{3}\,\mathrm{km}\,{{\rm{s}}}^{-1}$. This implies a ${w}_{\mathrm{cluster}}\,\sim {10}^{-5}$ by ideal gas arguments. Thus, the energy (as perceived by RW observers) of clusters cosmologically evolves as

$$\begin{eqnarray}&&{E}_{\mathrm{cluster}}\propto {\left(1+z\right)}^{3\times {10}^{-5}}.\end{eqnarray} \tag{ 115 }$$

How large is this effect? A typical length scale for a cluster is $\sim 1\,\mathrm{Mpc}$, giving a light-crossing time of $\sim 3\,\mathrm{Myr}$. In redshift, this could be between z_i = 1.001 and z_f = 1. During this time, the fractional shift in energy of the cluster is

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{E}_{\mathrm{cluster}}}{{E}_{\mathrm{cluster}}}={\left(\displaystyle \frac{1+1}{1+1.001}\right)}^{3\times {10}^{-5}}-1\simeq -{10}^{-8}.\end{eqnarray} \tag{ 116 }$$

For comparison, during a photon transit completing at z_f = 2, the cluster energy shifts less than 10⁻⁸. In other words, the correction to any integrated Sachs-Wolfe effect through clusters is unobservable. Similarly, corrections to gravitational lensing are unobservable.

What is the reciprocal effect on the zero-order expansion? The energy density in galaxy clusters is certainly bounded above by ${{\rm{\Omega }}}_{b}$. The cumulative adjustment to $\rho (a)$ from a conservative first cluster of z_i = 4 to z_f = 0 is then $\sim -{10}^{-6}$. The effect is unobservable at present and upcoming sensitivities.

4.1.3. GEneric Objects of Dark Energy (GEODEs)

As discussed in Section 3.4, GEODEs are explicit GR solutions, which schematically resemble the static de-Sitter sphere. Before gravitational-wave observations of ultrarelativistic object mergers, GEODEs were of theoretical interest because they are often free of physical singularities and horizons. In other words, they are regular solutions for gravitational-collapse remnants, which resolve the BH Information Paradox.

Consider a population of GEODEs positioned at fixed (comoving) coordinate positions. Recall that the equation of state of a de-Sitter sphere is

$$\begin{eqnarray}&&{w}_{\mathrm{dS}}=-1.\end{eqnarray} \tag{ 117 }$$

This equation of state is appropriate if the GEODE edge does not contribute significantly to the population-averaged ${ \mathcal P }$. GEODE material is maximally relativistic: it saturates the dominant energy condition. The material is also under tension, instead of under compression: the sign of w is inverted from all previously considered cases. This leads to very useful consequences.

The energy (as perceived by RW observers) of the GEODEs cosmologically evolves as

$$\begin{eqnarray}&&{E}_{\mathrm{dS}}\propto {\left(1+z\right)}^{-3}.\end{eqnarray} \tag{ 118 }$$

Unlike the cases previously considered, this shift is significant and acts to amplify the energy. In other words, GEODEs cosmologically blueshift. The effect is analogous to the photon redshift.

How large is this effect? Consider a $3{M}_{\odot }$ GEODE formed from stellar gravitational collapse at z = 1.5, near the peak in comoving star formation rate density (e.g., Madau & Dickinson 2014). Its mass, observed in a binary merger at z = 0.1, will be

$$\begin{eqnarray}&&{E}_{\mathrm{dS}}=3{M}_{\odot }{\left(\displaystyle \frac{1+1.5}{1+0.1}\right)}^{3}=35.2{M}_{\odot }.\end{eqnarray} \tag{ 119 }$$

This compares favorably with the masses observed by LIGO. A thorough exploration of the observational signatures of a remnant GEODE population, with respect to gravitational-wave observatories, is given by Croker et al. (2019). They find that the scenario is consistent with the LIGO GWTC-1 observed mass function. The cosmological blueshift of a GEODE can achieve the required masses with the standard common envelope binary formation channel. It does not require low-metallicity regions or prohibitively large progenitor stars.

What is the reciprocal effect on the zero-order expansion? Suppose 1% of all stellar material collapsed to Population III GEODEs, instead of traditional black holes,

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{collapse}}\sim 2\times {10}^{-4}.\end{eqnarray} \tag{ 120 }$$

For simplicity, assume that this happens at $z\sim 14$: after the Dark Ages but before reionization. The shift in energy of this population, observed at the present day, gives rise to

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{GEODE}}={{\rm{\Omega }}}_{\mathrm{collapse}}{\left(1+14\right)}^{3},\end{eqnarray} \tag{ 121 }$$

$$\begin{eqnarray}&&=\,0.675.\end{eqnarray} \tag{ 122 }$$

Note that ${\rho }_{{\rm{G}}{\rm{E}}{\rm{O}}{\rm{D}}{\rm{E}}}$ is essentially constant within Friedmann's equations. This follows because ${w}_{\mathrm{dS}}=-1$, so that Equation (105) gives no time evolution of the physical density. In other words, the GEODEs dilute in number ∝a⁻³ while gaining in mass ∝a³. The result is an apparent cosmological constant, which compares favorably with the observed ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$.

A first exploration of the observational signatures of Population III GEODEs, with respect to cosmological observables, is given by Croker (2018). They find that such a scenario is flexible enough to reproduce the observed late-time accelerated expansion, resolving the coincidence problem.

We have shown that non-relativistic material exhibits no observable shift. For relativistic objects, however, the energy shift can lead to measurable consequences. Neutron star (NS) material is highly relativistic (e.g., Abbott et al. 2018), with an object-averaged $0.05\lesssim {w}_{\mathrm{NS}}\lesssim 0.1$. Given that binary pulsar dynamics are often measured to exquisite precision, it may be possible to measure the cosmological shift in pulsar systems.

Since the cosmological shift for positive w_NS appears as an energy loss, one can use a result of Damour & Taylor (1991, Equation (4.1)) to estimate the shift in orbital period decay

$$\begin{eqnarray}&&{\rm{\Delta }}{\dot{P}}_{{\rm{b}},\cos }=+6{w}_{\mathrm{NS}}{{HP}}_{{\rm{b}}}.\end{eqnarray} \tag{ 123 }$$

Using P_b as reported by Weisberg et al. (2010, Table 3) for the Hulse–Taylor binary pulsar, we find

$$\begin{eqnarray}&&{\rm{\Delta }}{\dot{P}}_{{\rm{b}},\cos }=(2.7\pm 0.9)\times {10}^{-14},\end{eqnarray} \tag{ 124 }$$

at the present day. How does this shift compare to current measurements of the Hulse–Taylor system? From Weisberg et al. (2010, Equation (5) and Table 3), the kinematically corrected orbital period decay is

$$\begin{eqnarray}&&{\dot{P}}_{{\rm{b}}}-{\rm{\Delta }}{\dot{P}}_{{\rm{b}},\mathrm{gal}}=(-2.396\pm 0.005)\times {10}^{-12}.\end{eqnarray} \tag{ 125 }$$

The shift in orbital period decay at the present, due to cosmological evolution of NS energies, would appear to be $\sim 1 \%$.

The central value of Equation (125) agrees with the GR radiative loss prediction to $0.2 \%$. So a cosmological effect of the estimated magnitude has not been observed. Since Damour & Taylor (1991) assumed a flat spacetime, however, it follows from Equation (67) that the metric is known only to $O({{HP}}_{{\rm{b}}})$. In other words, the estimated shift given by Equation (123) is dominated by the error introduced by working under an asymptotically flat approximation within an RW cosmology. An investigation of binary pulsar systems with cosmologically evolving mass is the topic of future work.

The contribution to the cosmologically averaged quantities, and the necessary local evolution, is highly dependent on the object model. An instructive example is given by Dymnikova (1992). Dymnikova's object contains a very nearly Dark Energy interior, but a spatial average over her object gives exactly ${ \mathcal P }=0$. In her model, the "skin" acts like a vacuum vessel: positive pressures in the skin maintain the static character of the object. A population of Dymnikova's objects would behave as cosmological dust, and exhibit no cosmological shift.

At present, no known GEODE solution strongly rotates. As pointed out by J. D. Bjorken (2018, private discussion), vacuum cannot rotate because it has no privileged directions. Given clear evidence of Kerr exterior spacetimes, all of the spin must reside in material outside of, but extremely near, the de-Sitter core. GEODE solutions are often constructed to have a Schwarzschild exterior spacetime and thus cannot evolve cosmologically. Given the considerations of Section 3.1, and the consequences in Section 4, it seems essential to have adiabatically evolving object solutions that spatially asymptote to an arbitrary RW geometry. The construction of a realistic GEODE model is an open question that is beyond the scope of this paper.

Setting r = 0 gives

$$\begin{eqnarray}&&\displaystyle \frac{{ \mathcal P }(0)}{\rho }=-1+\displaystyle \frac{2}{3}\displaystyle \sum _{n=0}^{\infty }{\left[1-\displaystyle \frac{2\epsilon }{3}\right]}^{-n/2}\displaystyle \frac{1}{{3}^{n}}.\end{eqnarray} \tag{ 152 }$$

This expression satisfies the requirements for a uniformly convergent generalized binomial expansion:

$$\begin{eqnarray}&&-1+\displaystyle \frac{2}{3}\displaystyle \sum _{n=0}^{\infty }\displaystyle \sum _{k=0}^{\infty }{\epsilon }^{k}\left(\displaystyle \genfrac{}{}{0em}{}{n/2+k-1}{n/2-1}\right){\left(\displaystyle \frac{2}{3}\right)}^{k}\displaystyle \frac{1}{{3}^{n}}.\end{eqnarray} \tag{ 153 }$$

Since both sums are uniformly convergent, we commute terms:

$$\begin{eqnarray}&&-1+\displaystyle \frac{2}{3}\displaystyle \sum _{k=0}^{\infty }{\epsilon }^{k}\left[\displaystyle \sum _{n=0}^{\infty }\left(\displaystyle \genfrac{}{}{0em}{}{n/2+k-1}{n/2-1}\right){\left(\displaystyle \frac{2}{3}\right)}^{k}\displaystyle \frac{1}{{3}^{n}}\right].\end{eqnarray} \tag{ 154 }$$

By definition of the binomial coefficients,

$$\begin{eqnarray}&&\left(\displaystyle \genfrac{}{}{0em}{}{n/2+k-1}{n/2-1}\right)=\displaystyle \frac{1}{{2}^{k}k!}\displaystyle \prod _{j=0}^{k-1}\left(n+2j\right),\end{eqnarray} \tag{ 155 }$$

we have that

$$\begin{eqnarray}&&-1+\displaystyle \sum _{k=0}^{\infty }{\epsilon }^{k}\left[\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{2}{{3}^{(n+k+1)}k!}\displaystyle \prod _{j=0}^{k-1}\left(n+2j\right)\right].\end{eqnarray} \tag{ 156 }$$

To demonstrate that ${G}_{{ss}}^{(k)}(0)\leqslant 1$, we establish that c_k, defined as

$$\begin{eqnarray}&&{c}_{k}\equiv \displaystyle \sum _{n=0}^{\infty }{c}_{n}(k)\equiv \displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{2}{{3}^{(n+k+1)}k!}\displaystyle \prod _{j=0}^{k-1}\left(n+2j\right),\end{eqnarray} \tag{ 157 }$$

monotonically decreases in k. First, note that all c_n(k) are positive. We write the ratio of successive c_n(k), at fixed n

$$\begin{eqnarray}&&\displaystyle \frac{{c}_{n}(k+1)}{{c}_{n}(k)}=\displaystyle \frac{n+2k}{3(k+1)}.\end{eqnarray} \tag{ 158 }$$

Note that the ratio is unity when $n=k+3$. For terms with $n\lt k+3$, individual c_n(k) monotonically decrease in k by the ratio test. We now partition off these terms:

$$\begin{eqnarray}&&{c}_{k}\equiv \displaystyle \sum _{n=0}^{k+1}{c}_{n}(k)+\displaystyle \sum _{n=k+2}^{\infty }\displaystyle \frac{2}{{3}^{(n+k+1)}k!}\displaystyle \prod _{j=0}^{k-1}\left(n+2j\right).\end{eqnarray} \tag{ 159 }$$

We need not assert that the finite sum itself is monotonically decreasing in k. We are only bookkeeping individual terms in the sum over n that monotonically decrease in k. As k increases, the number of such terms increases. It suffices to show that the summed tail contribution is monotonically decreasing in k. Once this is established, then all contributions to c_k monotonically decrease in k. This will establish that the c_k monotonically decrease in k.

We now proceed to show that the tail monotonically decreases in k. First, re-index the infinite series,

$$\begin{eqnarray}&&{c}_{k}\equiv \displaystyle \sum _{n=0}^{k+1}{c}_{n}(k)+\displaystyle \sum _{m=0}^{\infty }\displaystyle \frac{2}{{3}^{(m+2k+3)}k!}\displaystyle \prod _{j=0}^{k-1}\left(m+k+2+2j\right),\end{eqnarray} \tag{ 160 }$$

and introduce notation for the coefficients of this infinite series

$$\begin{eqnarray}&&\displaystyle \sum _{m=0}^{\infty }{d}_{m}(k)\equiv \displaystyle \sum _{m=0}^{\infty }\displaystyle \frac{2}{{3}^{(m+2k+3)}k!}\displaystyle \prod _{j=0}^{k-1}\left(m+k+2+2j\right).\end{eqnarray} \tag{ 161 }$$

Expanding the above product gives

$$\begin{eqnarray}&&[m+(3k)][m+(3k-2)]...[m+(3k-2k-2)],\end{eqnarray} \tag{ 162 }$$

where we have grouped each multiplicand as a binomial in m and terms decrementing k. Note that there are ${2}^{k}$ grouped terms in its full expansion. For d_n(k) where $m\geqslant 3k$, ${2}^{k}{m}^{k}$ is then an upper bound on the product. The tail of this series contribution to Equation (160) is then

$$\begin{eqnarray}&&\displaystyle \sum _{m=3k}^{\infty }{d}_{m}(k)\lt \displaystyle \sum _{m=3k}^{\infty }\displaystyle \frac{{2}^{k+1}{m}^{k}}{{3}^{(m+2k+3)}k!}=\displaystyle \frac{{2}^{k+1}}{{3}^{(2k+3)}k!}\displaystyle \sum _{m=3k}^{\infty }\displaystyle \frac{{m}^{k}}{{3}^{m}}.\end{eqnarray} \tag{ 163 }$$

We bound the series with the integral

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{m=3k}^{\infty }\displaystyle \frac{{m}^{k}}{{3}^{m}}\lt \displaystyle \sum _{m=0}^{\infty }\displaystyle \frac{{m}^{k}}{{3}^{m}}\lt \displaystyle \frac{1}{{\left(\mathrm{ln}3\right)}^{k+1}}{\int }_{0}^{\infty }{p}^{k}{e}^{-p}\,{dp}\\ \quad =\ \displaystyle \frac{k!}{{\left(\mathrm{ln}3\right)}^{k+1}},\end{array}\end{eqnarray} \tag{ 164 }$$

and obtain

$$\begin{eqnarray}&&\displaystyle \sum _{m=3k}^{\infty }\displaystyle \frac{{2}^{k+1}{m}^{k}}{{3}^{(m+2k+3)}k!}\lt \displaystyle \frac{{2}^{k+1}}{{3}^{(2k+3)}{\left(\mathrm{ln}3\right)}^{k+1}}.\end{eqnarray} \tag{ 165 }$$

This expression is also monotonically decreasing in k.

The head of the series in Equation (161) contains terms d_m(k), where $m\lt 3k$. We establish that these d_m(k) monotonically decrease in m. The ratio of successive terms in m is

$$\begin{eqnarray}&&\displaystyle \frac{{d}_{m+1}(k)}{{d}_{m}(k)}=\displaystyle \frac{1}{3}\displaystyle \prod _{j=0}^{k-1}\left(1+\displaystyle \frac{1}{m+k+2+2j}\right).\end{eqnarray} \tag{ 166 }$$

Each term is largest when the denominator is smallest. This is achieved when $j=m=0$, giving an upper bound for the product

$$\begin{eqnarray}&&\displaystyle \frac{{d}_{m+1}(k)}{{d}_{m}(k)}\leqslant \displaystyle \frac{1}{3}{\left(1+\displaystyle \frac{1}{k+2}\right)}^{k}\lt \displaystyle \frac{e}{3}\lt 1.\end{eqnarray} \tag{ 167 }$$

We conclude the terms monotonically decrease in m for fixed k. The multiplicands in the head d_m(k) are largest when m = 0, giving the following upper bound

$$\begin{eqnarray}&&\displaystyle \sum _{m=0}^{3k-1}{d}_{m}(k)\lt 6k\displaystyle \frac{(3k)!!}{{3}^{2k+3}k!k!!}.\end{eqnarray} \tag{ 168 }$$

This expression is also monotonically decreasing in k.

We may (finally) write a monotonically decreasing bound for c_k

$$\begin{eqnarray}&&{c}_{k}\lt \displaystyle \sum _{n=0}^{k+1}{c}_{n}(k)+\displaystyle \frac{{2}^{k+1}}{{3}^{(2k+3)}{\left(\mathrm{ln}3\right)}^{k+1}}+6k\displaystyle \frac{(3k)!!}{{3}^{2k+3}k!k!!}.\end{eqnarray} \tag{ 169 }$$

To show that all ${c}_{k}\lt 1$, it now suffices to explicitly compute the k = 0 bound. Substitution of k = 0 into Equation (169) gives

$$\begin{eqnarray}&&{c}_{0}\lt \left[{c}_{0}(0)+{c}_{1}(0)+\displaystyle \frac{2}{{3}^{3}\mathrm{ln}3}\right]\lt \displaystyle \frac{26}{27}\lt 1.\end{eqnarray} \tag{ 170 }$$

Because ${c}_{k}\geqslant 0$ always, this completes the proof that $\left|{G}_{{ss}}^{(k)}(0)\right|\,\leqslant 1$ as required by Assumption 2.

As a check on our work, central coefficients can be explicitly computed through parametric differentiation with respect to of Equation (145)

$$\begin{eqnarray}&&{c}_{k}=\mathop{\mathrm{lim}}\limits_{\epsilon \to 0}\displaystyle \frac{1}{k!}\displaystyle \frac{{\partial }^{k}}{\partial {\epsilon }^{k}}{ \mathcal P }(0).\end{eqnarray} \tag{ 171 }$$

For example, evaluation at k = 1 gives

$$\begin{eqnarray}&&{c}_{1}=\displaystyle \frac{1}{6},\end{eqnarray} \tag{ 172 }$$

which agrees with Equation (157) evaluated at k = 1. As an additional check, the bound has been verified numerically through c₆₀.

Now that we have established that the coefficients of the central pressure are bounded, we show that the coefficients for $r\gt 0$ remain bounded. This is done readily with series long division. We will show that each coefficient at $r\gt 0$ is decreased relative to its value at r = 0 and that coefficients remain bounded below by zero.

To prepare for the algorithm, note that Equation (145) can be written using uniformly convergent binomial expansions as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{ \mathcal P }(r)}{\rho }-1\\ =\,\left(1+\displaystyle \sum _{j=1}^{\infty }{a}_{j}{\epsilon }^{j}\right)/\left\{1+\displaystyle \frac{1}{2}\displaystyle \sum _{k=1}^{\infty }{a}_{k}{\epsilon }^{k}\left[3-{\left(\displaystyle \frac{r}{R}\right)}^{2k}\right]\right\},\end{array}\end{eqnarray} \tag{ 173 }$$

where the a_n are defined as

$$\begin{eqnarray}&&{a}_{n}\equiv {\left(-1\right)}^{n}{\left(\displaystyle \frac{2}{3}\right)}^{n}\left(\displaystyle \genfrac{}{}{0em}{}{1/2}{n}\right).\end{eqnarray} \tag{ 174 }$$

The long division algorithm states that, if

$$\begin{eqnarray}&&\displaystyle \sum _{k=0}^{\infty }{c}_{k}{\epsilon }^{k}\equiv \displaystyle \sum _{j=0}^{\infty }{a}_{j}{\epsilon }^{j}/\displaystyle \sum _{{\ell }=0}^{\infty }{b}_{{\ell }}{\epsilon }^{{\ell }},\end{eqnarray} \tag{ 175 }$$

then

$$\begin{eqnarray}&&{c}_{0}=\displaystyle \frac{{a}_{0}}{{b}_{0}}\end{eqnarray} \tag{ 176 }$$

$$\begin{eqnarray}&&{c}_{k}=\displaystyle \frac{1}{{b}_{0}}\left({a}_{k}-\displaystyle \sum _{j=0}^{k-1}{c}_{j}{b}_{k-j}\right)\qquad (k\geqslant 1).\end{eqnarray} \tag{ 177 }$$

By inspection of Equation (173), we see that

$$\begin{eqnarray}&&{a}_{0}={b}_{0}={c}_{0}=1.\end{eqnarray} \tag{ 178 }$$

Now, note that

$$\begin{eqnarray}&&\left(\displaystyle \genfrac{}{}{0em}{}{1/2}{j}\right)\propto {\left(-1\right)}^{j+1}\qquad (j\geqslant 1).\end{eqnarray} \tag{ 179 }$$

By inspection of Equation (173), we may now write

$$\begin{eqnarray}&&{c}_{k}(r)=-\left|{a}_{k}\right|+\displaystyle \frac{1}{2}\displaystyle \sum _{j=0}^{k-1}{c}_{j}(r)\left|{a}_{k-j}\right|\left[3-{\left(\displaystyle \frac{r}{R}\right)}^{2(k-j)}\right].\end{eqnarray} \tag{ 180 }$$

Pull off the first term of the sum and group the a_k terms

$$\begin{eqnarray}\begin{array}{rcl}{c}_{k}(r) & = & \left|{a}_{k}\right|\left\{\displaystyle \frac{1}{2}{c}_{0}(r)\left[3-{\left(\displaystyle \frac{r}{R}\right)}^{2k}\right]-1\right\}\\ & & +\,\displaystyle \frac{1}{2}\displaystyle \sum _{j=1}^{k-1}{c}_{j}(r)\left|{a}_{k-j}\right|\left[3-{\left(\displaystyle \frac{r}{R}\right)}^{2(k-j)}\right].\end{array}\end{eqnarray} \tag{ 181 }$$

It is now clear that ${c}_{k}(r)\geqslant 0\Rightarrow {c}_{k+1}(r)\geqslant 0$, provided that

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{c}_{0}(r)\left[3-{\left(\displaystyle \frac{r}{R}\right)}^{2k}\right]\geqslant 1.\end{eqnarray} \tag{ 182 }$$

Using that ${c}_{0}=1$, it follows that the all c_k(r) are positive if

$$\begin{eqnarray}&&{\left(\displaystyle \frac{r}{R}\right)}^{2k}\leqslant 1,\end{eqnarray} \tag{ 183 }$$

which is always true interior to the sphere. By inspection of Equation (181), we can now conclude that each c_k(r) decreases monotonically in r. Since we have already bounded the central coefficients, it follows that $\left|{G}_{{ss}}^{(k)}(r)\right|\lt 1$ for all k within the constant-density sphere, as was to be shown.