Testing time-delayed cosmology

Motivated by the proposed time-delayed cosmology in the primordial inflationary era, we consider the application of the delayed Friedmann equation in the late-time Universe and explore some of its observable consequences. We study the background evolution predicted by the delayed Friedmann equation and determine the growth of Newtonian perturbations in this delayed background. We reveal smoking-gun imprints of time-delayed cosmology that can be traced to derivative discontinuities generic in delay differential equations. We show that a late-time cosmic delay is statistically consistent with Hubble expansion rate and growth data. Based on these observables, we compute a nonzero best estimate for the time delay parameter and find that the Bayesian evidence does not strongly rule out a late-time time delay but warrants the subject further study.


Introduction
The discovery of the cosmic microwave background stamped the Big Bang model as a canonical theory of cosmic evolution. But despite its success, the Big Bang model is afflicted by several problems. These include the flatness problem, which is the problem of why the Universe on very large scales is approximately spatially flat today and even flatter before [1]; the horizon problem, which is the problem of causally disconnected regions in the sky having the same temperature [2,3]; and the monopole problem, which refers to the absence of observed magnetic monopoles [4]. These problems are not pathologies of the Big Bang per se, but the model does not have the predictive power to solve them. Therefore, the need for an addendum to the Big Bang model is widely supported.
The mainstream resolution to these cosmic conundrums is inflation [5,6]. Inflation posits that the early Universe underwent exponential expansion. This accelerated expansion drove down the initial curvature of spacetime, locked in the uniformity of the Universe, and diluted the density of magnetic monopoles to negligible levels all at once. But despite the elegance of the theory, the usual implementation of inflation via scalar fields called inflatons comes with its own problems [7,8]. Inflaton models usually violate energy conditions [9,10], and the fundamental nature of inflatons also remains an open question [11,12]. These reasons continue to motivate the search for alternative mechanisms [13][14][15][16][17].
One such proposed mechanism is the time-delayed cosmology of Choudhury et al. [18]. In this proposal, the evolution of the energy density of the Universe as expressed in the Friedmann equation is delayed by a constant τ relative to expansion. While ad hoc and somewhat unnatural, the scheme does generate inflation without the aforementioned problems of inflation models. An exploration of its consequences is premised on the possibility of some nonlocal theories effectively generating time-delayed responses in gravitational dynamics [19][20][21], and on the richer dynamics afforded by time-delayed systems [22,23]. More broadly, it answers a general invitation to explore the potential role of delay differential equations in fundamental physics [24].
However, time-delayed cosmology has received scant attention from the community, perhaps largely owing to its detachment from fundamental theory. Not surprisingly, there are hardly any observational constraints on its key parameter, the time delay τ , which is generally just presumed to be of the order of the Planck time. One notable attempt was made in Ref. [25], where the matter power spectrum was calculated in an ad hoc general relativistic delay scheme. It also obtained an estimate for the delay parameter, though this was not based on the observational power spectrum data.
In this work, we adopt the stance that the time-delayed cosmology proposal can at least be empirically interesting and make initial steps towards filling the aforementioned constraint gap. In order to do this, we apply time-delayed cosmology to the late Universe where an optically invisible fluid, often dubbed dark energy, supersedes matter and radiation to source the observed late-time cosmic acceleration [26][27][28][29][30]. The substantial evidence for dark energy and the theoretical parallels between primordial inflation and dark energy make the application of time-delayed cosmology to the dark Universe today worth undertaking.
In particular, we consider the effects of the delayed Friedmann equation [Eq. (2.4)] at late times and determine the background evolution as well as the growth of Newtonian perturbations about this delayed background expansion. We calculate the Hubble expansion rate H(z) as well as two growth observables, the growth rate f (z) and f σ 8 (z). We show clear dependence of the predictions on the time delay parameter τ and estimate this parameter directly from observational data. This hints at the potential relevance of time-delayed cosmology in a cosmic era that has not been demonstrated until now.
In the next section, we briefly introduce important details of time-delayed cosmology. In Section 3, we discuss the background evolution. In Section 4, we set up the equations for the growth of matter perturbations and discuss the predictions of time-delayed cosmology. In Section 5, we perform a Markov chain Monte Carlo sampling and estimate the time delay parameter τ directly from the Hubble expansion rate H(z), the growth rate f (z), and f σ 8 (z) data. Finally, we conclude our work in Section 6. The code for reproducing the figures and calculations in this paper can be freely downloaded at [31].
Conventions. We work with the geometrized units c = 8πG = 1 and the mostly plus metric signature (−, +, +, +). A dot over a variable denotes differentiation with respect to the cosmic time t.

