# Imprints of decaying dark matter on cosmic voids

Earl Lester<EMAIL_ADDRESS>Krzysztof Bolejko
<EMAIL_ADDRESS>School of Natural Sciences, College of Sciences
and Engineering, University of Tasmania, Private Bag 37, Hobart TAS 7001,
Australia

###### Abstract

The Standard Cosmological Model assumes that more than 85% of matter is in the
form of collisionless and pressureless dark matter. Unstable decaying dark
matter has been proposed in the literature as an extension to the standard
cold dark matter model. In this paper we investigate a scenario when dark
matter decays and the resultant particle moves with respect to the dark
matter. A covariant hydrodynamical model is developed in which the decay is
modeled by the transfer of energy-momentum between two dark dust fluid
components. We parameterise the model in terms of the decay rate $\Gamma$ and
injection velocity $v_{i}$ of the resultant dark matter particles. We apply
the framework to study the evolution of cosmic voids which are environments
with low content of baryonic matter. Thus, unlike baryon-rich environments,
voids provide an opportunity to measure dark matter signals that are less
contaminated by complex baryonic processes. We find that the growth of S-type
voids is modified by the dark matter decay, leading to imprints at the present
day. This paper serves as a proof-of-concept that cosmic voids can be used to
study dark mater physics. We argue that future cosmological observations of
voids should focus on signs of reported features to either confirm or rule out
the decaying dark matter scenario. Lack of presence of reported features could
put constraints of the decay of dark matter in terms of $\Gamma>H_{0}^{-1}$
and $v_{i}<10$ km/s.

###### pacs:

98.80.-k, 98.80.Es, 98.80.Jk

##  I Introduction

The status of Dark Matter (DM) is as an open problem in modern physics. As DM
is known to form the majority of the universal matter content, its presence
plays a important part in cosmological modelling and observations. Yet,
constraints on the possible particle properties of DM via direct or indirect
detection experiments have been inconclusive Mack (2014), and DM particles
have never been directly detected in an experiment on Earth Craig (2013);
Akerib et al. (2014). Cosmology and astrophysics offer a laboratory of sorts
to explore the effects of particle physics, Dolgov and Zeldovich (1981).

In the standard cosmological model the source of the gravitational field is a
perfect fluid and dark matter is often modelled as a pressureless and
collisionless weakly interacting massive particle (WIMP) which in the light of
experimental data (or rather lack of thereof) will need to be replaced with
some other candidate Fox et al. (2014). However, it is well known that the
standard Cold Dark Matter (CDM) model seems to be inadequate on sub-galactic
scales ($<1$ Mpc) Boehm et al. (2002). This has motivated investigations into
modifications to the standard model such as Warm Dark Matter (WDM) which damps
the production of low-mass objects, Schaeffer and Silk (1988). Alternatively,
interactions in the dark sector have been investigated in recent years. Models
of DM interactions include annihilation, scattering or decay (DDM) and offer
simple extensions of the standard WIMP models which retain the desirable
properties of CDM while offering possible resolutions to known tensions and
discrepancies between theory and observation, Mack (2014); Boehm et al.
(2002); Wang et al. (2014). Strongly self-interacting modifications to CDM
include collisional DM suggested by Spergel and Steinhardt (2000) to account
for flat density profiles of galaxy cores but ruled out by Arabadjis et al.
(2002); Boehm and Schaeffer (2005); and self-annihilation with a quadratic
dependence on density is prevalent in galaxy cores, Mack (2014), but must be
carefully contrived to avoid complete self-annihilation in the early universe,
Craig and Davis (2001). Scattering processes such as DM-photon or DM-neutrino
interactions have been shown to have similar effects to warm DM models Boehm
et al. (2002) while also producing horizon-scale weak damping effects.

The effect of DM physics on cosmological structures offers a way of using
cosmological observations to gain insight into possible interactions in the
dark sector. Such observables include high redshift data such as cosmic dark
ages and cosmic microwave background (CMB), Slatyer and Wu (2017), or low-
redshift data such as the effects of dark matter on galactic halos Mack
(2014). However astrophysical signals such as, for example, DM annihilation
Mack (2014); Lacroix et al. (2016a) have been shown to be difficult to
distinguish from signals produced by baryonic processes Lacroix et al.
(2016b); Clark et al. (2018) making direct detection of the products of
annihilation at best inconclusive. Interactions could introduce additional
effects that might influence the overall evolution of the universe. These
effects, such as decay and scattering, could lead to energy or momentum flux,
or viscosity making the fluid description of dark matter imperfect Ellis
(1973); Ellis and Wainwright (1997); Ellis et al. (2012).

The focus of the present paper is decaying dark matter (DDM) which has a
linear density dependence and may decay into two- or many- bodies with
massless decay by-products, and may feature kick or injection velocities of
relativistic order, Blackadder and Koushiappas (2016). Hence DDM may produce
standard model particles, photons or exotic smaller mass dark particles or
dark radiation, Pandey et al. (2020); Ackerman et al. (2009), and energy-
injection by cascade effects, Slatyer and Wu (2017), offering possible means
of indirect detection. Unstable DDM is a natural development beyond the
canonical stable WIMP model.

One of the primary benefits of DDM is the early-time correspondence of the
model to CDM and consistency with the CMB anisotropy spectrum while late-time
tensions within the CDM framework such as the missing satellite problem are
offered potential resolutions by the long life-time instability of the DM
particle, Wang et al. (2014), – assuming a particle lifetime of order the
current age of the universe or longer. In this scenario, the observed DM
subhaloe underabundancy may be remedied by both decay into relativistic, Cen
(2001), or non-relativistic daughter particles, Wang et al. (2014). Similarly,
the correspondence of DDM with CDM at high-redshift offers possible
alleviation to, if not resolution of, the observed tensions in $H_{0}$ as
inferred from CMB and Type-1 Supernovae data, Vattis et al. (2019), as well as
the $\sigma_{8}$ matter fluctuations, Pandey et al. (2020). However, recent
work on relativistic single-species decay products have investigated the
similarity of the induced bulk viscosity and modified cosmic acceleration with
the negative pressure required of dark energy, Wilson et al. (2007), and the
signals constraining DDM physics may be degenerate to other model
modifications such a modified Newtonian dynamics, Peter et al. (2010a).

The present work proposes to investigate the possible decay of DM into other
dark particles via the effect this may have upon the growth of large-scale
structure, specifically cosmological voids. Voids are chosen as these offer an
unique environment with several efficacious features. The under-abundance of
galaxies in these avoided regions suggests an inhomogeneous DM under-density,
Peebles (2001), assuming a strong Light-to-Mass ratio, and therefore cosmic
voids will have less significant annihilation or self-interaction effects due
to the quadratic mass dependency of these. Furthermore, these regions have
minimal baryonic contamination so that scattering and baryonic feedback
processes may be negligible. Thus, in cosmic void regions the only remaining
particle interaction which may be introduced to modify DM and produce
significant signals, is decay.

As voids are the complementary structure of the filaments and sheets of the
cosmic web, realistic analysis cannot consider voids as isolated objects, van
de Weygaert (2014). The environmental dependence and hierarchical structure
formation of voids creates a complex morphology which may conveniently be
classified by the dominant evolutionary process determined by the surrounds:
either expansion or collapse, Sheth and Van De Weygaert (2004). Larger voids
are observed to generally fit the former scenario corresponding to a void-in-
void environmental configuration. The density profile of these voids rises
smoothly towards the mean background density and are classified as R-type in
recent literature. Alternatively, the density profiles of smaller, S-type
voids exhibit a compensating shell suggesting a void-in-cloud type
configuration, Ceccarelli et al. (2013). In both the uncompensated and
compensated cases, there is a general trend of spherical voids evolving
towards a bucket-shape density profile, van de Weygaert (2014). Furthermore,
the smoothing of asphericities as under-densities grow, Icke (1984), validates
the approximation of spherically symmetric voids. Though not explored here,
the ellipticity of voids is an important measure of tidal effects induced by
the surrounding environment and offers an significant probe into cosmology and
dark sector physics, Rezaei (2019, 2020). Despite the intimate coupling of
voids with one another and the overall environment, the standard first-
approximation of an isolated, spherically symmetric under-density is often
adopted; we also do this here since it offers an immediate quantitative and
qualitative appreciation of the effect of proposed dark matter decay upon void
evolution.

