Implications of Symmetry and Pressure in Friedmann Cosmology. III. Point Sources of Dark Energy that Tend toward Uniformity

Given that GEODEs are compact objects formed from matter, it is natural to expect that their bulk dynamics would track those of cold dark matter (CDM). Specifically, one would expect GEODEs to "clump" with CDM and eventually become a subpopulation of MAssive Compact Halo Objects (MACHOs). If formed at z ∼ 15, Population III GEODEs would blueshift into 10³–10⁵ M_⊙ GEODEs today and would mimic intermediate-mass black holes (IMBHs), on which there are observational constraints.

The IMBH mass range has been previously investigated as candidate objects for CDM halos. For compact objects in the IMBH mass range, the most powerful constraint comes from the existence of binary star systems within the Milky Way. Binney & Tremaine (2008, (8.65b)) give the timescale for catastrophic disruption of binary star systems as

$$\begin{eqnarray}&&{t}_{d,\mathrm{cat}}\simeq \displaystyle \frac{0.07}{{\rho }_{\mathrm{dis}}}\sqrt{\displaystyle \frac{{M}_{{\rm{b}}}}{{{GA}}^{3}}},\end{eqnarray} \tag{ 3 }$$

where ${\rho }_{\mathrm{dis}}$ is the mass density of disruptors, M_b is the mass of the binary system, and A is the semimajor axis³ of the binary orbit. Systems are considered vulnerable to disruption if this timescale is less than the lifetime of the Milky Way, roughly 10 Gyr. Due to the scaling 1/A^3/2, wide binaries are more easily disrupted.

To develop constraints on MACHOs, one studies halo wide binaries: those binary pairs that have a center-of-mass orbit that takes the system far from the galactic plane. Monroy-Rodríguez & Allen (2014, Table 3) report observations of several such stellar pairs. Adopting the pair with the largest semimajor axis and requiring ${t}_{d,\mathrm{cat}}$ to exceed the age of the universe, we find

$$\begin{eqnarray}&&{\rho }_{\mathrm{dis}}\lesssim 8\times {10}^{-24}\,\mathrm{kg}\,{{\rm{m}}}^{-3}.\end{eqnarray} \tag{ 4 }$$

While this is 10³ larger than ρ_cr, if Population III GEODEs behave as CDM, it is natural to expect that their local densities far exceed background values. It becomes clear that any viable pointlike DE explanation of ρ_Λ must dynamically keep Population III GEODE densities low in regions of CDM.

In this paper, we address the MACHO problem by studying the dynamical behavior of Population III GEODEs. The structure of this paper is as follows. In Section 2, we provide definitions necessary to pose Einstein's equations perturbatively. In Section 3, we review the zero-order (i.e., Friedmann) cosmological source appropriate for a GEODE species. In Section 4, we determine the first-order stress perturbations appropriate for a GEODE species. In Section 5, we determine the dynamical equations that govern the evolution of these stress perturbations. In Section 6, we present results from numerical solutions of these dynamical equations. In Section 7, we interpret these results physically and discuss relations between Kerr BHs and GEODEs. In Section 8, we briefly conclude.

Throughout this work, unless otherwise indicated, we use the following units: densities are in the critical density today ρ_cr, times are in the reciprocal Hubble rate today ${H}_{0}^{-1}$, and Fourier wavenumber k is in ${\mathrm{Mpc}}^{-1}$. We make every attempt to use the same notation as Hu (2004). We adopt cosmological parameters from the TT,TE,EE+lowE+lensing+BAO Planck 2018 cosmology, given by Aghanim et al. (2018, Table 2). For consistency with structure formation literature, all densities will be physical (as opposed to comoving). We will often refer to redshifts and scale factors as "times." No confusion will arise because these quantities evolve monotonically in any realistic Friedmann model.

We adopt the following spatially flat RW metric representation:

$$\begin{eqnarray}&&{g}_{\mu \nu }(\eta ,{\boldsymbol{x}}):= a{\left(\eta \right)}^{2}\left[{\eta }_{\mu \nu }+\epsilon {h}_{\mu \nu }^{(1)}(\eta ,{\boldsymbol{x}})+\cdots \right]\qquad \epsilon \lt 1.\end{eqnarray} \tag{ 5 }$$

This is the definition given by Croker & Weiner (2019, Equation (14)), abbreviated⁴ at first order. Croker & Weiner (2019, Sections 3.2–3.4) show that the coordinates η and ${\boldsymbol{x}}$ continue to make sense inside many ultrarelativistic compact objects. We will assume that they are well defined everywhere, and we will use only these coordinates. The metric representation in Equation (5) is given with respect to a class of coordinate frame fields spanned by nonnormal, and possibly nonorthogonal, bases {∂_μ}. We will call these frame fields "RW frame fields."

RW frame fields are useful because they are spanned by coordinate bases. The explicit prescription of coordinate functions then allows explicit computation of representations. Across the differentiable manifold ${ \mathcal M }$, RW frames differ at most by order quantities. Taken together, RW frames can be constructed to give representations that closely approximate nonrelativistic mechanics on large regions of spacetime. Note that first-order equations of motion are produced from variation of the ${h}_{\mu v}^{(1)}$. This means that source terms to Einstein's equations for the metric perturbation degrees of freedom must be tensor components in some RW frame field.

The stress tensor ${T}_{\mu \nu }$ is, by definition, proportional to the variation of the nongravitational Lagrange density with respect to the true (i.e., unapproximated) inverse metric. Unfortunately, we do not know the true metric, only that it uniquely exists. We may, however, leverage this uniqueness to give an invariant description of the true stress tensor.

The tensor T_μν can be regarded as a linear map from each tangent space at a point P to its dual space at P,

$$\begin{eqnarray}&&{T}_{\mu \nu }:{T}_{P}{ \mathcal M }\to {T}_{P}^{* }{ \mathcal M }.\end{eqnarray} \tag{ 6 }$$

In other words, "T_μν acting on an index up quantity produces an index down quantity." If we contract with the inverse metric at P, we produce a linear operator T on ${T}_{P}{ \mathcal M }$,

$$\begin{eqnarray}&&T:= {g}^{\mu \alpha }{T}_{\alpha \nu }:{T}_{P}{ \mathcal M }\to {T}_{P}{ \mathcal M }.\end{eqnarray} \tag{ 7 }$$

The operator T decomposes into the sum of a diagonalizable operator T_I and a nilpotent⁵ operator T_II. The diagonalizable operator T_I describes timelike stress, while the nilpotent operator T_II describes lightlike stress.

At any point P, T_I has one timelike and three spacelike eigenvectors. If the type I stress at P has no particular symmetries, as expected in general, these eigenvectors define a unique frame field $\left\{{\overline{e}}_{\mu }\right\}$ (up to permutations in order and sign). If there are symmetries at P, then there are repeat eigenvalues and many such frames at P. We will call these frame fields "eigenframe fields." We will distinguish quantities in an eigenframe field with an overline.