Time-delayed cosmology
We provide a short introduction to the foundations of time-delayed cosmology (Section 2.1) and discuss the method of steps for solving a delay differential equation (Section 2.2). We then describe the set-up of time-delayed cosmology at late times (Section 2.3).

Foundations
With the observational support for large-scale statistical homogeneity [32][33][34] and isotropy [35][36][37], the standard description of cosmic evolution is given by the following Friedmann equation: where a(t) is the scale factor which measures the expansion of the cosmos and ρ(t) is the energy density of the perfect fluid permeating the Universe. Assuming an equation of state of the form p(t) = ωρ(t) (where p is the fluid pressure and ω is a constant called the equation of state parameter) and solving the continuity equation, , (2.2) the Friedmann equation can be written as where ρ i,x is the initial energy density of the fluid denoted by x. A late Universe described by the Friedmann equation and filled with a mixture of cold (i.e. pressureless) dark matter (ω = 0) and a cosmological constant Λ (ω = −1) is referred to as the standard ΛCDM model. On the other hand, time-delayed cosmology is based on a delayed Friedmann equation where τ is some constant that, in the original application in the inflationary era, has units of Planck time t p ∼ O(10 −44 ) s. Delaying the energy density term in this way is completely ad hoc and not unique. A fundamental action is yet to be found to prescribe this time-delayed set-up. However, there are motivations for the delay. A nonlocal theory of (quantum) gravity may induce a delayed response on the Universe since nonlocality may imply time-smeared interactions [18]. For example, the Deser-Woodard models [38][39][40], which are inspired by quantum loop corrections, involve cosmological equations with retarded boundary conditions. A more explicit example is Ref. [41] in which nonlocality has resulted to equations of motion that are systems of delay differential equations. In the application to the late-time era, we take timedelay to be a phenomenological step to alternatively source cosmic acceleration. We shall not expound on the origins of the delay, rather we shall regard this delay as a phenomenological parameter.

The method of steps
The delayed Friedmann equation can be solved with the method of steps [22,23]. The essential idea of the method of steps is to replace the delayed term with a known solution a(t) so that the delay differential equation becomes ordinary within an interval that is then solvable with standard methods. Effectively, the solution to the delayed equation (as with any constant-delay differential equation) is a piecewise function with each composite solution defined on an interval of the size of the delay τ . The first composite solution, which is to be defined, is called an initial history. This is the equivalent of the initial value in ordinary differential equations. Figure 1 illustrates the method of steps.
For a power-law initial history φ(t) = t α defined on t ∈ [0, τ ), which may be due to a quantum gravity effect or a pre-Big Bang scenario, the delayed Friedmann equation admits inflation in the following interval: , a n,τ )

Get solution a n+1
Given DDE: · a = f(a, a τ ) t 0 + nτ t 0 + (n + 1)τ a τ := a(t − τ) Notation: a := a(t) Figure 1: The method of steps. Starting with a given delay differential equation, we define an initial history φ(t) on an interval at least the size of one delay unit, e.g. [t 0 , t 0 +τ ). On the succeeding interval [t 0 + τ, t 0 + 2τ ), we replace the delayed term with φ(t − τ ) and solve the ensuing ordinary differential equation, using φ(t 0 +τ ) as an initial value. We label the solution of this ordinary differential equation as a 1 (t). On the following interval [t 0 +2τ, t 0 +3τ ), again with a size of one delay unit, we replace the delayed term with a 1 (t − τ ) and again solve the ensuing ordinary differential equation for a 2 (t), with a 1 (t 0 +2τ ) as an initial value. We repeat this process for as many intervals as we like, using the previous solution a n (t) to replace the delayed term and solve the resulting ordinary differential equation for a n+1 (t). The piecewise function defined by a n (t)'s comprise the solution to the delay differential equation.
For succeeding times, the delayed equation has to be solved numerically. In this work, we use the ddeint Python package [42] for the numerical solutions. Clearly, in time-delayed cosmology, inflation can be naturally generated for a period of one delay unit. This period also seamlessly ends, thereby avoiding the "graceful exit" problem (see Figure 2).