As these cosmological voids are inhomogeneous and evolve non-linearly, Bolejko
(2006), the cosmology must move beyond the standard spatially homogeneous and
isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) model. Perturbative
approaches may be employed but require a stable background so cosmological
back-reaction may affect the accuracy of these models, and furthermore, gauge
issues may arise Ellis and Bruni (1989). The Lemaître-Tolman-Bondi (LTB) class
of exact spherically symmetric dust solutions to the Einstein equations are
inhomogeneous but isotropic and are appropriate for simplified void models as
pressure is generally negligible until shell-crossing occurs, Lemaître (1933,
1997); Tolman (1934, 1997); Bondi (1947). The Szekeres models are both
inhomogeneous and anisotropic, Krasiński (1997) and have been employed in the
modelling of cosmological voids. These methods are however constrained by the
requirement of being metric based. Recently, tetra-based methods have been
developed as an alternative to LTB models, Kim et al. (2019). The well-
established semi-tetrad approach of the hydrodynamical 1 + 3 covariant
methodology has an intuitive fluid flow interpretation, Ellis (1973), and when
taken to have local rotational symmetry (LRS) or silence (vanishing magnetic
Weyl and pressure gradients) offers a generalisation to th LTB models in
certain limits van Elst and Ellis (1996); Van Elst et al. (1997). In this
paper we adopt a fully general relativistic hydrodynamical 1 + 3 formulation
with LRS, to which further simplifications are imposed, see Section II.

The decay effect is proposed to be modeled by the one-way transfer of energy-
density via long-lifetime nuclear-like decay from a dust fluid comoving with
the fundamental observer (representative of the CDM, with baryonic matter
deemed insignificant and disregarded for the purposes of this model) into a
secondary non-comoving fluid. This is intended to represent a simple
exponential decay of the massive parent particle into a massive daughter
particle.

The average relative spatial velocity of the second non-comoving dust
generated by the decay induces an effective imperfect fluid EM tensor in the
frame of the fundamental observer. Tilted models wherein imperfect fluid
effects arise have been studied in the literature, see e.g. Coley (1990). The
presence of these imperfect terms corresponding to heat flux, pressure and
viscosity induces non-FLRW effects such as shear or Weyl curvature (which are
identically zero for the FLRW models), which in turn modifies the overall
evolution of the cosmology and structure therein. However, for small peculiar
velocities the EM tensor approaches that of a dust with heat conduction,
investigated earlier in Goode (1986) where an explicit form for the space-time
metric was determined, and more recently in Sebastián and Sussman (2021),
wherein the momentum-transfer was interpreted as the effect of a non-comoving
fluid with non-relativistic peculiar velocity, in accordance with the present
work.

The structure of this paper is as follows. In Section II, the tilted two-
fluid, 1 + 3 covariant model is developed which leads to a governing system of
PDEs for a set of scalars once LRS is imposed. This system of equations is
then employed in Section III to model spherically symmetric, isolated voids of
the two types discussed above. The results of numerical analysis of this
system are presented and discussed, with an emphasis on the possibility of
using these results as an indirect means of constraining DDM. Limitations of
our model are also discussed, and the article then concludes in Section IV.

##  II Methods

###  II.1 The 1+3 covariant framework

The tilted 1 + 3 covariant Lagrangian formalism is employed. The model is
chosen primarily owing to its direct physical interpretation of the kinematic
variables adopted, and, secondly, for its proven capability of addressing non-
linear evolution of which voids are a weak example, Maartens et al. (1999). We
also remark on its prevalence in the literature of relativistic cosmology.

Following e.g. Ellis (1973); Ellis and Van Elst (1999); Tsagas et al. (2008),
flow-orthogonal projections are made to the spacetime by projection into the
instantaneous rest space of an observer using

$h^{ab}=g^{ab}+u^{a}u^{b}$

where $g^{ab}$ is the metric tensor and $u^{a}$ the comoving 4-velocity. The
gradient of the 4-velocity field of the comoving observer is decomposed into
vorticity, $\omega_{ab}$, shear $\sigma_{ab}$, expansion $\Theta$ and
acceleration $A_{a}$ such that

$u_{a;b}=\omega_{ab}+\sigma_{ab}+\dfrac{1}{3}\Theta
h_{ab}-A_{a}n_{b}\,\text{.}$ (1)

As mentioned in the introduction, some DDM models postulate the production of
standard matter particles which may be directly detected. However, in the
present investigation, we wish to determine possible means of indirect
detection of DDM via the imprint on cosmological voids, if any. Thus, for the
purposes of this preliminary work, it is sufficient to assume that the
secondary fluid produced via the decay is another dark fluid.

This work proposes a model involving two dust components: the first is
comoving with the fundamental observer, has density $\rho$ and a normalised
comoving 4-velocity field $u^{a}=\delta^{a}_{t}$, and decays with a long-life
time into the second non-comoving resultant fluid with density $\epsilon$ and
4-velocity $v^{a}$. The energy-momentum tensor for the total fluid is taken to
be the sum of these dust components such that
$T^{ab}=\Sigma_{(i)}T^{ab}_{(i)}$ and hence,

$T^{ab}=\rho u^{a}u^{b}+\epsilon v^{a}v^{b}\,\text{.}$ (2)

The non-comoving velocity is tilted with respect to the first by the Boost

$v^{a}=\gamma(u^{a}+V^{a})$ (3)

with Lorentz factor $\gamma=(1-V_{a}V^{a})^{-1/2}$ and spacelike peculiar
velocity $V^{a}$, King and Ellis (1973).

Hence, an observer moving with $u^{a}$, i.e. comoving with the first fluid
with density $\rho$, measures the total density to be

$\mu=T^{ab}u_{a}u_{b}=\rho+\eta,$

where $\eta=\gamma^{2}\epsilon$, the total pressure to be

$p=\frac{1}{3}h_{ab}T^{ab}=\frac{1}{3}\eta h_{ab}V^{a}V^{b},$

the heat flux to be

$q^{a}=-(\mu u^{a}+u_{a}T^{ab})=\eta V^{a},$

and the anisotropic stress as

$\pi_{ab}=T_{cd}h^{c}_{\langle a}h^{d}_{b\rangle}.$

Thus, a comoving observer does not discern two dust components, but rather the
source of the gravitational field is apparently a single imperfect fluid,
which may be described by

$T^{ab}=\mu u^{a}u^{b}+ph^{ab}+2q^{(a}u^{b)}+\pi^{ab}.$ (4)

The assumption of LRS takes there to exist locally an axis of symmetry or
preferred spatial direction induced by some properties of the flow (e.g. shear
or vorticity), such that at each event the spacetime exhibits partial
isotropy, Goode and Wainwright (1986). Observations become independent of
spatial dimensions orthogonal to the preferred spatial direction, Singh et al.
(2016, 2018). Hence, all spacelike covariant vector fields in the LRS
spacetime are proportional to the preferred spacelike unit direction vector.
Let this preferred spacelike unit vector be denoted by $z^{a}$ and note that

$z^{a}z_{a}=1,\quad u^{a}z_{a}=0.$

Two new variables may be introduced to describe the evolution of $z^{a}$. The
magnitude of the divergence $\alpha$ is taken to be

$\alpha=D_{a}z^{a},$

where $D^{a}X=h^{ab}X_{;b}$ is the projected covariant derivative. The
magnitude of the spatial rotation, or twist, $\tau$, is determined by

$\tau=-z_{a}\epsilon^{abcd}u_{d}z_{c;b},$

where $\epsilon^{abcd}$ is the Levi-Cita symbol. Note that it follows that,
van Elst and Ellis (1996),