Eigenframe fields are useful because they give an invariant characterization of the stress tensor. In any eigenframe field, the representation of T_I is diagonal:

$$\begin{eqnarray}&&\bar{{T}_{\mu \nu }^{{\rm{I}}}}={\rm{diag}}\left[\bar{\rho },\bar{p[1]},\bar{p[2]},\bar{p[3]}\right],\end{eqnarray} \tag{ 8 }$$

where the energy density $\overline{\rho }$ and principal pressures $\overline{p[{\ell }]}$ are the eigenvalues of T_I at P. In other words, the components of the stress tensor in an eigenframe are invariant, physical quantities. Suppose the timelike eigenvalue is unique. Then the timelike eigenvector ${\overline{e}}_{0}$ becomes the velocity of the type I stress at P, if viewed from any noneigenframe at P. From the perspective of eigenframe fields, material that can be at rest is at rest.

Astrophysical sources will contribute energy densities and principal pressures that greatly exceed the background averages. More so, these large deviations from the background only occur at the positions of the sources. Typically then, $\delta {\rho }^{{\boldsymbol{x}}}$, $\delta {p}^{{\boldsymbol{x}}}$, and ${\left({{\rm{\Pi }}}_{k}^{j}\right)}^{{\boldsymbol{x}}}$ are regarded as continuum fluid approximations. In the continuum fluid approximation (e.g., Landau & Lifshitz 1959, Section 1), one coarse-grains the source into "fluid elements." The fluid elements must be large enough so that a bulk description is useful on the desired scales. The fluid elements must also be small enough so that four-momentum is distributed within the elements faster than between them. An average is then performed over each fluid element volume to produce density and pressure fields.

Unfortunately, no quantitative prescription is given for any part of this procedure. In other words, the ambiguity resolved by Croker & Weiner (2019) at zero order reemerges in the first-order theory. A spatial average must be taken, but no length scale over which to average, or Jacobian to use during integration, is unambiguously prescribed.

In Fourier space, however, this pathology is bypassed. In practice, the Fourier expansion is always truncated at some finite k_max. Upon inversion of the transform, this leads to first-order real-space sources that are low-pass filtered below $2\pi /{k}_{\max }$. All of the resulting dynamical variables now smoothly track the distribution in real space of the pointlike objects. In fact, there is a measured scale 2π/k_nl (Dodelson 2003, Section 7.1.1), above which the observed matter power spectrum agrees with linear theory predictions. This scale defines a continuum fluid approximation. It is the first-order analog to the zero-order homogeneity scale ${{ \mathcal V }}^{1/3}$. In summary,

• 1.  
  a well-defined, first-order theory requires the Fourier transforms of Equations (19), not the Fourier transforms of the continuum fluid approximations to Equations (19).

In particular, pressures interior to all objects will get "smeared" into all modes below k_max. For ensembles of ultrarelativistic objects, peculiar flows at first order can then become influenced by these objects' internal dynamics. This is the same phenomenon observed at zero order, but replayed at first order.

We have used Fourier transforms to produce a well-defined continuum fluid approximation from the actual first-order sources to Einstein's equations. D'Eath (1976, Section 3) proves that solutions to Einstein's equations for perturbed RW spaces with spatial curvature ${{\rm{\Omega }}}_{K}:= 0$ do not generally exist unless perturbations decay sufficiently fast at infinity. This pathology implies that the Fourier transforms of Equation (19) may diverge.

The standard approach (e.g., Peebles 1980, Section 27) regards the universe as "box periodic." This procedure changes the global topology of the universe to that of a three-torus. We adopt this convention and define spatial Fourier transforms on a three-torus with radii much larger than the zero-order homogeneity scale ${{ \mathcal V }}^{1/3}$.

In this section, we define the GEODE ensemble in terms of a local model and use it to compute the background equation of state −1 + χ. In order that our definition be invariant, we must work in an eigenframe field. Regard $\overline{{\rho }_{G}^{{\boldsymbol{x}}}}$ and $\overline{{p}_{G}^{{\boldsymbol{x}}}}$ as being sourced by a population of N GEODEs at positions ${{\boldsymbol{x}}}_{j}(\eta )$:

$$\begin{eqnarray}\begin{array}{rcl}\overline{{\rho }_{G}^{{\boldsymbol{x}}}} & := & \displaystyle \sum _{j}^{N}\left[u(\eta )+\overline{{\kappa }_{0}}(\eta ,{\boldsymbol{x}}-{{\boldsymbol{x}}}_{j})\right]{{\bf{1}}}_{{{ \mathcal B }}_{j}}\\ \overline{{p}_{G}^{{\boldsymbol{x}}}} & := & \displaystyle \sum _{j}^{N}\left[-u(\eta )+\overline{{\kappa }_{1}}(\eta ,{\boldsymbol{x}}-{{\boldsymbol{x}}}_{j})\right]{{\bf{1}}}_{{{ \mathcal B }}_{j}}.\end{array}\end{eqnarray} \tag{ 22 }$$

Each GEODE is defined on a support of radius⁶ R centered at ${{\boldsymbol{x}}}_{j}(\eta )$,