Application to late times
At late times, the phenomenon of interest is cosmic acceleration due to dark energy. Because of the parallels between primordial inflation and late-time cosmic acceleration, the application of the delayed Friedmann equation at late times is worth considering. Furthermore, it will be easier to place constraints on time-delayed cosmology if we can show its impact on the expansion era.
In this application, we phenomenologically regard dark energy as a mixture of a cosmological constant and a time delay. The delayed Friedmann equation is therefore of the form where H(t) :=ȧ(t)/a(t) is the Hubble parameter, H i is the initial Hubble parameter value, Figure 3 shows that the delayed Friedmann equation can also accommodate a late-time cosmic acceleration.
Note that in the original implementation of the delay modification in Ref. [18], as well as in Ref. [25], the delay τ is assumed to have units of Planck time t p , being the relevant time scale for inflation. Because the delay is very small, it would have no impact on late-time observables, which is why the computed power spectrum in Ref. [25] expectedly appears to be in excellent agreement with observations. In this application, we allow the delay to be large since the relevant time scale t c is also large; we will integrate in the matter-dominated  i , with H i being the Hubble parameter in the matter-dominated era. We find the value of H i using the following relation for a constant dark energy density where H 0 and Ω 0,Λ are the Hubble constant and the present value of the cosmological constant density parameter, respectively. The latest Planck 2018 estimates are H 0 = 67.4 ± 0.5 km s −1 Mpc −1 and Ω 0,Λ = 0.6889 ± 0.0056 [30]. Solving for H i in Equation 2.8, we get Starting our integrations in the time when Ω i,Λ = 10 −6 (and Ω i,m = 1 − Ω i,Λ ), this results to a time scale of t c = 0.0175 ± 0.0001 Gyr, where here and throughout the paper we have used the SOAD package [43] for (asymmetric) error propagation. The units of the time delay will be in terms of this time scale t c .