$\dot{z}^{a}=0,\quad z^{a}_{;a}=\alpha+A^{a}z_{a},$

where $\dot{X}=u^{a}X_{;a}$.

From the preferred spacelike unit vector $z^{a}$ is defined the symmetric
tensor

$e_{ab}=\frac{1}{2}(3z_{a}z_{b}-h_{ab}),$

which has the properties, van Elst and Ellis (1996), that

$\dot{e}_{ab}=0,\quad D_{b}e^{ab}=\frac{3}{2}\alpha z^{a},\quad
e^{a}_{b}e^{b}_{a}=\frac{3}{2}.$

By LRS, in order to reduce the above equations to a scalar system, following
van Elst and Ellis (1996) and Bolejko (2017), all spacelike covariant vector
fields may be represented as

$A^{a}=Az^{a},\quad V^{a}=Vz^{a},\quad\omega^{a}=\omega z^{a}$

from which it follows that

$q^{a}=Qz^{a},\quad\pi^{ab}=\Pi e^{ab}$

and

$p=\frac{1}{3}\eta V^{2},\quad Q=\eta V,\quad\Pi=\frac{2}{3}\eta V^{2}.$

In a similar fashion as above, LRS allows all the irreducible kinematic
components of the shear and Electric and Magnetic parts of the Weyl tensors to
be written in terms of scalars such as

$\sigma^{ab}=\Sigma e^{ab},\quad E^{ab}=We^{ab},\quad H^{ab}=Be^{ab}.$

The substitution of the above into the 1 + 3 covariant governing equations
(see e.g. Tsagas et al. (2008), and, adopting, for this section, units such
that $c=8\pi G$), produces the following system of propagation equations:

$\displaystyle\dot{\Theta}=-\dfrac{1}{3}\Theta^{2}-\dfrac{1}{2}(\mu+3p)-\frac{3}{2}\Sigma^{2}$
$\displaystyle\quad\quad+3\omega^{2}+A^{{}^{\prime}}+\alpha
A+A^{2}+\Lambda\,\text{,}$ (5)
$\displaystyle\dot{\Sigma}=-\dfrac{2}{3}\Theta\Sigma-\frac{1}{2}\Sigma^{2}-\frac{2}{3}\omega^{2}-W+\frac{1}{2}\Pi+\frac{2}{3}A^{{}^{\prime}}-\frac{1}{3}A\alpha+A^{2}\,\text{,}$
(6) $\displaystyle\dot{\omega}=-\frac{2}{3}\Theta\omega+\Sigma\omega+\tau
A\,\text{,}$ (7) $\displaystyle\dot{W}+\frac{1}{2}\dot{\Pi}=-\Theta
W-\dfrac{1}{2}(\mu+p)\Sigma-\frac{1}{6}\Theta\Pi+\frac{3}{2}\Sigma W$
$\displaystyle\quad\quad-\frac{1}{4}\Sigma\Pi+\frac{1}{6}Q\alpha-\frac{1}{3}Q^{{}^{\prime}}\,\text{,}$
(8) $\displaystyle\dot{B}=-\Theta B+\frac{3}{2}\tau
W-\frac{3}{4}\tau\Pi+\frac{3}{2}B\Sigma-\omega Q\,\text{,}$ (9)
$\displaystyle\dot{\mu}=-\Theta(\mu+p)-\frac{3}{2}\Pi\Sigma-Q^{{}^{\prime}}-Q\alpha-2QA\,\text{,}$
(10)
$\displaystyle\dot{Q}=-(\Sigma+\frac{4}{3}\Theta)Q-\frac{3}{2}\Pi\alpha-\Pi^{{}^{\prime}}-p^{{}^{\prime}}-\Pi
A-(\mu+p)A\,\text{,}$ (11)

where $X^{\prime}=z^{a}X_{;a}$. Moreover, we acquire the constraints:

$\displaystyle
Q=\dfrac{2}{3}\Theta^{{}^{\prime}}-\Sigma^{{}^{\prime}}-\frac{3}{2}\alpha\Sigma-\tau\omega\,\text{,}$
(12) $\displaystyle B=-\frac{3}{2}\tau\Sigma+\omega(2A-\alpha)\,\text{,}$ (13)
$\displaystyle
W^{{}^{\prime}}=-\frac{3}{2}W\alpha+\frac{1}{3}\mu^{{}^{\prime}}-\frac{1}{2}\Pi^{{}^{\prime}}-\frac{3}{4}\Pi\alpha-\frac{1}{3}\Theta
Q$ $\displaystyle\quad\quad+\frac{1}{2}\Sigma Q-3B\omega\,\text{,}$ (14)
$\displaystyle B^{{}^{\prime}}=-\frac{3}{2}\alpha
B+(\mu+p)\omega+\frac{1}{2}\tau Q+3\omega W-\frac{1}{2}\Pi\omega\,\text{,}$
(15) $\displaystyle\omega^{{}^{\prime}}=(A-\alpha)\omega\,\text{.}$ (16)

###  II.2 Pressureless, irrotational and dissipative fluids

In order to further simplify the system the stringent conditions of
irrotationality i.e. vanishing vorticity, $\omega=0$, and twist, $\tau=0$, are
assumed. The former ensures the existence of a well-defined cosmological time.
It follows via the constraint equation (13) that a twist-free, irrotational
flow is necessarily purely gravito-electric such that $B=0$.

The vanishing of the effective pressure scalars $p$ and $\Pi$ and heat scalar
$Q$ immediately follows in the limit of $V\to 0$. In this situation, the model
reduces to the well-know spherically symmetric LTB cosmology. As mentioned in
the introduction, exact homogeneous solutions featuring dissipative dust were
found and analysed in Goode (1986). Early research explored inhomogeneous dust
models with heat flow, generalising the subclasses of inhomogeneous and
anisotropic Szekeres solutions de Lima and de Garcia Maia (1985). However, as
recently pointed out by Nájera and Sussman (2020), the heat flow of
dissipative dust models is inappropriate for standard CDM and late-time
universe models, and is analogous to the energy flux induced in the observer
frame by a non-comoving, non-relativistic peculiar velocity. The two-fluid
tilted model developed here has produced an effectively imperfect total EM
tensor.

Observational analyses of large-scale peculiar velocities suggest non-
relativistic orders. It therefore seems valid to assume a non-relativistic,
non-comoving peculiar velocity of the second fluid such that $\gamma\approx
1$, and the effective pressure scalars $p$ and $\Pi$ become second-order and
negligible. The effective heat scalar $Q$, however, is not. The tractability
of the mathematical model is further increased by assuming that the
covariantly projected dot and prime derivatives of the effective pressures are
also first order and negligible. This is justified as the initial velocity
distribution is taken to be homogeneous, and hence, all terms featuring $p$ or
$\Pi$ vanish.

The above assumptions lead to the set of propagation equations:

$\displaystyle\dot{\Theta}=-\dfrac{1}{3}\Theta^{2}-\dfrac{1}{2}\mu-\frac{3}{2}\Sigma^{2}+A\alpha+A^{2}+\Lambda\,\text{,}$
(17)
$\displaystyle\dot{\Sigma}=-\dfrac{2}{3}\Theta\Sigma-\frac{1}{2}\Sigma^{2}-W-\frac{1}{3}A\alpha+\frac{2}{3}A^{2}\,\text{,}$
(18) $\displaystyle\dot{W}=-\Theta W-\dfrac{1}{2}\mu\Sigma+\frac{3}{2}\Sigma
W+\frac{1}{6}Q\alpha-\frac{1}{3}Q^{\prime}-\frac{2}{3}AQ\,\text{,}$ (19)
$\displaystyle\dot{\mu}=-\Theta\mu-Q^{\prime}-Q\alpha-2AQ\,\text{,}$ (20)
$\displaystyle\dot{Q}=-(\Sigma+\frac{4}{3}\Theta)Q-\mu A\,\text{,}$ (21)

and the constraints:

$\displaystyle
Q=\dfrac{2}{3}\Theta^{\prime}-\Sigma^{\prime}-\frac{3}{2}\alpha\Sigma\,\text{,}$
(22) $\displaystyle
W^{{}^{\prime}}=-\frac{3}{2}W\alpha+\frac{1}{3}\mu^{\prime}-\frac{1}{3}\Theta
Q+\frac{1}{2}\Sigma Q\,\text{.}$ (23)

Therefore, to close this system an evolution equation for the unknown spatial
divergence $\alpha$ is required. This is determined by application of the
Voss-Weyl divergence formula to the comoving 4-velocity $u^{a}$, and to the
preferred spatial direction $z^{a}$, giving, respectively

$\Theta=\frac{\dot{\nu}}{\nu}\quad\text{
and}\quad\alpha=\frac{\nu^{\prime}}{\nu}\,\text{,}$ (24)

where $\nu=\sqrt{-g}$, with $g$ representing the metric determinant.

Hence, the system formed by Equations (17) - (24), serves as a closed set of
evolution equations which, given appropriate initial conditions, may be
numerically marched forward in time. The constraints may be helpful as an
alternate substantiation of the results of evolving the system.

In addition, one can define the scale factor $a$

$a=\left(\frac{\nu}{\nu_{i}}\right)^{1/3}.$ (25)

In general the scale factor defined in this way is a function of both time and
position but in the FLRW limit the scale factor becomes a function of time
only, i.e. $a=a(t)$ and $a^{\prime}=0$.

###  II.3 FLRW limit

In order to facilitate a comparison with the standard model and a discussion
of the modifications introduced, briefly consider necessary and sufficient
conditions for a spacetime to be FLRW Krasiński (1997):

  1. 1.

the metric obeys the Einstein equations with a perfect fluid source; and

  2. 2.

the velocity field of the perfect fluid source has zero rotation, shear and
acceleration.

The first condition requires

$Q=0,$

while the second condition gives

$\Sigma=0,{\rm~{}~{}~{}and~{}~{}~{}}A=0.$

Since these requirements must be fulfilled at all times, it follows from (18)
that

$W=0.$

This is expected as the FLRW models are conformally flat (i.e. the Weyl
curvature vanishes). These conditions, together with (14) require that
$\mu^{\prime}=0$, which in turn via (20) implies that $\Theta^{\prime}=0$,
meaning all gradients vanish, including $a^{\prime}=0$.

Using (24) and (25) it follows that the Hubble parameter

$\frac{\dot{a}}{a}\equiv H=3\Theta.$ (26)

This implies that the Raychaudhuri equation (17) reduces to the 2nd Friedmann
equation, which together with the continuity equation (20) yields the 1st
Friedmann equation

$\frac{\dot{a}^{2}}{a^{2}}=-\frac{Kc^{2}}{a^{2}}+\frac{8\pi
G}{3}\rho+\frac{1}{3}\Lambda c^{2},$ (27)

which is usually conveniently rewritten at the present instant ($t=t_{0}$) in
terms of the $\Omega-$parameters

$\Omega_{M}+\Omega_{\Lambda}+\Omega_{K}=1,$ (28)

where

$\displaystyle\Omega_{K}=-\frac{K}{H_{0}^{2}a_{0}^{2}},$ (29)
$\displaystyle\Omega_{M}=\frac{8\pi G}{3H_{0}^{2}}\,\rho_{0},$ (30)

and

$\displaystyle\Omega_{\Lambda}=\frac{1}{3H_{0}^{2}}\,\Lambda.$ (31)

Cosmological parameters determined by the 2015 Planck observations are adopted
so that $\Omega_{K}=0$, $\Omega_{M}=0.308$, $\Omega_{\Lambda}=0.692$ and
$H_{0}=67.81$ km s-1 Mpc-1, Planck Collaboration et al. (2016).

###  II.4 Initial conditions

The model begins from the last scattering instant (CMB), the redshift of which
is taken to be $z_{CMB}=1090$. At the initial instant, it is assumed there
exists only a single, comoving dust fluid. Then $V=0$ and the system is
approximated by an exact spherically symmetric solution of the Einstein
equations for the gravitational field sourced by dust only, i.e. the LTB
model. In this, the fluid quantities are

$\displaystyle\Theta=\frac{\dot{R}^{\prime}}{R^{\prime}}+2\frac{\dot{R}}{R},$
(32)
$\displaystyle\Sigma=-\frac{1}{3}\left(\frac{\dot{R}^{\prime}}{R^{\prime}}-\frac{\dot{R}}{R}\right),$
(33) $\displaystyle{\cal
W}=\frac{1}{3}\left(\frac{M}{R^{3}}-\frac{M^{\prime}}{R^{2}R^{\prime}}\right),$
(34) $\displaystyle\mu=\frac{2M^{\prime}}{R^{2}R^{\prime}},$ (35)
$\displaystyle Q=0.$ (36)

The evolution equations can be reduced to a single ODE Bolejko et al. (2009)

$\dot{R}^{2}=-K(r)+\frac{2M(r)}{R}+\frac{1}{3}\Lambda R^{2}.$ (37)

Thus, to specify the LTB model, and in our case the initial condition, one is
required to provide the two functions: $M(R)$ and $\dot{R}(R)$.

The two void configurations discussed in the introduction will be explored in
this work by specifying initial conditions for these types separately. In both
cases it is assumed that the radial coordinate coincides initially with the
function $R$, i.e.

$\displaystyle
r=\left.R\right|_{t_{i}}~{}~{}\Rightarrow~{}~{}\left.R^{\prime}\right|_{t_{i}}=1$
(38) $\displaystyle t_{i}=t_{CMB}.$ (39)

The above is equivalent to defining the scale factor to be $a_{i}=1$ when
$t=t_{i}$. The standard convention is to set the scale factor at present to
unity and put $a(z)=1/(1+z_{CMB})$. However, this scaling does not have any
physical meaning, so does not influence the evolution of the system, and it is
more convenient to set the scale factor to unity at the CMB (when the model is
almost homogeneous) than at the present time (when the system is far from
FLRW).

Due to the system being not FLRW the evolution of the radius $R$ does not
follow a simple FLRW scaling so that $R\neq ar=r(1+z)/(1+z_{CMB})$. However,
for the structures considered here (cosmic voids and present-day scales of
Mpc), when calculating the physical distance from the origin to a considered
point, this is approximated by

$R(r,t)\approx\int\limits_{0}^{r}\,{\rm dr}\,a(t,r),$ (40)

where the value of the scale factor $a(t,r)$ follows from (25). The above
approximation neglects the contribution from the spatial curvature $K$. While
the function $K$ is known at the initial instant, dissipate processes mean
that it is not a function of $r$ (as in LTB models) but evolves in time
Bolejko et al. (2011).

It is important to emphasise that the value of $R$ does not affect the
evolution of the system governed by (17) – (21). Rather, it is only used for
plotting purposes to represent graphically the physical quantities, i.e.
rather than plotting fluid variable (i.e. $\mu,\Theta,\Sigma,W$, or $Q$) as a
function of r (in Kpc), in Sec. III these quantities are plotted in terms of
$R$, as given by (40).

Figure 1: Density profile at $z=0$ of the evolved S-type void with (blue solid
lines) and without (red dashed lines) decay, for various combinations of the
parameters $\Gamma$ and $v_{i}$. This panes clearly illustrate the injection-
velocity-dependent growth of novel secondary structure at the edge of the
under-dense region.

####  II.4.1 S-type voids

The S-type voids are surrounded by a compensating shell, i.e. the mass of the
system outside the void is approximately the same as the mass of a region of a
similar size but otherwise homogeneous. For the S-type voids, the initial
profile for the function $M$ is taken to be

$M(r)=\frac{1}{6}\rho_{CMB}\left[1+\frac{1}{2}m_{0}\left(1-\tanh\frac{r-r_{0}}{2\Delta
r}\right)\right]r^{3},$ (41)