$$\begin{eqnarray}{{\bf{1}}}_{{{ \mathcal B }}_{j}}({\boldsymbol{x}}):= \left\{\begin{array}{ll}1 & \left|{\boldsymbol{x}}-{{\boldsymbol{x}}}_{j}(\eta )\right|\leqslant R\\ 0 & \left|{\boldsymbol{x}}-{{\boldsymbol{x}}}_{j}(\eta )\right|\gt R,\end{array}\right.\end{eqnarray} \tag{ 23 }$$

and separated into a de Sitter portion with uniform energy density u and a superposed "crust." The crust is characterized by an energy density $\overline{{\kappa }_{0}}$ and an isotropic pressure $\overline{{\kappa }_{1}}$.

From Equation (16), it follows that

$$\begin{eqnarray}&&\displaystyle \frac{1}{{ \mathcal V }}\displaystyle \sum _{j}^{N}{\int }_{{{ \mathcal B }}_{j}}\left[(-1+\chi )(u+\overline{{\kappa }_{0}})+u-\overline{{\kappa }_{1}}\right]\,{d}^{3}x=0.\end{eqnarray} \tag{ 24 }$$

All GEODEs are presumed identical, so

$$\begin{eqnarray}&&\displaystyle \frac{N{ \mathcal B }}{{ \mathcal V }}\left[(-1+\chi )(u+{\left\langle \overline{{\kappa }_{0}}\right\rangle }_{{ \mathcal B }})+u-{\left\langle \overline{{\kappa }_{1}}\right\rangle }_{{ \mathcal B }}\right]=0.\end{eqnarray} \tag{ 25 }$$

Rearranging, and assuming that ${\left\langle \overline{{\kappa }_{0}}\right\rangle }_{{ \mathcal B }}\ll u$, the binomial expansion gives

$$\begin{eqnarray}&&\chi =\displaystyle \frac{{\left\langle \overline{{\kappa }_{1}}+\overline{{\kappa }_{0}}\right\rangle }_{{ \mathcal B }}}{u}+\cdots .\end{eqnarray} \tag{ 26 }$$

At leading order, χ is the GEODE volume-averaged sum of the crust energy density and pressure, scaled by the GEODE interior DE density.

In Fourier space, Hu (2004, Equation (106)) defines the squared nonadiabatic sound speed ${c}_{G}^{2}$ as

$$\begin{eqnarray}&&{c}_{G}^{2}({\boldsymbol{k}},\eta ):=\displaystyle \frac{\delta {p}_{G}({\boldsymbol{k}},\eta )}{\delta {\rho }_{G}({\boldsymbol{k}},\eta )}{| }_{{\rm{GEODE}}{\rm{rest}}}.\end{eqnarray} \tag{ 27 }$$

For point objects, the "rest" gauge for a particular constituent comoves⁷ with the peculiar flow of that constituent. GEODEs are extended objects, so the "rest" gauge perceives mass centers to be at rest. These RW frame fields are not, however, fully adapted to individual GEODEs. They will still perceive intrinsic GEODE motions, such as spin. For this reason, we do not use the word "rest." Instead, we will refer to this gauge as the GEODE "flow gauge."

From Equations (19) and (22), the Fourier transforms $\delta {p}_{G}$ and $\delta {\rho }_{G}$ will track the GEODE distribution in space at modes below ∼1/R, due to the presence of many ${{\bf{1}}}_{{{ \mathcal B }}_{j}}$. This remains true in the flow gauge, so ${c}_{G}^{2}$ is constant in space for modes below ∼1/R. We use this to compute ${c}_{G}^{2}$ directly from the Fourier transforms. For k ≪ 2π/R, we have that

$$\begin{eqnarray}&&\int {e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}\left({p}_{G}^{{\boldsymbol{x}}}-{p}_{G}\right)\,{d}^{3}x={c}_{G}^{2}\int {e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}\left({\rho }_{G}^{{\boldsymbol{x}}}-{\rho }_{G}\right)\,{d}^{3}x.\end{eqnarray} \tag{ 28 }$$

At ${\boldsymbol{k}}=0$, the computation is just a spatial average and the fluctuations vanish, independently of ${c}_{G}^{2}$. Consider ${\boldsymbol{k}}\ne 0$. By orthogonality, all spatially independent contributions vanish:

$$\begin{eqnarray}&&\int {e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}{p}_{G}^{{\boldsymbol{x}}}\,{d}^{3}x={c}_{G}^{2}\int {e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}{\rho }_{G}^{{\boldsymbol{x}}}\,{d}^{3}x.\end{eqnarray} \tag{ 29 }$$

Substitution of Equations (22), formally transformed into the flow gauge, gives

$$\begin{eqnarray}&&\int {e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}\displaystyle \sum _{j}^{N}\left[{\kappa }_{1}-u-{c}_{G}^{2}({\kappa }_{0}+u)\right]{{\bf{1}}}_{{{ \mathcal B }}_{j}}\,{d}^{3}x=0.\end{eqnarray} \tag{ 30 }$$

The pure de Sitter contributions are invariant because pure DE is a multiple of the metric. We may use the supports of the individual GEODEs to reexpress the transform as a sum over the ensemble

$$\begin{eqnarray}&&\displaystyle \sum _{j}^{N}{\int }_{{{ \mathcal B }}_{j}}{e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}}\left[{\kappa }_{1}-u-{c}_{G}^{2}\left({\kappa }_{0}+u\right)\right]\,{d}^{3}x=0.\end{eqnarray} \tag{ 31 }$$

From each term, factor out the constant phase $\exp (i{\boldsymbol{k}}\cdot {{\boldsymbol{x}}}_{j})$:

$$\begin{eqnarray}&&\displaystyle \sum _{j}^{N}{e}^{i{\boldsymbol{k}}\cdot {{\boldsymbol{x}}}_{j}}{\int }_{{{ \mathcal B }}_{j}}{e}^{i{\boldsymbol{k}}\cdot ({\boldsymbol{x}}-{{\boldsymbol{x}}}_{j})}\left[{\kappa }_{1}-u-{c}_{G}^{2}\left({\kappa }_{0}+u\right)\right]\,{d}^{3}x=0.\end{eqnarray} \tag{ 32 }$$

On each support, $| {\boldsymbol{k}}| \ll 1/R$, so the exponential under the integral is essentially 1,

$$\begin{eqnarray}&&\displaystyle \sum _{j}^{N}{e}^{i{\boldsymbol{k}}\cdot {{\boldsymbol{x}}}_{j}}{\int }_{{{ \mathcal B }}_{j}}{\kappa }_{1}-u-{c}_{G}^{2}\left({\kappa }_{0}+u\right)\,{d}^{3}x=0.\end{eqnarray} \tag{ 33 }$$

Performing the integrations gives

$$\begin{eqnarray}&&{ \mathcal B }\displaystyle \sum _{j}^{N}{e}^{i{\boldsymbol{k}}\cdot {{\boldsymbol{x}}}_{j}}\left[{\left\langle {\kappa }_{1}\right\rangle }_{{{ \mathcal B }}_{j}}-u-{c}_{G}^{2}\left({\left\langle {\kappa }_{0}\right\rangle }_{{{ \mathcal B }}_{j}}+u\right)\right]=0.\end{eqnarray} \tag{ 34 }$$

Because all GEODEs are presumed identical, the term in square braces can be factored out:

$$\begin{eqnarray}&&\left[{\left\langle {\kappa }_{1}\right\rangle }_{{ \mathcal B }}-u-{c}_{G}^{2}\left({\left\langle {\kappa }_{0}\right\rangle }_{{ \mathcal B }}+u\right)\right]\displaystyle \sum _{j}^{N}{e}^{i{\boldsymbol{k}}\cdot {{\boldsymbol{x}}}_{j}}=0.\end{eqnarray} \tag{ 35 }$$

If the sum happens to vanish at the fixed ${\boldsymbol{k}}$ we have chosen, no constraint on ${c}_{G}^{2}$ can be established. In Appendix A, we prove that the sum of phases is nonzero almost everywhere in ${\boldsymbol{k}}$ for any distribution of N GEODEs in space. So if the sum vanishes, there always exists $\tau \in {{\mathbb{R}}}^{3}$ with $| \tau |$ arbitrarily small so that the sum does not vanish at ${\boldsymbol{k}}+\tau$. In this sense, we may always divide by the sum to arrive at the desired relation,

$$\begin{eqnarray}&&{c}_{G}^{2}=\displaystyle \frac{-u+{\left\langle {\kappa }_{1}\right\rangle }_{{ \mathcal B }}}{u+{\left\langle {\kappa }_{0}\right\rangle }_{{ \mathcal B }}}.\end{eqnarray} \tag{ 36 }$$