Hubble expansion
To obtain the background evolution from Equation 2.7, we must specify an initial history. Throughout this paper, we assume a power-law initial history of the form φ(t) ∼ t α for the delayed Friedmann equation. We have checked numerically that the observables we are interested in in this paper do not strongly depend on the parameter α (see Figure 12) in redshifts that are currently accessible to us and especially for reasonable values of α (that is, for α ≈ 2/3 which refers to the canonical matter-era solution). Furthermore, although α gains a stronger effect at very large redshifts in terms of affecting the magnitude of the observables and for very large time delays (H i τ > 10), the general shape of the observables are determined by the time delay parameter and not by α. For these reasons and combined with the fact that an α = 2/3 constitutes a more natural initial history (taking after the canonical t 2/3 matter-era solution), we choose to fix the value of α to 2/3 in the following calculations instead of taking it as a free parameter. Due to Figure 3, we can already expect that the background evolution of time-delayed cosmology closely follows that of ΛCDM. Indeed, if we look at the predictions for the apparent magnitude m(z) of supernovae in Figure 4, we can see that time-delayed cosmological predictions are virtually indistinguishable from the ΛCDM prediction. This is the case even when considering delays on the order of H i τ ∼ 10 and when considering larger redshifts. The difference between ΛCDM and time-delayed cosmology is revealed when we look at the Hubble expansion rate H(z). Figure 5 shows the evolution of the Hubble expansion rate H(z) for a fixed H 0 . The dashed red curve shows the prediction of the standard ΛCDM model and the black curves are the predictions of time-delayed cosmology. Predictions due to delays that are larger than H i τ = 1 already notably deviate from the ΛCDM prediction at redshifts z > 5. Notice that the predictions appear to start at different redshifts. This is the case for all the redshift plots in this paper. This happens because different delays affect the scale factor evolution differently, which is then used to obtain the redshift. However, we have made sure that all quantities start out with the same initial condition at the same starting integration time.
A striking observation in Figure 5 are kinks (encircled in blue) or points at which the predictions of time-delayed cosmology change sharply. In fact, the derivative of H(z) at any of the kinks is undefined. This is an expected feature and an artifact of delay differential equation models. It is well known that discontinuities propagate in the derivatives of the solution to delay differential equations [22,23]. At the start of integration, the first derivative is discontinuous. One delay unit afterwards, the discontinuity propagates in the second derivative. Since our definition of the Hubble expansion rate involvesȧ(t), we expect to see the discontinuity in the second derivativeä(t) at a certain point (the kinks).
These kinks in the Hubble expansion rate already provide an upper bound on the time delay without further statistical analysis. In Figure 6  magnitude H i τ ≈ 19 can already be ruled out due to the presence of a kink that the data clearly does not accommodate. As the value of the time delay is increased, the kinks in the Hubble expansion rate are revealed at smaller and smaller redshifts. This is also true if we consider other observables. Therefore, all time delays H i τ > 19 are also ruled out.
Of course, we do not expect real, physical quantities of the Universe to exhibit these discontinuities. But if our Universe is correctly modeled by a delayed Friedmann equation, the abrupt transitions in H(z) serve as a generic smoking gun that make them empirically interesting. There may be fundamental reasons behind these discontinuities. For example, the improved Deser-Woodard model [38] was shown to have a discontinuous evolution of matter perturbation [44]. A possible cause of the discontinuity is a strong nonlocal effect.
For time delays with magnitude H i τ < 19, the smoking-gun imprints of time-delayed cosmology on the background evolution are only revealed at large redshifts for which data is still unavailable. When we look at redshifts z < 2 ( Figure 6), we can see that timedelayed cosmology closely follows ΛCDM even for delays on the order of H i τ ∼ 10. This shows that the key time delay parameter does not have to be of the order of Planck time as originally envisaged in order to fit observational data. Even large cosmic delays appear to be viable. An unfortunate consequence of this is that the Hubble expansion data is unable to distinguish time-delayed cosmology from ΛCDM. To observe the difference, we look to Newtonian perturbations.

Newtonian perturbations
In this work, we choose the growth of Newtonian matter perturbations as an additional probe of time-delayed cosmology. In particular, we are interested in the growth rate f (z) and another observable f σ 8 (z), where z is the redshift, that are both dependent on the amplitude of perturbations. The former quantity is the speed of growth of perturbations in the Universe with respect to the cosmic expansion, and the latter is essentially the growth rate scaled by the evolving root-mean-square of matter perturbations. These observational probes of large-scale structures have been used to distinguish between modified gravity theories and the standard ΛCDM model [46][47][48][49][50][51][52][53][54]. As we shall see, these will also be useful for obtaining constraints on time-delayed cosmology.