with $m_{0}=-0.006$, $\Delta r=1.3$ [Kpc], and $r_{0}=6.0$ [Kpc],
$\rho_{CMB}=(1+z_{CMB})^{3}\Omega_{M}3H_{0}^{2}/(8\pi G)$ and
$z_{CMB}=1090.0$. The function $K$ is fixed by the condition of a uniform age
of the universe $t_{B}=0$, so that it corresponds to existence of pure growing
modes Bolejko et al. (2009); Sussman (2013). This set of parameters leads to a
void that at its central part has a density contrast of $\delta_{0}=-0.9$ (cf.
Fig. 1).

####  II.4.2 R-type

A void of R-type has a continuously rising density profile, without a
compensating shell. Such a shape may be modelled with a pure Gaussian density
fluctuation, and then it follows from (35), that the function $M(r)$ when
$t=t_{i}$, is

$M(r)=\frac{1}{2}\rho_{CMB}\int\limits_{0}^{r}{\rm
d}\tilde{r}\,\tilde{r}^{2}\left[1+m_{0}\exp\left(-\frac{\tilde{r}^{2}}{\sigma^{2}}\right)\right],$
(42)

where $\rho_{CMB}$ is as above, and $m_{0}=-0.0033$, $\sigma=24$ [Kpc]. As
above, the function $K$ is fixed by the condition of pure growing modes
perturbations (Bolejko et al., 2009; Sussman, 2013). This set of parameters
leads to a void that at its central part has a density contrast of
$\delta_{0}=-0.8$ (cf. Fig. 2).

Figure 2: Density profile at $z=0$ of the evolved R-type void with (blue solid
lines) and without (red dashed lines) decay, for various combinations of the
parameters $\Gamma$ and $v_{i}$. Unlike in the S-type case, there is only
apparent a slight shallowing of the centre of the R-type under-density. Figure
3: Present day form of the non-dimensionalised $H$ (see eq. (26)), shear
scalar $\Sigma$, and electric Weyl scalar $W$ for the S-type void with
parameters $\Gamma=1/H_{0}$ and $v_{i}=7$ km/s. There is an insignificant
difference in the scaled Hubble parameter between the decaying and non-
decaying models, while significant changes in the latter two parameters.

###  II.5 Dark Matter physics

The gravitational field in the present model is sourced by DM only which is
assumed to consist of two fluids: the original comoving dark dust fluid that
decays into a secondary dark dust fluid which is non-comoving. In order to
account for the interaction, the conservation equations for each fluid must be
satisfied in the form

$T^{ab}_{(i)\,;b}=I^{a}_{(i)}$ (43)

where $I^{a}_{(i)}$ is the interaction term such that the total fluid is
conserved, $\Sigma_{(i)}I^{a}_{(i)}=0$, i.e. $I_{(2)}=-I_{(1)}$.

For this model of decay the interaction term is assumed to have the following
form

$I_{(1)}^{a}=-\Gamma\rho(u^{a}+w^{a})\,\text{.}$ (44)

where $w^{a}=v_{i}z^{a}$ is a vector that describes initial velocity of the
second fluid immediately after decay.

The temporal component of (43) is a typical decay law. Indeed, using
$\rho=N/V$ and $\dot{V}/V=\Theta$, this equation can be re-written as

$\dot{N}=-\Gamma N,$

which is a typical decay law of exponential (or nuclear) decay. The spatial
component of (43) allows us to evaluate the acceleration $A$ since

$\rho A-z_{a}I^{a}_{(1)}=0~{}~{}\Rightarrow~{}~{}A=\Gamma v_{i},$ (45)

where $v_{i}$ is the initial (or the injection or kick) velocity of the
particles of the second fluid with respect to the decaying fluid. In practice
$v_{i}$ is the average velocity of the ‘injected’ particles. Therefore, this
model contains two parameters: the decay rate $\Gamma$ and and injection
velocity $v_{i}$.

Figure 4: Present-day form of the non-dimensionalised $H$ (see eq. (26)),
shear scalar $\Sigma$, and electric Weyl scalar $W$ for the R-type void with
parameters $\Gamma=1/H_{0}$ and $v_{i}=7$ km/s. Figure 5: Evolution of density
profile for S-void at different redshift and for different values of the decay
rate $\Gamma$. Top panels show the void without the feature present, while
bottom panels show voids with the feature present, thus allowing to link the
lowest redshift where no feature is present with the properties of dark
matter. Due to the degeneracy between $\Gamma$ and $v_{i}$ in order to use
observations of voids at different redshifts to constrain the decay rate a
further information about $v_{i}$ (or equivalently particle mass ratio and
energy released at the decay) is required.

Investigations of self-gravitating halos of nearly-degenerate, non-
relativistic DDM have produced observational constraints that rule out the
combination of decay lifetimes $\Gamma<H_{0}^{-1}$ Gyr and injection
velocities of $v_{i}>20$ km/s Peter et al. (2010a). Analysis of the observed
mass-concentration and cluster mass function has constrained the DDM
parameter-space of decay lifetimes and kick velocities for numerical
simulations of DDM, suggesting an allowed region of lifetimes from a few times
longer than the age of the universe to less than $100H_{0}^{-1}$ and kick
velocities $1\,\text{km/s}<v_{k}<100\,\text{km/s}$, Peter et al. (2010b).
Other indirect constraints on the parameters of generic DDM models have been
developed via the luminosity distance of supernovae standard candles which has
a cosmological and energy-density dependence, Blackadder and Koushiappas
(2014). These results suggest that in the case that 90% of the mass is
transferred to the massive daughter particle, the decay lifetime has a lower
limit of $8.4$ Gyr with 95% confidence. Alternatively, when 15% is
transferred, the lower limit on the lifetime is 25 Gyr, Blackadder and
Koushiappas (2016). Other observational constraints on DDM models have been
studied in the literature by means such as cosmological weak lensing effects
Wang and Zentner (2012) or the influence on reionisation via photons emitted
during the decay process Mapelli et al. (2006). In view of this preceding
research, we propose to investigate a range of decay parameters within similar
regions. We investigate the region of parameter space relating to decay
lifetimes of between one and a few times the age of the universe, and small
non-relativistic injection velocities.

##  III Results

Starting from the initial conditions outlined above, the governing system of
propagation equations, (17) – (21), was time-marched forward to the present
day, $z=0$. The results of the evolution of the S-type void is presented in
Fig. 1 and the R-type void in Fig. 2. The density profiles shown in these
figures are normalised by the present day background density given by eq.
(30). Figure 1 illustrates the evolution of the S-type void over a range of
decay-rates $\Gamma$ and injection velocities $v_{i}$. Overlaid in each plot
as a red dashed line are the results of the evolution of the density profile
without decay, $\Gamma=0$, and thus no non-comoving second fluid nor effective
heat/momentum flow (i.e. such a model reduces to the LTB model).

The results summarised in Fig. 1 indicate that the heat flow associated with
the decay leads to a formation of a new feature within a void. This structure
may be interpreted as an additional ditch-like underdensity adjacent to the
edge of the void. The formation of this secondary structure can occur either
in the case of a short decay rate (i.e. the decay starts earlier in the past)
and low injection velocity or, alternatively, with a longer decay rate (i.e.
decay starts at a later stage of evolution) and larger injection velocity.
This is due to the fact that there is an interplay between the two parameters
$\Gamma$ and $v_{i}$ wherein the decrease in one may be compensated by the
increase of the other. This degeneracy follows from eq. (45), which shows that
the two parameters are both encapsulated in the effect of the scalar
acceleration term.