In order to express this parameter in terms of intrinsic GEODE properties, we must express the flow gauge stress components in terms of eigenvalues. Pick any GEODE. The flow frames near this GEODE perceive only higher-order center-of-mass motion. Pick a flow frame at some point P in the GEODE crust. Suppose that this flow frame perceives the crust in motion, spinning at velocity $\beta (\eta ,{\boldsymbol{x}})$ in some direction denoted $\hat{\phi }$. The transformation from an eigenframe at P to the flow frame at P is detailed in Appendix B. After transformation, the crust energy density becomes

$$\begin{eqnarray}&&{\kappa }_{0}={\gamma }^{2}\left(\overline{{\kappa }_{0}}+{\beta }^{2}\overline{{\kappa }_{1}}\right)+O(\epsilon ),\end{eqnarray} \tag{ 37 }$$

while the crust principal pressures become

$$\begin{eqnarray}&&{p}_{G}^{{\boldsymbol{x}}}[\hat{\phi }]={\gamma }^{2}\left(\overline{{\kappa }_{1}}+{\beta }^{2}\overline{{\kappa }_{0}}\right)+O(\epsilon )\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{p}_{G}^{{\boldsymbol{x}}}[\hat{r}]=\overline{{\kappa }_{1}}+O(\epsilon )\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{p}_{G}^{{\boldsymbol{x}}}[\hat{\theta }]=\overline{{\kappa }_{1}}+O(\epsilon ).\end{eqnarray} \tag{ 40 }$$

The frames are related by a boost, with order corrections. Because ${c}_{G}^{2}$ relates two order- quantities, these corrections are to be dropped. By definition, κ₁ is the isotropic pressure in the flow gauge, so

$$\begin{eqnarray}&&{\kappa }_{1}=\displaystyle \frac{1}{3}\left(2\overline{{\kappa }_{1}}+\displaystyle \frac{\overline{{\kappa }_{1}}+{\beta }^{2}\overline{{\kappa }_{0}}}{1-{\beta }^{2}}\right).\end{eqnarray} \tag{ 41 }$$

We may now express ${c}_{G}^{2}$ in terms of GEODE intrinsic properties. If the GEODE crust rotates rigidly, then functions of β(η) have no position dependence and commute through spatial averages. Substituting Equations (37) and (41) into Equation (36) gives

$$\begin{eqnarray}&&{c}_{G}^{2}=\displaystyle \frac{1}{3}\left[\displaystyle \frac{2{\left\langle \overline{{\kappa }_{1}}\right\rangle }_{{ \mathcal B }}+{\gamma }^{2}\left({\left\langle \overline{{\kappa }_{1}}\right\rangle }_{{ \mathcal B }}+{\beta }^{2}{\left\langle \overline{{\kappa }_{0}}\right\rangle }_{{ \mathcal B }}\right)-3u}{{\gamma }^{2}\left({\left\langle \overline{{\kappa }_{0}}\right\rangle }_{{ \mathcal B }}+{\beta }^{2}{\left\langle \overline{{\kappa }_{1}}\right\rangle }_{{ \mathcal B }}\right)+u}\right].\end{eqnarray} \tag{ 42 }$$

It is instructive to consider the two limits

$$\begin{eqnarray}{c}_{G}^{2}=\left\{\begin{array}{ll}-1+\chi +... & \beta \to 0\\ 1/3 & \beta \to 1.\end{array}\right.\end{eqnarray} \tag{ 43 }$$

The nonspinning β → 0 limit sets the squared sound speed of density perturbations equal to the background equation of state. As Hu (2004, Section 5.2) remarks, a DE model with this behavior will undergo accelerated clumping. This is consistent with expectations from Croker et al. (2020, Section 3). Pointlike GEODEs at small scales, without spin, aggregate more quickly due to blueshifting mass.

The rapidly spinning β → 1 limit can be recognized as perceiving the crust as a nearly null flux, with equation of state 1/3. Like radiation overdensities, a GEODE population with very high spin will tend to disperse at first order.

The dynamics for the velocity field of any decoupled constituent are given in an arbitrary gauge by Hu (2004, Equation (19)),

$$\begin{eqnarray}&&\begin{array}{l}\left[\displaystyle \frac{d}{d\eta }+4\displaystyle \frac{\dot{a}}{a}\right]({\rho }_{G}+{p}_{G})\displaystyle \frac{{v}_{G}-B}{k}\\ \quad =\,\delta {p}_{G}-\displaystyle \frac{2}{3}{p}_{G}{{\rm{\Pi }}}_{G}+({\rho }_{G}+{p}_{G}){\rm{\Psi }}.\end{array}\end{eqnarray} \tag{ 46 }$$

Here B is the vector shift, which is defined to be 0 in Newtonian gauge, and we have set Hu's RW spatial curvature K = 0. Substituting the background equation of state from Equation (16) gives

$$\begin{eqnarray}&&\begin{array}{l}\left[\displaystyle \frac{d}{d\eta }+4\displaystyle \frac{\dot{a}}{a}\right]\chi {\rho }_{G}\displaystyle \frac{{v}_{G}-B}{k}\\ \quad =\,\delta {p}_{G}+\displaystyle \frac{2}{3}(1-\chi ){\rho }_{G}{{\rm{\Pi }}}_{G}+\chi {\rho }_{G}{\rm{\Psi }}.\end{array}\end{eqnarray} \tag{ 47 }$$

A dynamical model for the peculiar velocity field requires the derivative term ${\dot{v}}_{G}$, so we divide by χ:

$$\begin{eqnarray}&&\begin{array}{l}\left[\displaystyle \frac{d}{d\eta }+4\displaystyle \frac{\dot{a}}{a}\right]{\rho }_{G}\displaystyle \frac{{v}_{G}-B}{k}={\chi }^{-1}\delta {p}_{G}\\ \quad +\,\displaystyle \frac{2}{3}({\chi }^{-1}-1){\rho }_{G}{{\rm{\Pi }}}_{G}+{\rho }_{G}{\rm{\Psi }}.\end{array}\end{eqnarray} \tag{ 48 }$$

To model a population of individual DE point objects, however, we need to take χ → 0. Evidently, the dynamics cannot be perturbatively consistent in the χ → 0 limit unless

$$\begin{eqnarray}&&\delta {p}_{G}=-\displaystyle \frac{2}{3}{\rho }_{G}{{\rm{\Pi }}}_{G}.\end{eqnarray} \tag{ 49 }$$

An algebraic relation between δp_G and Π_G is unsurprising, given Equation (27) and that k ≪ 1/R. It is remarkable, however, that Equation (49) contains no explicit ${\boldsymbol{k}}$ dependence. The proportionality holds at all modes, not just modes below ∼1/R. In other words, viable GEODE interiors must have nonvanishing anisotropic stress. This result is consistent with the conclusions of Cattoen et al. (2005), who found that a wide class of GEODEs must feature anisotropic stress. For k ≪ 1/R, the GEODE anisotropic stress perturbation is determined by combining Equation (49) with Hu (2004, Equation (107)),

$$\begin{eqnarray}&&{{\rm{\Pi }}}_{G}=-\displaystyle \frac{3}{2}\left[{c}_{G}^{2}{\delta }_{G}+3\chi {Ha}\displaystyle \frac{{v}_{G}-B}{{{kk}}_{\mathrm{fac}}}({c}_{G}^{2}+1-\chi )\right].\end{eqnarray} \tag{ 50 }$$

We are now in position to define the dynamics of the GEODE peculiar velocity field. Substitution of Equation (49) back into Equation (46) gives

$$\begin{eqnarray}&&\displaystyle \frac{\chi {\rho }_{G}}{k}\left[{\dot{v}}_{G}=\dot{B}+(3\chi -4)({v}_{G}-B)\displaystyle \frac{\dot{a}}{a}+k{\rm{\Psi }}-\displaystyle \frac{2}{3}k{{\rm{\Pi }}}_{G}\right].\end{eqnarray} \tag{ 51 }$$

All quantities outside the square braces are nonzero, and we divide them off:

$$\begin{eqnarray}&&{\dot{v}}_{G}=\dot{B}+(3\chi -4)({v}_{G}-B)\displaystyle \frac{\dot{a}}{a}+k{\rm{\Psi }}-\displaystyle \frac{2}{3}k{{\rm{\Pi }}}_{G}.\end{eqnarray} \tag{ 52 }$$

This equation is true for an arbitrary gauge, so the dynamics of the Population III GEODE velocity field are well defined for arbitrarily small χ, provided that Equation (49) is satisfied. Fixing to Newtonian gauge, expressing derivatives in terms of the scale factor, scaling units, and rearranging give

$$\begin{eqnarray}&&\displaystyle \frac{{{dv}}_{G}}{{da}}=(3\chi -4)\displaystyle \frac{{v}_{G}}{a}+\displaystyle \frac{{{kk}}_{\mathrm{fac}}}{{{Ha}}^{2}}\,\left({\rm{\Psi }}-\displaystyle \frac{2}{3}{{\rm{\Pi }}}_{G}\right).\end{eqnarray} \tag{ 53 }$$

Note that GEODEs with χ → 0 have the largest possible "Hubble drag" consistent with the DEC.

The model we have constructed has three constant parameters: the redshift of burst formation z_G, the deviation from perfect de Sitter interior χ > 0, and the nonadiabatic sound speed ${c}_{G}^{2}$. The time of burst formation z_G affects the fraction of baryons Ξ required to produce the observed Ω_Λ. Planck determines a baryon density parameter at recombination:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{b}{h}^{2}=(2.2236\pm 0.015)\times {10}^{-2}\quad (68 \% \ \mathrm{confidence}).\end{eqnarray} \tag{ 54 }$$

GEODEs, however, form at much later times. Independent late-time measurements of Ω_b by Macquart et al. (2020) with fast radio bursts are consistent with the measured value at recombination, with 40% uncertainties.

We assume that the epoch of Population III star formation is 8 ≤ z ≤ 20. In our cosmology, this occurs over ∼0.5 Gyr. Inoue et al. (2014, Figure 2) use the extragalactic background light (EBL) to place constraints on the Population III star formation rate (SFR). Based on this figure, an averaged SFR density for z ≤ 20 is $\sim 5\times {10}^{-2}{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}$. Combining these estimates gives the following density in Population III stars,

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{III}}\sim 2\times {10}^{-4},\end{eqnarray} \tag{ 55 }$$