Set-up
The growth of Newtonian matter perturbations is given by [46,47] where δ(t) := δρ m (t)/ρ m (t) is the density contrast quantifying the inhomogeneity of the universe and ρ m (t) is the background energy density of (dark) matter. This fluctuation equation is valid for sub-horizon perturbations, i.e., λ/a(t) H(t) −1 , where λ is the comoving mode wavelength of the density perturbation. It is convenient to rewrite and solve this equation in terms of the scale factor a or the redshift z (using the relation z = (1/a) − 1). Once δ(a) is obtained, the two observables of interest can be easily calculated using the following definitions: In what follows, we solve the fluctuation equation where we have replaced the background matter energy density with its Hubble function equivalent using the (delayed) Friedmann equation. The observables of interest in time are In particular, the predictions for time-delayed cosmology decreases initially before increasing and eventually peaking. There are also kinks (encircled in blue) in the predictions.
then given by where t 0 denotes the present day. In addition to assuming an initial history of the form φ(t) ∼ t 2/3 for reasons we mentioned before, we also set the canonical a(t) ∼ δ(t) ∼ t 2/3 solution as an initial condition for the perturbation equation. This means that we integrate deep in the matter-dominated era up to the present dark energy-dominated era. In comparing the results with the ΛCDM model, we use the latest value of σ 8 given by Planck: σ 8 = 0.811 ± 0.006 [30]. We note that we are using the standard (i.e. non-delayed) perturbation equation here instead of a new delayed perturbation equation. Absent a fundamental action for timedelayed cosmology, this is unfortunately the best that one can do, short of proposing further ad hoc prescriptions about how the delay directly affects perturbations. Here, we adopt the conservative position that time-delayed cosmology manifests itself only through the cosmic expansion, i.e. the Hubble function. We find that this conservative modification is enough to see interesting consequences of time-delayed cosmology without having to develop an actionbased delayed perturbation theory. Figure 7 shows the plot of the growth rate f (z). The dashed red curve shows the prediction of the standard ΛCDM model, whereas the black curves are the predictions of time-delayed cosmology at different values of the time delay parameter τ . Immediately, we can see that time-delayed cosmology makes very different predictions for f (z) even for "intermediate" (i.e. H i τ ∼ 10) values of the time delay parameter. The growth rate of a delayed Universe dips in the matter-dominated era (i.e. z > 1) and then peaks later on before dark energy finally suppresses it for good (see Figure 9). On the other hand, if the delays are "small" (i.e. H i τ 1), then these dips and peaks are weak or not visible at all and the growth rate is virtually indistinguishable from the prediction of ΛCDM.

Growth rate f (z)
The characteristic decreasing of the growth rate predictions earlier on implies that, for a certain period, the delayed Universe was expanding faster than the perturbations were growing. We can see this clearly in Figure 8a. In Figure 8b, we can also see that the rate of expansion is initially greater than the rate of growth of the fluctuation. Combined together, these two scenarios suppress the growth rate at early times. But after some time, the growth rate predictions start to increase after decreasing. Notice that the transition to this increasing phase is very abrupt. Since our definition of the growth rate also involvesȧ(t), we expect to see kinks just as we saw in the background evolution. Interestingly, the growth rate becomes greater than unity at a certain point, implying that in the delayed Universe, perturbations will eventually grow faster than the Universe is expanding. This is also clear in Figures 8a  and 8b. Later on, however, dark energy starts to dominate and the growth rate is eventually driven down. Figure 9 shows a closer look at the growth rate up to redshift z ∼ 3. Here, we can see that the kinks in the growth rate evolution can also provide an upper bound. While the uncertainties of growth rate data at z > 2.5 are very large, it is safe to say that H i τ ≈ 18 is already very unlikely to be viable. On the other hand, time delays with magnitude H i τ 10 do appear viable. Figure 10 shows the plot of f σ 8 (z). Again, the standard ΛCDM prediction is shown in dashed red, and the black curves are the predictions of time-delayed cosmology at different values of the time delay parameter. Similar to the scenario with the growth rate, time-delayed cosmology models with intermediate time delay parameter values predict f σ 8 (z) evolutions that are different from ΛCDM. The f σ 8 (z) of the delayed Universe starts off smaller than but eventually surpasses the standard prediction. When the cosmological constant becomes more important than matter during the dark energy-dominated era, time-delayed cosmology and ΛCDM follow each other in similar evolutions. Naturally, we also find that smaller delays lead to f σ 8 (z) predictions that are indistinguishable from ΛCDM. As with the growth rate, we also get a kink because the definition of f σ 8 (t) also includesȧ(t). Notice that in Figure  10 the predictions start out at different magnitudes. Note that all calculations started out with the same initial condition. The differences in the initial value in these plots are due to the normalizing constant δ(t = t 0 ), which is of course different for different models. Figure 11 is a closer look at f σ 8 (z) up to redshift z = 2. In this case, we do not see any kink even for a time delay H i τ ≈ 18. However, it is clear that a time delay H i τ ≈ 18 is already unlikely since its predicted evolution already misses plenty of data points. Time delays H i τ 10 do however lead to predictions that are already distinct from ΛCDM while also appearing viable.