A selection of total density profiles is displayed in Figure 2 for an
uncompensated R-type profile and various parameters. In the case of the R-type
voids, apart from a very slight shallowing of the centre (in comparison to the
case without decay), there is no apparent impact on the mass distribution
inside the void. The reason why R-type voids do not exhibit distinctive
features is related to the fact the dark matter decays effects are second-
order, which is more clearly evident in Fig. 3 and 4. Initially (at the CMB
instant), the models considered in this paper, are very close to the FLRW
model, the evolution of which is driven by density and the expansion rate only
(cf. Sec. II.3). Cosmic voids are formed due to the growth of perturbations
imposed on the FLRW background Bolejko et al. (2005); Bolejko (2006, 2007).
Perturbations in the expansion rate and density are first-order whereas, as
seen from eqs. (17) – (21), both shear $\Sigma$ and Weyl curvature $W$ are
clearly second-order. Figures 3 and 4 show the present day form of the
expansion rate, shear, and Weyl curvature. The expansion rate is expressed in
terms of the Hubble parameter $H=\Theta/3$. For both types of voids (S-type
and R-type) the perturbation in the expansion rate is large but still very
similar to the LTB case (i.e. the inhomogeneous model but without presence of
heat flux). However, the largest difference between these models is seen in
the plots of shear $\Sigma$ and Weyl curvature $W$. In the case of the S-type
models the shear $\Sigma$ and Weyl curvature $W$ are clearly distinct and
large. In contrast, while for the R-type models the shear $\Sigma$ and Weyl
curvature $W$ are distinct from the LTB case they are now small, so their
contribution to the evolution of the system is slight. The results presented
in Figs. 3 and 4 indicate that the existence of the additional ditch-like
underdensity forming next to the edge of the void is a second-order effect
associated with decay-related heat flux. The process of forming this
substructure is dynamical and interrelated: as seen from eqs. (17) – (21), the
density evolution is affected by the divergence of the heat flux $Q$, which in
turn is influenced by shear $\Sigma$; the shear itself depends on $W$, which
is a function the divergence of heat flux. Thus if the gradients are small,
the growth of the heat flux and heat-flux related growth of shear and Weyl is
limited – this is the case of the R-type voids. On the other hand, in S-type
voids the gradients are large and their presence enhances the growth of the
shear and $W$, and leads to the formation of the additional ditch-like
underdensity.

Figure 6: Present day velocity of the second fluid, $v=Q/\eta$, for different
values of the injection velocity $v_{i}$ and $\Gamma=1/H_{0}$. For the S-type
void, the approximation of small velocity breaks for $v_{i}>7$ km/s.  Figure
7: Reconstruction of a density profile based on the number counts of tracer
within a void. The solid line shows the present-day density profile of an
S-type void with $\Gamma=1/H_{0}$, $v_{i}=7$ km/s. The decrease in density
near the edge of the void is $\Delta=0.16$. The number of points required to
obtain presented distribution is 1400 points. Assuming that there are
approximately 5 galaxies per void, this corresponds to stacking of
approximately 280 voids of similar size and spherical shape.

The extra sub-structure opens new avenues for using cosmic voids in studies of
the physics of DM. Detection of such a feature would suggest the presence of
heat flux associated with the decay of DM. In addition, one could target voids
at different redshifts to detect when the secondary structure starts is
initiated. Figure 5 shows a series of snap-shots of the evolution of an S-type
density profile at a few different redshift and different decay rates
$\Gamma$. As the models with shorter decay rate develop this new feature
sooner, this structure will be apparent at a higher redshift. Thus by
targeting voids at different redshift one could, in principle, attempt to
constrain various particle properties such as the half-life
$t_{1/2}=\Gamma^{-1}\ln(2)$, the average decay injection velocity $v_{i}$, or
even particle masses. However, the results presented in Figure 5, should only
be treated a proof-of-concept.

First of all, as mentioned above, there is a degeneracy between the decay rate
$\Gamma$ and injection velocity $v_{i}$. Figure 5 shows snapshots for
$v_{i}=18$ km/s. With a longer decay rate, one could also obtain the feature
at higher redshifts if the injection velocity is sufficiently large. Thus,
without additional constraints on properties of dark matter particles (such as
masses and energy released during the decay) one may not be able to
disentangle constraints on $\Gamma$ from $v_{i}$.

Secondly, and most importantly, the model considered in this paper is built on
a framework that assumes that the relative velocity of one fluid with respect
to another is small. This neglects contributions from pressure and pressure
gradients. Thus, the framework breaks once the velocity of one fluid with
respect to another is sufficiently large.

In the linear regime of small perturbations, the relative velocity of fluids
decays. It is only when second-order effects (such as shear and Weyl
curvature) become non-negligible that the relative velocity starts to increase
Gaspar et al. (2019). Note that in the cases considered here, the decay-
related heat flux is sensitive to the injection velocity. If the injection
velocity is too large, the growth of the heat flux enhances the growth of the
shear and $W$, which in turn leads to a rapid growth of the relative velocity
and to a failure of the framework implemented in this paper. This phenomenon
is illustrated in Fig. 6. If the injection velocity is larger than $v_{i}=7.0$
km/s the model breaks at the present day. If the injection velocity is larger
then the model fails earlier. For the examples presented in Fig. 5: the model
with $\Gamma=0.5H_{0}^{-1}$ and $v_{i}=18$ km/s fails approximately at
$z=0.8$, while the model with $\Gamma=H_{0}^{-1}$ and $v_{i}=18$ km/s stops at
$z=0.4$. While the results presented here open a new possibility of testing
dark matter physics, i.e. by observing this secondary structure at certain
red-shifts, the framework requires further development in order to
quantitatively allow observational data to constrain properties of dark matter
physics.

Such a structure could be observed either directly, by mapping the galaxy
distribution inside cosmic voids, or, indirectly, by investigating light
propagation through cosmic voids. Due to a limited number of galaxies inside
cosmic voids (a few counts per void), it is not possible to use the galaxy
counts within a single void to detect this structure. Rather, one requires a
large number of voids for stacking the signal. By stacking voids of similar
size and shape, one can increase the number of tracers inside a “stacked
void”. Such a procedure is presented in Fig. 7 which shows an estimate of the
(minimal) number of tracers needed to detect a decrease of number of galaxies
near the edge of a void.

The solid line in Fig. 7 shows the present-day density profile of an S-type
void with $\Gamma=1/H_{0}$, $v_{i}=7$ km/s. This density profile was used to
generate a random distribution of tracers within a galaxy. The points
represent the average density obtained by taking a ratio of galaxies $N$
within a given radius-bin to average number of galaxies $\bar{N}$ (i.e.
assuming galaxies are uniformly distributed, with number density $\bar{n}$,
the expected number of galaxies within the radius-bin is $\bar{N}=\bar{n}V$,
where $V$ is the volume of a particular distance-bin). The error bars
represent 95% confidence interval which was inferred from Monte Carlo
simulations: using the density profile (solid line) random number of galaxies
within a void was generated; the density then followed
$\rho/\rho_{0}=N/\bar{N}$; repeating the process large number of times and
excluding lower 2.5% and upper 2.5% estimates, resulted in the 95% confidence
interval. The required number of tracers for the case presented in Fig. 7 is
1400. So, assuming, 5 galaxies per void Pan et al. (2012); Sutter et al.
(2014), this means that at least approximately 300 voids of similar size and
spherical shape would be required to detect the decrease in density
$\Delta=0.16$.

The decrease $\Delta$ depends on the parameters of void and properties of DDM
(such as $\Gamma$ and $v_{i}$). The larger the decay rate or the smaller the
injection velocity, the smaller value of $\Delta$, and hence more voids would
be required to detect this feature. This is demonstrated in Fig. 8 which shows
that when $\Gamma>1.05H_{0}^{-1}$ and $v_{i}<7$ km/s the decrease is lower
than $\Delta=0.05$ and that requires at least a few tens of thousand of voids
(of similar size and shape). Given the fact that for larger redshift (for tens
of thousands of voids, one requires higher redshift to probe a large volume)
the amplitude of $\Delta$ is lower this shows that a direct method of
detecting this feature based on the number counts of galaxies within voids is
not effective. For example, the void catalog of SDSS DR12 BOSS111The Sloan
Digital Sky Survey, https://www.sdss.org/ consists of 10,643 voids Mao et al.
(2017), discovered between $z=0.2$ and $0.7$. While impressive, it still seems
insufficient to detect the discussed feature (indeed cf. Fig. 6 in Mao et al.
(2017)). A potentially better strategy for detecting this feature is based on
the indirect method of gravitational lensing Krause et al. (2013); Ema et al.
(2021); Hossen et al. (2021) which could also be combined with the Doppler
lensing method Bolejko et al. (2013); Bacon et al. (2014); Hossen et al.
(2021). The detectability issue will be explored in future work.