expressed in units of the critical density today. No Population III star has been observed, so we regard the entirety of this fraction as collapsing and accreting onto GEODEs. Because the total baryon fraction is Ω_b ∼ 0.05, this means

$$\begin{eqnarray}&&{\rm{\Xi }}\sim 4\times {10}^{-3}.\end{eqnarray} \tag{ 56 }$$

We combine this collapse fraction with the measured χ from Aghanim et al. (2018, Equation (51)),

$$\begin{eqnarray}&&\chi \leqslant 0.05\qquad (95 \% \ \mathrm{confidence}),\end{eqnarray} \tag{ 57 }$$

and use Equation (18) to estimate an approximate burst time,

$$\begin{eqnarray}{z}_{G}\sim \left\{\begin{array}{ll}15 & \chi \to 0.05\\ 14 & \chi \to 0.\end{array}\right.\end{eqnarray} \tag{ 58 }$$

As expected, this value lies in the middle of our assumed Population III formation range.

We have ignored the details of accretion and dynamic spin-up in our simplified treatment. The approximate burst formation time given in Equation (58) is consistent with EBL constraints and regenerates the correct Ω_Λ. It need not be appropriate for inferring GEODE distribution in space. This is because GEODEs need not initially form with high spin and so may initially clump. There may be a delay between when GEODEs form and when GEODEs reach sufficient spin to dynamically repel. For dynamical studies, we will consider burst formation times over the entire Population III range,

$$\begin{eqnarray}&&8\leqslant {z}_{G}\leqslant 20.\end{eqnarray} \tag{ 59 }$$

In order that the GEODE population resist collapse, ${c}_{G}^{2}\gt 0$. From Equation (42), this requirement becomes

$$\begin{eqnarray}&&2{\left\langle \overline{{\kappa }_{1}}\right\rangle }_{{ \mathcal B }}+{\gamma }^{2}\left({\left\langle \overline{{\kappa }_{1}}\right\rangle }_{{ \mathcal B }}+{\beta }^{2}{\left\langle \overline{{\kappa }_{0}}\right\rangle }_{{ \mathcal B }}\right)\gt 3u.\end{eqnarray} \tag{ 60 }$$

Essentially, the boosted pressure within the crust must exceed the energy density of the de Sitter region by a factor of 3. We can develop a lower bound by assuming a dusty crust $\left\langle \overline{{\kappa }_{1}}\right\rangle := 0$,

$$\begin{eqnarray}&&{\gamma }^{2}{\beta }^{2}\chi \gt 3.\end{eqnarray} \tag{ 61 }$$

Evaluation with Equation (57) gives

$$\begin{eqnarray}&&\gamma \gtrsim 8.\end{eqnarray} \tag{ 62 }$$

In this regime, we may determine a constraint between ${c}_{G}^{2}$, χ, and γ:

$$\begin{eqnarray}&&{c}_{G}^{2}\simeq \displaystyle \frac{1}{3}-\displaystyle \frac{4}{3{\gamma }^{2}\chi }.\end{eqnarray} \tag{ 63 }$$

A discussion of plausible γ and its relation to the Kerr solution can be found in Section 7.2. We will consider the following range of nonadiabatic squared sound speeds,

$$\begin{eqnarray}&&0\lt {c}_{G}^{2}\lt \displaystyle \frac{1}{3}.\end{eqnarray} \tag{ 64 }$$

We solved the initial value problem presented in Section 5, Appendices C, and D with the odeint package in Python 3.5. Initial conditions were determined using the CLASS Einstein–Boltzmann code (Blas et al. 2011). We verified the correct configuration of CLASS, and our analysis tool chain, by regenerating Planck's reported ${\sigma }_{8}$ to <0.5σ from directly output Fourier mode amplitudes. We verified our Einstein code in two ways under ΛCDM assumptions. We regenerated an indistinguishable σ₈ by integrating from z_G with initial conditions generated by CLASS. We also regenerated σ₈ to ∼5% by integrating from primordial times, regarding baryons as CDM. In a Population III GEODE scenario, Ω_Λ = 0 until the formation of GEODEs at z_G. Despite this, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$ set to the typical value at primordial times within the CLASS code can have only a negligible effect on all initial conditions used in this study.

We consider two burst formation scenarios bracketing our assumed epoch of Population III star formation: z_G = 20 and z_G = 8. Figure 1 displays an early-burst formation at z_G = 20 for parameter values at log spacing. Initially, GEODEs start with spectra identical to baryons. For all parameter values, these spectra damp substantially during evolution. GEODE linear power for z < 2 is always suppressed by at least 10⁵, relative to the peak linear CDM power ∼10³. Upon approach to the present day, depending on parameter space, GEODE power can become further suppressed by 10³. In other words, early-burst GEODEs with γ ≳ 8 tend toward uniformity. They achieve near uniformity before the epoch of galaxy formation, and density contrast continues to damp during galaxy evolution.

Zoom In Zoom Out Reset image size

Figure 2 displays a late-burst formation at z_G = 8 for parameter values at log spacing. We regard this scenario as a dynamical proxy to an earlier burst formation but with a spin-up delay. Initially, GEODEs start with spectra identical to baryons. Again, for all parameter values, the spectra damp substantially during evolution. GEODE linear power for z < 2 is always suppressed by at least 10³, relative to peak linear CDM power ∼10³. Due to a later burst, oscillations about uniformity can persist through the epoch of galaxy formation until the present day. In other words, late-burst GEODEs with γ ≳ 8 also tend toward uniformity but can exhibit residual structure.

Zoom In Zoom Out Reset image size

The first-order analysis we have performed cannot make strong statements about the distribution of Population III GEODEs above the nonlinearity scale k_nl. For a late-burst scenario, possible structure can remain in the GEODE distribution. It is then of interest to ask under what circumstances do GEODEs become anticorrelated at the smallest scales accessible to linear theory. Consider anticorrelation achieved around the epoch of galaxy formation that persists until the present day. In such scenarios, it becomes even less likely that the nonlinear evolution of the GEODE density contrast would build significant overdensities on scales ≳k_nl.

To investigate this question, we compute the GEODE-CDM cross-correlation coefficient:

$$\begin{eqnarray}&&\xi (z):= \displaystyle \frac{{\displaystyle \sum }_{k}{\delta }_{G}(z,k){\delta }_{\mathrm{CDM}}(z,k)}{\sqrt{\left({\displaystyle \sum }_{k}{\delta }_{G}{\left(z,k\right)}^{2}\right)\left({\displaystyle \sum }_{k}{\delta }_{\mathrm{CDM}}{\left(z,k\right)}^{2}\right)}}.\end{eqnarray} \tag{ 65 }$$

We restrict our attention to the high end of the linear theory $0.01\leqslant k\leqslant {k}_{\mathrm{nl}}$. In Figure 3, we display smoothed values of ξ(0) determined from 400 uniform samples of parameter space for a burst formation at z_G = 8. Much of the parameter space produces GEODE populations that are weakly correlated or anticorrelated with CDM at the present day. Overlaid on this figure is also the requisite spin boost γ required to achieve the combination of χ and ${c}_{G}^{2}$. In Figure 4, we magnify the anticorrelated region of Figure 3. Overplotted contours represent regions of parameter space where the GEODE density contrast has been anticorrelated with CDM on average since z < 2. Evidently, regions of parameter space exist where GEODEs have been separate from CDM for ∼10 Gyr and remain so today.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All dynamics can be understood as damped oscillations about a uniformly distributed equilibrium. From Equation (27), the overpressure is ${c}_{G}^{2}$ times the density contrast. An appreciable and positive ${c}_{G}^{2}$ will cause Population III GEODEs to disperse from their progenitor baryonic overdensities. If they become underdense, gravitational attraction with CDM brings them back together. A large ${c}_{G}^{2}$ means greater pressures, and so larger amplitude oscillation and longer time to damp. Inspection of Equation (C4) reveals the effect of χ. A large χ implies less Hubble friction and greater sensitivity to the gravitational environment, and so more prone to oscillation. A large z_G implies earlier relaxation to low GEODE power. At small z_G, matter overdensities have had more time to grow, so GEODE contrast will be more prone to oscillation.

MACHO constraints in the IMBH range are consistent with a density of disruptors 10³ρ_cr. We have shown that high-spin Population III GEODEs disperse and tend toward a uniform distribution with density ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$. In other words, Population III GEODEs are consistent with MACHO constraints over large ranges of parameter space.

In Section 4, we found that GEODE peculiar flow has two regimes, depending on spin. For γ ≲ 8, GEODEs undergo accelerated clustering. For γ ≳ 8, GEODEs effectively repel each other. To place this value for γ in context, contemporary terrestrial accelerators, such as the Large Hadron Collider (LHC; e.g., Evans & Bryant 2008), readily achieve proton beams with γ > 10³. The radius of the LHC proton ring is ∼4.3 km, which is roughly the Schwarzschild radius of a 1 M_⊙ BH. Astrophysically, however, Hessels et al. (2006) report the fastest spinning pulsar has surface γ = 1.03.

Because GEODEs mimic BHs, it is of interest to relate γ to the Kerr dimensionless spin parameter a_*. A reasonable definition of tangential velocity for a Kerr hole is somewhat arcane. Smarr (1973, Equation (20)) gives a definition that may be expressed in terms of the Kerr dimensionless spin via Smarr (1973, Equations (11)). Converting to contemporary⁹ notation, this relation becomes

$$\begin{eqnarray}&&\beta =\displaystyle \frac{\sqrt{1-\sqrt{1-{a}_{* }^{2}}}}{\sqrt{1+\sqrt{1-{a}_{* }^{2}}}}\qquad \gamma =\displaystyle \frac{1}{\sqrt{2}}\sqrt{1+\displaystyle \frac{1}{\sqrt{1-{a}_{* }^{2}}}}.\end{eqnarray} \tag{ 66 }$$

In order to achieve γ ≳ 8, we require a Kerr dimensionless spin

$$\begin{eqnarray}&&{a}_{* }\gt 0.99997.\end{eqnarray} \tag{ 67 }$$

The largest Kerr spin parameter measured by the LIGO Scientific Collaboration et al. (2019, Table III) is GW 170729, with 0.68 < a_* < 0.88, which becomes γ < 2.2. Walton et al. (2016) report 0.93 ≲ a_* ≲ 0.96 for the accreting X-ray binary BH Cygnus X-1, which becomes γ ≲ 3.2. If these objects are also GEODEs, these data suggest that GEODEs do not initially have sufficient spin to enter the repulsive regime.

Repulsive behavior is required for GEODEs to be a viable late-time pointlike DE candidate. While this behavior is also superficially required for stability of the dynamical equations, β may be time dependent within the model we have considered. Assume a conservative γ(z_G) ∼ 2, based on spin parameters reported by the LIGO Scientific Collaboration et al. (2019). In this regime, GEODE peculiar flow will undergo accelerated collapse.

Population III GEODEs form at an epoch with much higher density and from much larger progenitor stars. Given angular-momentum-conserving accretion due to GEODE blueshift and frequent mergers while in the accelerated collapse regime, Population III GEODEs with high spin are plausible. In fact, Bardeen (1970, Equation (4)) computes that ${a}_{* }\to 1$ as quickly as

$$\begin{eqnarray}&&m\to {m}_{i}\left\{1+3\left[\sqrt{{\sin }^{-1}\left(\displaystyle \frac{2}{3}\right)}-{\sin }^{-1}\left(\displaystyle \frac{1}{3}\right)\right]\right\},\end{eqnarray} \tag{ 68 }$$

where m_i is the Kerr mass at a_* = 0. In ideal settings, a 4 M_⊙ remnant with a_* = 0 will achieve a_* = 1 at $\sim 10{M}_{\odot }$, if reached via accretion.

It must be emphasized that the Kerr solution is asymptotically flat and has a mass parameter that does not blueshift. Following the discussion in Croker & Weiner (2019, Section 3.1), Equation (68) can only be usefully applied in scenarios where accretion occurs much more rapidly than cosmological blueshift. The above factor in mass gain occurs from blueshift alone at

$$\begin{eqnarray}&&{z}_{f}\simeq {\left(\displaystyle \frac{2}{3}\right)}^{1/3}{z}_{G}.\end{eqnarray} \tag{ 69 }$$

For z_G = 15, blueshift has contributed the same mass required from accretion by z_f = 13, or 60 Myr later. This is equivalent to a constant accretion rate of $\sim 2.5\times {10}^{-7}{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Accretion must be significantly larger than this, e.g., $\gtrsim \,{10}^{-5}{M}_{\odot }\,{\mathrm{yr}}^{-1}$, to apply Equation (68). This would imply GEODE γ ≫ 8 in as little as ∼1 Myr. In this regime, the burst formation approximation we have studied is justified.

Under slower accretion rates, the appropriate spin-up time is unclear. This justifies studying the dynamics of a z_G = 8 burst scenario, with the implicit understanding that the baryon budget is adjusted with respect to an earlier formation z_G ≳ 15. The detailed study of the Population III GEODE scenario with time-dependent $\beta (\eta )$ is the subject of a future paper.

The eigenframe fields are definitive: they are physical, and so do not depend on arbitrary choices of "the user," such as the metric ansatz. A key virtue of eigenframe fields is that they are orthonormal. In eigenframe fields, the unapproximated metric has the canonical Minkowski representation

$$\begin{eqnarray}&&{g}_{\mu \nu }=\mathrm{diag}(-1,1,1,1).\end{eqnarray} \tag{ B1 }$$

Pick any point P. All subsequent discussion will apply to frames at this point, so we will stop explicitly saying "at P." In an eigenframe, the metric representation is canonical Minkowski. If we perform a Lorentz transformation, the metric representation does not change. In general, we are no longer in an eigenframe, because our frame is now spacetime rotated relative to the eigenframe. We are, however, still in an orthonormal frame (as evidenced by the canonical metric representation). In general, every orthonormal frame is related to every other orthonormal frame by some Lorentz transformation. We will leverage this fact to connect the eigenframe to the RW frame.

Our first goal is to construct an orthonormal frame from any RW frame. We will call this frame the "associated vierbein frame." The procedure is as follows:

• 1.  
  Pick any RW frame for which the representation is defined (e.g., like Equation (5)).
• 2.  
  The lowered metric representation gives the RW frame basis overlaps:

  $$\begin{eqnarray}{g}_{\mu \nu }=\left[\begin{array}{ccc}{\partial }_{0}\cdot {\partial }_{0} & \cdots & {\partial }_{0}\cdot {\partial }_{3}\\ \vdots & \ddots & \vdots \\ {\partial }_{3}\cdot {\partial }_{0} & \cdots & {\partial }_{3}\cdot {\partial }_{3}\end{array}\right].\end{eqnarray} \tag{ B2 }$$

  These inner products can be used in the Gram-Schmidt procedure to determine an orthonormal basis {v_μ}, in terms of the RW frame coordinate basis {∂_μ}. This can always be done because the metric is semi-Riemannian and assumed to be well defined everywhere on ${ \mathcal M }$.
• 3.  
  Arrange the orthonormalized basis vectors (expressed with components from the RW frame representation) as rows in a table of numbers ${V}_{\sigma }^{\mu }$.

The change of basis matrix ${V}_{\sigma }^{\mu }$ is sometimes called the vierbein. In our setting, ${V}_{\sigma }^{\mu }$ transforms from the RW frame $\{{\partial }_{\mu }\}$ into the associated vierbein frame {v_μ}. The qualifier "associated" arises because the construction procedure takes as input some RW frame. Different RW frames will have different associated vierbein frames, in general.

Now that we have an associated vierbein frame, we may complete the transformation from the eigenframe to the RW frame. Because the associated vierbein frame is orthonormal, it is always related to the eigenframe by a Lorentz transformation. Let Λ be a boost representing the four-velocity from the eigenframe to the associated vierbein frame. Let V represent the associated vierbein and A be its inverse. Then,

$$\begin{eqnarray}&&{T}_{\mu \nu }={A}_{\mu }^{\lambda }{{\rm{\Lambda }}}_{\lambda }^{\sigma }\overline{{T}_{\rho \sigma }}{{\rm{\Lambda }}}_{\tau }^{\rho }{A}_{\nu }^{\tau }.\end{eqnarray} \tag{ B3 }$$

It can be verified that this procedure regenerates Bardeen (1980, Equations (2.17)–(2.20)) under three assumptions: the boost velocity is defined to be order , the source is cleaved into a background and an order perturbation, and one computes with the vierbein frame associated with Equation (5).

Source term representations in the RW frame fields can now be determined without explicit prescription of coordinate bases for the strong sources. No other knowledge is required except that the coordinates of the RW frame field remain well defined on open sets containing the source. In this sense, our procedure is a relativistic generalization of the procedure outlined by Bardeen (1980, Section IIA).

We now have sufficient machinery to determine the transformation required in Section 4.4. To enter the associated vierbein frame from the eigenframe, two separate boosts must be performed. Because boosts are closed under composition, Equation (B3) continues to apply. The first boost introduces the tangential velocity due to spin at P. The second boost introduces the peculiar center -of-mass-motion.

There is never any Hubble flow contribution, because all RW frame fields are built from comoving coordinates. The second boost is required because the "flow gauge" is designed to remove energy flux ${T}_{k}^{0}$ in RW frame fields that originally perceived it. Two separate boosts are required because vector four-velocity addition is incorrect in special relativity. The two successive boosts introduce a Thomas rotation and, in general, do not commute. This does not affect the isotropic pressure or energy density.

Finally, the inverse vierbein is applied to enter the flow frame. Clearly, this can only introduce order- adjustments. This procedure has been symbolically verified, in general, to produce ${T}_{\,k}^{0}$ with only spin contributions, as expected. The diagonal components combine as indicated in Equations (37) and (41).

The GEODE density contrast is defined as

$$\begin{eqnarray}&&{\delta }_{G}:= \displaystyle \frac{\delta {\rho }_{G}}{{\rho }_{G}},\end{eqnarray} \tag{ C1 }$$

where $\delta {\rho }_{G}$ is the GEODE density perturbation. This quantity is −1 in true vacuum and 0 in regions where the local density is equal to the background. The continuity equation at first order for a decoupled component is given by Hu (2004, Section 3.3, Equation (33)):

$$\begin{eqnarray}&&\left(\displaystyle \frac{d}{d\eta }+3\displaystyle \frac{\dot{a}}{a}\right)\delta {\rho }_{G}+3\displaystyle \frac{\dot{a}}{a}\delta {p}_{G}=-({\rho }_{G}+{p}_{G})({{kv}}_{G}+3\dot{{\rm{\Phi }}}).\end{eqnarray} \tag{ C2 }$$

Substituting the GEODE background evolution from Equation (17), the anisotropic stress from Equation (49), and the definition of the density contrast gives

$$\begin{eqnarray}&&\displaystyle \frac{d{\delta }_{G}}{d\eta }+3(1-\chi )\displaystyle \frac{\dot{a}}{a}{\delta }_{G}-2\displaystyle \frac{\dot{a}}{a}{{\rm{\Pi }}}_{G}=-\chi ({{kv}}_{G}+3\dot{{\rm{\Phi }}})\end{eqnarray} \tag{ C3 }$$

after evaluation of the derivative. Switching to scale factor and scaling gives the desired dynamical equation,

$$\begin{eqnarray}&&\displaystyle \frac{d{\delta }_{G}}{{da}}=3(\chi -1)\displaystyle \frac{{\delta }_{G}}{a}+\displaystyle \frac{2{{\rm{\Pi }}}_{G}}{a}-\chi \left(\displaystyle \frac{{{kk}}_{\mathrm{fac}}}{{{Ha}}^{2}}{v}_{G}+3\displaystyle \frac{d{\rm{\Phi }}}{{da}}\right).\end{eqnarray} \tag{ C4 }$$

We adopt the standard CDM evolution equations and regard baryons at z < 45 as essentially dust. This means that baryons and CDM will have the same evolution equations; their dynamics will only differ by their initial conditions set at z_G. We will now use the subscript m for matter.

The velocity evolution equation for matter is equal to Equation (53), with the subscript G replaced with m, χ → 1, and ${{\rm{\Pi }}}_{m}:= 0$:

$$\begin{eqnarray}&&\displaystyle \frac{{{dv}}_{m}}{{da}}=-\displaystyle \frac{{v}_{m}}{a}+\displaystyle \frac{{{kk}}_{\mathrm{fac}}}{{{Ha}}^{2}}{\rm{\Psi }}.\end{eqnarray} \tag{ C5 }$$

Again, we will compute with the density contrast

$$\begin{eqnarray}&&{\delta }_{m}:= \displaystyle \frac{\delta {\rho }_{m}}{{\rho }_{m}}.\end{eqnarray} \tag{ C6 }$$

Substitution of this definition into the continuity Equation (C2) gives

$$\begin{eqnarray}&&{\dot{\delta }}_{m}=-{{kv}}_{m}-3\dot{{\rm{\Phi }}}.\end{eqnarray} \tag{ C7 }$$

Switching to scale factor and rescaling units give the desired dynamical equation:

$$\begin{eqnarray}&&\displaystyle \frac{d{\delta }_{m}}{{da}}=-\displaystyle \frac{{{kk}}_{\mathrm{fac}}{v}_{m}}{{{Ha}}^{2}}-3\displaystyle \frac{d{\rm{\Phi }}}{{da}}.\end{eqnarray} \tag{ C8 }$$

Note that this equation agrees with Equation (C4) when $\chi := 1$ and ${{\rm{\Pi }}}_{G}:= 0$. This completes the specification of the material bulk degrees of freedom describing the CDM and baryon-tracking fluids at late times.

The Einstein equations are entirely determined via the action principle, the solution ansatz, and the choice of gauge. For ease of comparison with existing literature, we will express the dynamics in terms of the Newtonian scalar curvature Φ. Combine Hu (2004, Equations 17(a), (c)) and fix the gauge to Newtonian to find

$$\begin{eqnarray}&&{k}^{2}{\rm{\Phi }}=4\pi {{Ga}}^{2}{\rho }_{\mathrm{cr}}\displaystyle \sum _{J}\delta {\rho }_{J}+3\displaystyle \frac{\dot{a}}{a}\left(\displaystyle \frac{\dot{a}}{a}{\rm{\Psi }}-\dot{{\rm{\Phi }}}\right),\end{eqnarray} \tag{ C9 }$$

where we define the dimensionless sum of density perturbations as

$$\begin{eqnarray}&&\displaystyle \sum _{J}\delta {\rho }_{J}:= {{\rm{\Omega }}}_{m}{a}^{-3}{\delta }_{m}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}{a}^{-3\chi }{\delta }_{G}.\end{eqnarray} \tag{ C10 }$$

Hu (2004, Equation (17)) gives the algebraic relation between the gravitational potentials and the anisotropic stress perturbation:

$$\begin{eqnarray}&&{k}^{2}({\rm{\Psi }}+{\rm{\Phi }})=-8\pi {{Ga}}^{2}{p}_{G}{{\rm{\Pi }}}_{G},\end{eqnarray} \tag{ C11 }$$

where we have fixed to Newtonian gauge. Note that GEODEs are the only relevant source of anisotropic stress. Rearranging, substituting the background evolution, and scaling give

$$\begin{eqnarray}&&{\rm{\Psi }}=\displaystyle \frac{3(1-\chi ){a}^{2-3\chi }{{\rm{\Omega }}}_{{\rm{\Lambda }}}{{\rm{\Pi }}}_{G}}{{k}^{2}{k}_{\mathrm{fac}}^{2}}-{\rm{\Phi }}.\end{eqnarray} \tag{ C12 }$$

Rearranging, switching to scale factor, and scaling the remaining terms give the desired dynamical equation for the Newtonian scalar curvature:

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Phi }}}{{da}}=-\displaystyle \frac{{{ak}}^{2}{k}_{\mathrm{fac}}^{2}{\rm{\Phi }}}{3{H}^{2}{a}^{4}}+\displaystyle \frac{{a}^{3}{\displaystyle \sum }_{J}\delta {\rho }_{J}}{2{H}^{2}{a}^{4}}+\displaystyle \frac{{\rm{\Psi }}}{a}.\end{eqnarray} \tag{ C13 }$$

This completes the specification of the scalar gravitational degrees of freedom.