Delay estimate
We have already obtained a strict upper bound for the time delay but to achieve a best estimate, we confront our numerical solutions for H(z), f (z), and f σ 8   where P (p|d, m) (the posterior) is the probability distribution of the parameters p given d and m, P (d|p, m) (the likelihood) is the probability of getting the data d given p and m, P (p|m) (the prior) is the probability of the parameters p according to our prior beliefs, and finally P (d|m) (the evidence) is a normalizing constant that, as we shall see, is important for model comparison.
We consider a likelihood L given by   [53] and [55].

Parameter
Prior where N is the number of data points, µ obs i is an observational data point, µ th i is a predicted value, and σ i is the observational error. Aside from the time delay, we also take the Hubble constant H 0 and σ 8 as free parameters to be estimated. Again, we fix α to 2/3 because α has a weak effect on the observables at small redshifts (see Figure 12) and because this value of α gives the canonical matter-era solution.
Our priors are shown in Table 1. We intentionally choose priors defined over wide ranges so as to avoid inadvertently cutting the posterior short. We note that since our priors are uniform, our arbitrary cutoffs for the priors do not affect the value of the best estimates of the parameters so long as the priors include these best estimates in their ranges. Since each of our prior is defined over a wide range, the best estimate for a parameter is guaranteed to be within the prior for that parameter.
Note that we do not consider negative delays (i.e. H i τ < 0). A negative delay means that the delayed Friedmann equation is advanced in time rather than retarded. In which case, we must provide future information instead of an initial history. The solution then would be the past evolution. Since we are interested in the predictions of time-delayed cosmology in the late Universe, the time delay must be strictly positive.
We use the PyMultiNest [56] and GetDist [57] Python packages to sample the posteriors via MCMC and post-process the resulting MCMC chains. We consider the Hubble expansion rate data compiled in Ref. [45], the growth rate data compiled in Refs. [53] and [55], and f σ 8 (z) data compiled in Ref. [50]. In what follows, we choose to report the median estimate which is more robust to outliers as compared with the mean. We have checked however that the median estimates below are not too different from the mean estimates and overlap with them within 1σ. Figure 13 shows the posterior distributions for the time delay parameter τ and the Hubble constant H 0 using Hubble expansion rate data alone. The median estimate for the time delay is H i τ = 5.59 +4.89 −3.86 or τ = 0.098 +0.085 −0.068 Gyr. Meanwhile, the median estimate for the Hubble constant is H 0 = 72.02 +1.13 −1.07 km s −1 Mpc −1 . Notably, the credible interval of the time delay estimate is rather large and this is the case for all the data we consider in this work. This may be attributed to two things. Firstly, the uncertainties of the observational data points are themselves large. And secondly, from Figure 6, we can't expect a sharply peaked posterior with a narrow credible interval because the predictions of time-delayed cosmology for varying time-delays are very similar. Figure 14 shows the posterior distribution for the time delay parameter τ using the growth rate dataset alone. Upon sampling, we find that the median estimate for the time delay is H i τ = 6.10 +4.04 −4.10 or τ = 0.106 +0.072 −0.072 Gyr. The estimate for τ has notably increased and we also find that the mass of the posterior distribution has moved to a nonzero time delay. On the other hand, Figure 15 shows the posterior distributions for the time delay parameter τ and σ 8 using the f σ 8 (z) dataset. The median estimate for the time delay is H i τ = 8.58 +4.04 −5.30 or τ = 0.150 +0.068 −0.091 Gyr. The median estimate for σ 8 is 0.77 ± 0.01. Notice in this case that the estimate for the time delay has gotten much larger. Figure 16 shows the posterior distributions when we combine the growth rate and f σ 8 (z) datasets. We find that the median estimates for the parameters are H i τ = 7.26 +3. 35 −4.57 or τ = 0.125 +0.060 −0.080 Gyr and σ 8 = 0.77 ± 0.01. What these results show is that growth observables or perturbations consistently prefer nonzero values of the time delay parameter, especially f σ 8 data. We can see this not only in the median estimate but also in the mode.
To arrive at a best estimate, we combine the background and growth datasets. Figure  17 shows     the Bayes factor which is roughly the Bayesian equivalent of the p−value used for classical (frequentist) hypothesis testing. Given data d and two models m 1 and m 2 , the preference for m 1 over m 2 in light of d is quantified by the Bayes factor B 12 defined as which is simply the ratio of the evidence of m 1 to the evidence of m 2 . This definition assumes that both models are equally probable before accounting for the data. This is a fair assumption in our case since this is the first time that a time delay is even being considered in late-time cosmology and we do not have prior information whether time-delayed cosmology is preferred over ΛCDM.
Letting m 1 denote ΛCDM and m 2 denote time-delayed cosmology, we compute the Bayes factor based on evidences calculated using the Hubble expansion rate data alone (ln B 12 = 0.446 ± 0.111), the combined growth data (ln B 12 = 0.100 ± 0.092), and finally the combined background and growth data (ln B 12 = 0.175 ± 0.142). Following the criteria in Ref. [58], we find that regardless of the data considered, the Bayes factor indicates a statistical preference against time-delayed cosmology in favor of ΛCDM that is not worth more than a bare mention (an odds in favor of ΛCDM less than 3:1). In other words, no conclusion can be drawn as to which model is favored.