Figure 8:  The fractional decrease in density near the edge of the voids
$\Delta$ as a function of DDM parameters $\Gamma$ and $v_{i}$. Colors
correspond to a minimal number of voids required to detect such a feature. We
emphasise this should be treated as a simple estimate obtained using a direct
method based on the number counts of galaxies within cosmic voids. This was
done by stacking voids of the same size and shape. Given that real voids are
of different size and shape as well as are subject to intrinsic variations in
density profiles, the actual number of required voids is most likely much
higher.

##  IV Conclusion

The aim of this paper has been to investigate the effects of DDM on cosmic
voids. We focused on cosmic voids as they are pristine environments that are
less prone to contamination by complex baryonic physics than to dense
environments such as galaxy clusters. By targeting voids one hopes to
disentangle the observational signatures of dark matter from those of baryonic
matter. To this end, we have developed a novel method of modelling the decay
of dark matter. In this model the comoving CDM component decays into another
dark particle with a non-relativistic velocity, thereby producing a secondary
fluid which has non-comoving average spatial velocity induced by the velocity
of the produced daughter particles. The effect of this non-comoving secondary
fluid on the evolution of cosmological voids has been investigated.

Recent work, Gaspar et al. (2019); Tsagas (2010), has presented thorough
analyses of the effect of non-comoving components upon cosmological and
structure formation. The present article has attempted to further these works
by suggesting unstable DDM as a mechanism for producing non-comoving fluid
components. Rather than the non-comoving velocity distribution, we have
explored the effect of the DDM particle parameters on void evolution. The
results of this analysis suggest that for S-type voids there is the formation
of a significant secondary substructure induced by the non-comoving DM
component produced via the decay.

With further sophistication, the model proposed and presented here may offer a
contribution towards constraining the unknown particle properties of DM, in
conjunction with constraints determined by other means. Our method could be
used to further constrain or discover DDM signatures. For example, if no
secondary substructure is observed than one could rule out combination of
decay lifetimes $\Gamma<H_{0}^{-1}$ Gyr and injection velocities of $v_{i}>10$
km/s. Our methods will allow us to put tighter constraints on DDM. As a
preliminary investigation, our goal has been to proffer cosmological voids as
a possible means of indirect detection of DDM, with no attempts made to
constrain DDM particle candidates, nor to satisfy the possible thermodynamic
relations which would be required for physicality, Coley (1989). The present
paper gives a proof-of-concept that demonstrates that cosmological voids can
offer an indirect, observational means to constrain properties of DM.

## Acknowledgements

The authors wish to acknowledge and thank Andrew Bassom, Larry Forbes, and
Stephen Walters for their support and discussions. KB acknowledges support
from the Australian Research Council through the Future Fellowship
FT140101270.