Conclusion
This paper has made initial steps towards confronting the predictions of time-delayed cosmology with data. We have applied the delayed Friedmann equation in the late-time Universe and chose the Hubble expansion rate H(z) and Newtonian matter perturbations as our observational probes. We obtained the predictions for the late-time background evolution and the growth data. In calculating the growth observables, we have used the standard perturbation equation and assumed that the effects of time-delayed cosmology enter through the background expansion only. We find that the conservative assumptions we have made are sufficient to reveal smoking-gun imprints of the phenomenological time delay. These imprints can be credited to the propagation of discontinuities inherent in the solutions of delay differential equations. This is the first time that the effects of these discontinuities have been demonstrated in this model.
We showed that for "intermediate" (i.e. τ ∼ 0.175 ± 0.001 Gyr) values of the time delay parameter, time-delayed cosmology already makes different predictions as compared to ΛCDM, especially when looking at growth observables. The difference can be seen at redshifts that are currently accessible to us. Our best estimate of the key time delay parameter is τ = 0.113 +0.060 −0.069 Gyr using the combined Hubble expansion rate and growth datasets. Our calculation shows that the key time delay parameter does not have to be in orders of Planck time as originally proposed in order to be consistent with observations. We also calculated the Bayes factor and find no conclusive evidence in favor of ΛCDM against time-delayed cosmology. To our knowledge, this study is the first systematic attempt to place a statistical and data-driven constraint on time-delayed cosmology as applied to late times.
While we are mainly interested in the late Universe, one can take our calculated value of the time delay and ask what it means for an inflationary time-delayed cosmology. The biggest implication of a time delay as large as our estimated value is that inflation will last for millions of years. This goes against the usual estimate of the period of inflation lasting for a tiny fraction of a second which is based on certain assumptions on initial conditions, e.g. inflation started at around 10 −36 s after the initial singularity with a large Hubble parameter. If we relax these assumptions and allow a pre-Big Bang or ekpyrotic scenario, then it is not immediately evident how such a long period of inflation would be troublesome. It is the number of e-folds, and not the length, of inflation that is important.
Our results show that time-delayed cosmology is not only interesting theoretically, but it can also hold up against the standard ΛCDM model when confronted with currently available background and growth data. Our work provides a data-driven motivation to further study this phenomenological model. Future large-scale structure surveys [59,60] and high redshift distance indicators such as proposed standardizable candles (quasars [61] and gamma ray bursts [62]) and standard sirens [63,64] can be expected to further constrain the time delay, if not rule it out completely should the kink inherent to a time-delayed solution not be observed. We leave the search for a fundamental action that leads to a time-delayed cosmology for future work.