## References

  * Mack (2014) K. J. Mack, Monthly Notices of the Royal Astronomical Society 439, 2728 (2014).
  * Craig (2013) N. Craig, arXiv e-prints arXiv:1309.0528 (2013), [arXiv:1309.0528].
  * Akerib et al. (2014) D. S. Akerib, H. M. Araújo, X. Bai, A. J. Bailey, J. Balajthy, S. Bedikian, E. Bernard, A. Bernstein, A. Bolozdynya, A. Bradley, et al., Phys. Rev. Lett. 112, 091303 (2014), [arXiv:1310.8214].
  * Dolgov and Zeldovich (1981) A. Dolgov and Y. B. Zeldovich, Reviews of Modern Physics 53, 1 (1981).
  * Fox et al. (2014) P. J. Fox, G. Jung, P. Sorensen, and N. Weiner, Phys. Rev. D 89, 103526 (2014), URL https://link.aps.org/doi/10.1103/PhysRevD.89.103526.
  * Boehm et al. (2002) C. Boehm, A. Riazuelo, S. H. Hansen, and R. Schaeffer, Physical Review D 66, 083505 (2002).
  * Schaeffer and Silk (1988) R. Schaeffer and J. Silk, The Astrophysical Journal 332, 1 (1988).
  * Wang et al. (2014) M.-Y. Wang, A. H. Peter, L. E. Strigari, A. R. Zentner, B. Arant, S. Garrison-Kimmel, and M. Rocha, Monthly Notices of the Royal Astronomical Society 445, 614 (2014).
  * Spergel and Steinhardt (2000) D. N. Spergel and P. J. Steinhardt, Physical Review Letters 84, 3760 (2000).
  * Arabadjis et al. (2002) J. Arabadjis, M. Bautz, and G. Garmire, The Astrophysical Journal 572, 66 (2002).
  * Boehm and Schaeffer (2005) C. Boehm and R. Schaeffer, Astronomy & Astrophysics 438, 419 (2005).
  * Craig and Davis (2001) M. W. Craig and M. Davis, New Astronomy 6, 425 (2001).
  * Slatyer and Wu (2017) T. R. Slatyer and C.-L. Wu, Physical Review D 95, 023010 (2017).
  * Mack (2014) K. J. Mack, M.N.R.A.S. 439, 2728 (2014), [arXiv:1309.7783].
  * Lacroix et al. (2016a) T. Lacroix, O. Macias, C. Gordon, P. Panci, C. BÅ‘hm, and J. Silk, Phys. Rev. D 93, 103004 (2016a), [arXiv:1512.01846].
  * Lacroix et al. (2016b) T. Lacroix, J. Silk, E. Moulin, and C. BÅ‘hm, Phys. Rev. D 94, 123008 (2016b), [arXiv:1603.05977].
  * Clark et al. (2018) H. A. Clark, P. Scott, R. Trotta, and G. F. Lewis, JCAP 2018, 060 (2018), [arXiv:1612.01539].
  * Ellis (1973) G. F. Ellis, Cargese lectures in physics 6, 1 (1973).
  * Ellis and Wainwright (1997) G. F. R. Ellis and J. Wainwright, eds., _Dynamical Systems in Cosmology_ (Cambridge University Press, 1997).
  * Ellis et al. (2012) G. F. Ellis, R. Maartens, and M. A. MacCallum, _Relativistic cosmology_ (Cambridge University Press, 2012).
  * Blackadder and Koushiappas (2016) G. Blackadder and S. M. Koushiappas, Physical Review D 93, 023510 (2016).
  * Pandey et al. (2020) K. L. Pandey, T. Karwal, and S. Das, Journal of Cosmology and Astroparticle Physics 2020, 026 (2020).
  * Ackerman et al. (2009) L. Ackerman, M. R. Buckley, S. M. Carroll, and M. Kamionkowski, Physical Review D 79, 023519 (2009).
  * Cen (2001) R. Cen, The Astrophysical Journal Letters 546, L77 (2001).
  * Vattis et al. (2019) K. Vattis, S. M. Koushiappas, and A. Loeb, Physical Review D 99, 121302 (2019).
  * Wilson et al. (2007) J. R. Wilson, G. J. Mathews, and G. M. Fuller, Physical Review D 75, 043521 (2007).
  * Peter et al. (2010a) A. Peter, C. Moody, and M. Kamionkowski, in _American Astronomical Society Meeting Abstracts# 215_ (2010a), vol. 215, pp. 375–01.
  * Peebles (2001) P. Peebles, The Astrophysical Journal 557, 495 (2001).
  * van de Weygaert (2014) R. van de Weygaert, Proceedings of the International Astronomical Union 11, 493 (2014).
  * Sheth and Van De Weygaert (2004) R. K. Sheth and R. Van De Weygaert, Monthly Notices of the Royal Astronomical Society 350, 517 (2004).
  * Ceccarelli et al. (2013) L. Ceccarelli, D. Paz, M. Lares, N. Padilla, and D. G. Lambas, Monthly Notices of the Royal Astronomical Society 434, 1435 (2013).
  * Icke (1984) V. Icke, Monthly Notices of the Royal Astronomical Society 206, 1P (1984).
  * Rezaei (2019) Z. Rezaei, Monthly Notices of the Royal Astronomical Society 487, 2614 (2019).
  * Rezaei (2020) Z. Rezaei, The Astrophysical Journal 902, 102 (2020).
  * Bolejko (2006) K. Bolejko, Monthly Notices of the Royal Astronomical Society 370, 924 (2006).
  * Ellis and Bruni (1989) G. F. Ellis and M. Bruni, Physical Review D 40, 1804 (1989).
  * Lemaître (1933) A. G. Lemaître, Annales de la Société Scientifique de Bruxelles 53 (1933).
  * Lemaître (1997) A. G. Lemaître, General Relativity and Gravitation 29, 641 (1997), ISSN 1572-9532, URL http://dx.doi.org/10.1023/A:1018855621348.
  * Tolman (1934) R. C. Tolman, Proc. Natl. Acad. Sci. U.S.A. 20, 169 (1934).
  * Tolman (1997) R. C. Tolman, General Relativity and Gravitation 29, 935 (1997), ISSN 0001-7701, URL http://dx.doi.org/10.1023/A%3A1018891418565.
  * Bondi (1947) H. Bondi, Mon. Not. R. Astron. Soc. 107, 410 (1947).
  * Krasiński (1997) A. Krasiński, _Inhomogeneous Cosmological Models_ (Cambridge University Press, 1997).
  * Kim et al. (2019) D. Y. Kim, A. N. Lasenby, and M. P. Hobson, Monthly Notices of the Royal Astronomical Society 488, 4081 (2019).
  * van Elst and Ellis (1996) H. van Elst and G. F. Ellis, Classical and Quantum Gravity 13, 1099 (1996).
  * Van Elst et al. (1997) H. Van Elst, C. Uggla, W. M. Lesame, G. F. Ellis, and R. Maartens, Classical and Quantum Gravity 14, 1151 (1997).
  * Coley (1990) A. A. Coley, General relativity and gravitation 22, 3 (1990).
  * Goode (1986) S. W. Goode, Classical and Quantum Gravity 3, 1247 (1986).
  * Sebastián and Sussman (2021) N. Sebastián and R. A. Sussman, The European Physical Journal. C, Particles and Fields. 81 (2021).
  * Maartens et al. (1999) R. Maartens, T. Gebbie, and G. F. Ellis, Physical Review D 59, 083506 (1999).
  * Ellis and Van Elst (1999) G. F. Ellis and H. Van Elst, in _Theoretical and Observational Cosmology_ (Springer, 1999), pp. 1–116.
  * Tsagas et al. (2008) C. G. Tsagas, A. Challinor, and R. Maartens, Physics Reports 465, 61 (2008).
  * King and Ellis (1973) A. R. King and G. F. R. Ellis, Communications in Mathematical Physics 31, 209 (1973).
  * Goode and Wainwright (1986) S. Goode and J. Wainwright, General relativity and gravitation 18, 315 (1986).
  * Singh et al. (2016) S. Singh, G. F. R. Ellis, R. Goswami, and S. D. Maharaj, Phys. Rev. D 94, 104040 (2016), [arXiv:1609.03288].
  * Singh et al. (2018) S. Singh, R. Goswami, and S. D. Maharaj, arXiv e-prints arXiv:1805.05664 (2018), [arXiv:1805.05664].
  * Bolejko (2017) K. Bolejko, Classical and Quantum Gravity 35, 024003 (2017).
  * de Lima and de Garcia Maia (1985) J. A. S. de Lima and M. R. de Garcia Maia, Physics Letters A 110, 366 (1985).
  * Nájera and Sussman (2020) S. Nájera and R. A. Sussman, arXiv preprint arXiv:2011.11192 (2020).
  * Planck Collaboration et al. (2016) Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, J. G. Bartlett, et al., Astron. Astroph. 594, A13 (2016), [arXiv:1502.01589].
  * Bolejko et al. (2009) K. Bolejko, A. Krasiński, C. Hellaby, and M.-N. Célérier, _Structures in the Universe by Exact Methods: Formation, Evolution, Interactions_ (Cambridge University Press, 2009).
  * Bolejko et al. (2011) K. Bolejko, M.-N. Célérier, and A. Krasiński, Classical and Quantum Gravity 28, 164002 (2011), [arXiv:1102.1449].
  * Sussman (2013) R. A. Sussman, Class. Quant. Grav. 30, 235001 (2013), [arXiv:1305.3683].
  * Peter et al. (2010b) A. H. Peter, C. E. Moody, and M. Kamionkowski, Physical Review D 81, 103501 (2010b).
  * Blackadder and Koushiappas (2014) G. Blackadder and S. M. Koushiappas, Physical Review D 90, 103527 (2014).
  * Wang and Zentner (2012) M.-Y. Wang and A. R. Zentner, Physical Review D 85, 043514 (2012).
  * Mapelli et al. (2006) M. Mapelli, A. Ferrara, and E. Pierpaoli, Monthly Notices of the Royal Astronomical Society 369, 1719 (2006).
  * Bolejko et al. (2005) K. Bolejko, A. Krasiński, and C. Hellaby, M.N.R.A.S. 362, 213 (2005), [arXiv:gr-qc/0411126].
  * Bolejko (2006) K. Bolejko, M.N.R.A.S. 370, 924 (2006), [arXiv:astro-ph/0503356].
  * Bolejko (2007) K. Bolejko, Phys. Rev. D 75, 043508 (2007), [arXiv:astro-ph/0610292].
  * Gaspar et al. (2019) I. D. Gaspar, J. C. Hidalgo, and R. A. Sussman, The European Physical Journal C 79, 106 (2019).
  * Tsagas (2010) C. G. Tsagas, Monthly Notices of the Royal Astronomical Society 405, 503 (2010).
  * Pan et al. (2012) D. C. Pan, M. S. Vogeley, F. Hoyle, Y.-Y. Choi, and C. Park, M.N.R.A.S. 421, 926 (2012), [arXiv:1103.4156].
  * Sutter et al. (2014) P. M. Sutter, G. Lavaux, B. D. Wandelt, D. H. Weinberg, M. S. Warren, and A. Pisani, M.N.R.A.S. 442, 3127 (2014), [arXiv:1310.7155].
  * Mao et al. (2017) Q. Mao, A. A. Berlind, R. J. Scherrer, M. C. Neyrinck, R. Scoccimarro, J. L. Tinker, C. K. McBride, D. P. Schneider, K. Pan, D. Bizyaev, et al., Astrophys. J. 835, 161 (2017), [arXiv:1602.02771].
  * Krause et al. (2013) E. Krause, T.-C. Chang, O. Doré, and K. Umetsu, Astrophys.J.Lett. 762, L20 (2013), [arXiv:1210.2446].
  * Ema et al. (2021) S. A. Ema, R. M. Hossen, K. Bolejko, and G. F. Lewis, accepted to MNRAS (2021), [arXiv:2110.14089].
  * Hossen et al. (2021) R. M. Hossen, S. A. Ema, K. Bolejko, and G. F. Lewis, submitted to MNRAS (2021).
  * Bolejko et al. (2013) K. Bolejko, C. Clarkson, R. Maartens, D. Bacon, N. Meures, and E. Beynon, Phys. Rev. Lett. 110, 021302 (2013), [arXiv:1209.3142].
  * Bacon et al. (2014) D. J. Bacon, S. Andrianomena, C. Clarkson, K. Bolejko, and R. Maartens, M.N.R.A.S. 443, 1900 (2014), [arXiv:1401.3694].
  * Coley (1989) A. Coley, Physics Letters A 137, 235 (1989).