Observational constraints on cosmological future singularities

In this work we consider a family of cosmological models featuring future singularities. This type of cosmological evolution is typical of dark energy models with an equation of state violating some of the standard energy conditions (e.g. the null energy condition). Such a kind of behavior, widely studied in the literature, may arise in cosmologies with phantom fields, theories of modified gravity or models with interacting dark matter/dark energy. We briefly review the physical consequences of these cosmological evolution regarding geodesic completeness and the divergence of tidal forces in order to emphasize under which circumstances the singularities in some cosmological quantities correspond to actual singular spacetimes. We then introduce several phenomenological parameterizations of the Hubble expansion rate to model different singularities existing in the literature and use SN Ia, BAO and H(z) data to constrain how far in the future the singularity needs to be (under some reasonable assumptions on the behavior of the Hubble factor). We show that, for our family of parameterizations, the lower bound for the singularity time cannot be smaller than about 1.2 times the age of the universe, what roughly speaking means ∼2.8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sim }2.8$$\end{document} Gyrs from the present time.


I. INTRODUCTION
The standard model of cosmology together with the inflationary paradigm provide an accurate description of the universe, although it requires the presence of three unknown ingredients, namely: Dark matter, dark energy and the inflaton field.The last two share the property of being introduced in order to support phases of accelerating expansion.Moreover, while the inflaton accounts for the first instants of life of our universe, dark energy should determine its final fate as the component that will eventually dominate.If dark energy turns out to be simply a cosmological constant, then we are doomed to an asymptotically de Sitter universe in the future.The situation is much more subtle when dynamical dark energy or modified gravity is brought in as possible explanations for the late time accelerated expansion (for a review about dark energy models, see [1]).In some cases, dark energy is ascribed to a so-called phantom fluid, i.e., a fluid satisfying ρ + p < 0 and, thus, violating the Null Energy Condition (NEC) [2].For a set of minimally coupled scalar fields, this condition implies the presence of, at least, a Laplacian instability in the inhomogeneous perturbations, although this can be resolved by allowing non-minimal couplings (see for instance [3]).Moreover, such kind of behavior can be also a consequence of a modification of General Relativity instead of a fluid with a non-standard equation of state [4].In any case, the phantom behavior may affect the background evolution giving rise to a future singularity occurring at a finite time where the scale factor diverges.Nevertheless, note that some models with violations of the null energy condition do not drive the universe to a singularity but to regular scenarios that may affect the local structures, known as little Rip, Pseudo-Rip and Little Sibling [5][6][7].
The described singular behaviour is actually shared by many dynamical dark energy models and modified gravity scenarios, where divergences in different cosmological parameters at a finite time can appear.The nature of the future singularities may differ among the different scenarios and they can be classified according to the cosmological parameters that diverge.An alternative way of classifying the future singularities is by means of the derivative of the scale factor that diverges.This classification is very useful because it helps understanding the severity of the different types of singularities (for a classification of cosmological singularities, see Ref. [8,9]).At this respect, it is worth reminding that a singular spacetime is characterized by the incompleteness of the geodesics [10].Since the geodesic equations are linear in the connection, it will contain, at most, first derivatives of the metric.Thus, the geodesics will be regular as long as the metric is continuous at the singularity.For a cosmological model, this will mean that the scale factor should remain finite at the singularity, even if divergences in the Hubble expansion rate or its derivatives are present.This type of behaviour has been recently used in [11] in order to replace the Big Bang singularity with a milder one that can be trespassed by the geodesics.
Another useful equation in order to characterize the strength of a singularity is the geodesics deviation equation.That equation essentially determines the tidal forces suffered by two infinitesimally close geodesics and it depends on the curvature of the spacetime.This means that tidal forces are sensitive to singularities which do not necessarily affect the completeness of the geodesics.
Again in a cosmological context, if the scale factor remains regular, but the Hubble rate diverges, it is possible to have a regular geodesic congruence with divergent tidal forces.Some criteria based on the behaviour of the Riemann tensor as we approach the singularity exist in the literature to decide whether the singularity is strong or weak, being the Tipler [12] and Krolak [13] conditions two widely used ones.
Regardless the physical consequences of having a future singularity at a finite time, a natural question to ask is how close a given type of singularity is to us [14].This is the analogous of asking about the age of the universe, determined by our distance to the original Big Bang singularity.Nevertheless, in the same way as we do not expect the Big Bang singularity to exist actually, but rather being regularized by some quantum effects, high curvature corrections to Einstein's gravity or even by varying physical constants [15], we do not expect the future singularities to be physical, at least the strongest types where physical quantities diverge [16].However, it will be useful to have some estimation on how close to us a given singularity can be and, therefore, have an idea of how far in the future we could extrapolate a model with a certain type of future singularity.It is important to notice that an effective equation of state for dark energy w < −1 is within the confidence regions of observational data [17] so the possibility of having a future singularity is plausible.Moreover, such models have also received attention because of some theoretical implications, since possible quantum effects close to the singularity become important.We know that General Relativity is to be regarded as an effective field theory whose strong coupling scale is, in the most optimistic scenario, at the Planck scale.Thus, knowing at which time the singularity is essentially reached will give us also an idea of until when we can keep using General Relativity as an effective field theory.
The purpose of the present work is precisely to draw such an estimation in a fairly model independent framework.An important difficulty arising here with respect to the Big Bang case is that, while in that case we have control on the different phases that the universe has gone trough from the initial singularity until today, for the future singularity we cannot know what the future phases will be.Thus, we need to make some assumptions to eventually determine how close the singularity can be.In order to achieve this, we will use some classes of phenomenological parameterizations for the Hubble expansion rate as proxies for a universe with a transition from a matter dominated era to a dark energy phase leading to a future singularity.We will then confront them to SN Ia, BAO and H(z) data to obtain the time of the singularity.Obviously, there could be transient phases that could delay the singularity, but this will not concern us since we are actually interested in obtaining a general lower bound for a future singularity.
The paper is organized as follows: section II is devoted to a brief review about future cosmological singularities.
In section III, the parameterizations of the Hubble rate which are analyzed in the paper are introduced.Then, the observational data used to fit the models is described in section IV.Finally, section V is devoted to the results discussions.

II. FUTURE COSMOLOGICAL SINGULARITIES
Assuming a homogeneous and isotropic universe at large scales, in compliance with the cosmological principle, the metric is given by the Friedmann-Lemaître-Robertson-Walker (FLRW) line element which is expressed as follows where we have assumed spatial flatness.Within General Relativity and assuming a perfect fluid as matter source, the gravitational equations can be written as Here ρ and p are the energy and pressure densities respectively of the perfect fluid, while H = ȧ a is the Hubble parameter.These equations are enough to describe the background cosmological evolution once the matter content of the universe is specified.In addition to these equations, the Bianchi identities allow to obtain the continuity equation ρ + 3H(1 + w)ρ = 0 with w ≡ p/ρ = −1 the equation of state (EoS) parameter.In cosmological scenarios based on non-standard fields, general fluids of modified gravity, several singularities have been found to appear at a future finite cosmic time.The different types of finite late-time singularities can be classified according to the divergent cosmological quantity at the singularity as follows (see Refs. [8,9], • Type I ("Big Rip singularity"): For t → t s , a → ∞ and ρ → ∞, |p| → ∞.Time-like geodesics are incomplete [14,18].
• Type IV ("Generalized Sudden singularity"): For t → t s , a → a s and ρ → ρ s , p → p s but higher derivatives of Hubble parameter diverge.They are weak singularities [21].
Type III singularities were shown to appear naturally in traditional vector-tensor theories of gravity [23], while type II singularities can appear in a novel class of vector field theories that arise in generalized Weyl geometries [24].In addition, there are other scenarios where no quantity diverges at a finite time but at infinity, namely the "Little Rip" [5], "Pseudo-Rip" [6] and "Little Sibling" [7].The above classification is useful since it groups together different models exhibiting a background evolution where some cosmological quantity meets a divergence in the future.The fact that some given quantities might have a divergence is usually regarded as a nondesirable feature to have in a regular spacetime.However, a regular spacetime is only defined in terms of its geodesic completeness.Thus, a spacetime with a curvature divergence can be regular as long as the geodesics can smoothly go through the divergence.Hence, cosmological models with some of the divergences in the above classification do not need to correspond to singular spacetimes and, consequently, singular future universes.In order to study whether the different singularities correspond to a geodesically incomplete spacetime we will consider the geodesic equations given by where λ is some affine parameter (proper time for instance for non-null geodesics) and Γ µ αβ are the corresponding Christoffel symbols.This equation already shows that it is the connection which determines the smoothness of the geodesics.In general, the solutions of the differential equations will be better behaved than the coefficients of the equations, so it is plausible to have a divergence in the connection with the geodesics remaining well-defined.It is also important to notice that the curvature contains derivatives of the connection and, therefore, there can be situations with curvature divergences, but where the connection (and consequently the geodesics) are perfectly regular.We will illustrate this below for some specific cases.The relevant case for the cosmological evolution is a spacetime described by the FLRW metric.In that case, the geodesic equations read1 These equations can be easily integrated.We start by rewriting the Hubble parameter in terms of the affine parameter as Then, we can rewrite Eq. ( 5) as d dλ a 2 dx i dλ = 0 (7) which can be immediately integrated to obtain with u i 0 some integration constants.We can then use this solution into Eq.( 4) to obtain dt dλ where C 0 is another integration constant.We thus see that the geodesics will be regular (with a well-defined tangent vector) as long as the scale factor remains regular.If the scale factor does not diverge and is nonvanishing (so the metric is regular) the 4-velocities of the geodesics remain regular and the spacetime will be said to be non-singular.If the scale factor diverges at some point, then the geodesics stop there and cannot go through it.As we have discussed above, it is important to notice that the geodesics are insensitive to divergences in the expansion rate H or its derivatives if they do not correspond to a singular behavior of the scale factor.This will be the case of the types II, III and IV singularities in the above classification where the scale factor remains finite while all the divergences only appear in its derivatives.So far we have discussed the singularities from the point of view of geodesic completeness.Another class of criteria that is useful to study the presence of a singular physical behaviour is the geodesic deviation equation, which allows to infer the potential existence of divergences in tidal forces.The corresponding equations depend on the Riemann tensor, which explicitly contains the Hubble expansion rate and its first time-derivative so that it is, in principle, sensitive to divergences that do not affect the geodesics themselves.Two common criteria to classify these divergences are the so-called strong curvature divergences in the Tipler and Krolak sense, which are respectively characterized by the following integrals: with u i the 4-velocity of the geodesic approaching the singularity.Here again we see that divergences in the curvature do not necessarily lead to a physical singularity because integrals of a given function are generally better behaved than the function itself.Thus, even if the spacetime contains a curvature divergence, it can remain regular according to the above criteria.The physical reason roots in the fact that the geodesic deviation equation measures the infinitesimal deviation, i.e., the tidal force between infinitesimally closed geodesics.However, extended physical objects have a finite physical volume and the above criteria precisely give the conditions for a finite volume to remain finite when going through the singularity.On the other hand, if the tidal forces are strong enough such that the volume shrinks to zero, the singularity is said to be strong.

III. THE MODELS
In this section we will describe the parameterizations that we will use for the subsequent confrontation to observational data.We emphasize that we intend to establish a general lower bound for the time of the future singularity t s .Since we are dealing with future singularities occurring at a finite proper time, it is reasonable to perform our parameterizations in terms of proper time.Moreover, as we have discussed, the severity of the different types of singularities is essentially determined by whether the scale factor or any of its time-derivatives presents a divergence.Therefore, the natural cosmological quantity to parameterize is the scale factor.However, for convenience when confronting to SN Ia and BAO, it will be more appropriate to parameterize the Hubble expansion rate directly.By doing this, we also avoid the ambiguity in the normalization of the scale factor.
As commented in the introduction, the main difficulty with respect to constraining the time of the Big Bang is that, while we have an accurate knowledge about the past history of the universe so we can robustly compute such a time, the future evolution of the universe is completely unknown.Because of that we need to make some relatively strong assumptions on our parameterizations.First of all, we want to have an approximate matter dominated phase at early times; we will use low-redshift (z ≤ 2) cosmological data so that by early time we actually mean well inside the matter domination epoch, but much later than equality and decoupling times, i.e., redshifts 10 ≤ z ≤ 1000.In order to comply with this requirements we propose to use the following form for the Hubble expansion rate: where t s is the time when one of the above singularities occur and the function F (t, t s ) is assumed to be negligible for t ≪ t 0 with t 0 the present time, such that H(t ≪ t 0 ) ≃2 3t as it corresponds for a matter dominated universe.In terms of the scale factor, this translates into a parameterization of the form This matter dominated phase will then be matched to an evolution with a future time singularity, i.e., F (t, t s ) is a function that either itself or some of its time-derivatives diverges at t = t s .We will assume that this divergence originates from the fact that the differential equation of the underlying physical model presents a regular singular point so that the solution near the singularity can be expressed as a Frobenius series.We will further assume that the transition from the matter era is sufficiently fast so that the dominant term of the series rapidly takes over.This will not affect our goal of obtaining a lower bound for the time of the transition since making the transition slower typically delays the appearance of the divergence.Since we are looking for a future divergence where a given derivative of the scale factor diverges while the lower derivatives remain finite, a reasonable Ansatz for F (t, t s ) is some half-integer power.With these considerations in mind, we have chosen the specific parameterizations summarized in Table .I together with their main properties 2 .All the models contain two parameters characterizing the time of the singularity t s and an additional parameter n that regulates the time of the transition from matter domination.Notice that all the parameterizations share the property of containing a latetime de Sitter evolution when the time of the singularity is sent to the asymptotic future3 t s → ∞.However, it is important to notice that the existence of a matter phase at early times matching a de Sitter universe in the asymptotic future does not necessarily mean that the evolution mimics that of a ΛCDM model, because the transition era between the two phases may be completely different.
In fact, it is not difficult to see that none of our parameterizations contains ΛCDM within its parameter space.Although our parameterizations do not rely on a specific gravitational theory, in order to make contact with previous literature and give a physical intuition of what kind of theoretical models might be described by our parameterizations, we will now consider some explicit cases.If we assume General Relativity for the gravitational interaction, we can define an effective equation of state for the universe as Thus, we can interpret our parameterizations in terms of an effective equation of state parameter for the content of the universe, assuming a perfect fluid form and a barotropic equation of state.Since at early times our parameterizations give H ≃ 2/(3t), we recover a matter dominated universe with w eff ≃ 0 as it should.
On the other hand, the models considered in this manuscript can alternatively be interpreted as the result of the interaction between dust and dark energy.In such a case, the continuity equations for both fluids would be coupled, and could be expressed by [8], Here Q(t) accounts for the energy exchange between both fluids.By combining equations (15) with the Friedmann equation, the following expression for Q(t) is obtained: Hence, by assuming a particular Ansatz H = H(t), the corresponding interacting term Q(t) is obtained.Moreover, note that other type of non-interacting models can lead to the concerned models of Table I by means of non-standard EoS's for dark energy that effectively stand for modifications of the Hilbert-Einstein action, viscosity terms (see Ref. [27]).Nevertheless, the later case may lead to negative energy densities and other undesirable consequences, as shown below for our parameterizations of the Hubble parameter.Finally, we can also mention that a given background evolution for the scale factor can be mapped onto a scalar field theory by suitable choice of the action.In the spirit of the effective field theory of dark energy, we can think of the time coordinate as corresponding to a foliation of the spacetime according to the scalar field, where the unitary gauge has been chosen.Then, a natural interpretation of the future singularity would be a point where the scalar field meets a pathology in its evolution as dictated by the field equations.

IV. DATA
The analysis has been performed using three different standard cosmological tools.They are at low redshift (z 2), because we are not interested in changing early time evolution and we assume that a possible signature for "future" evolution toward a singularity, if any, is detectable now or, at least, in the recent past only.
Just for sake of clarity and computational motivations, all the models we propose are written in terms of the dimensionless variable x = t/t 0 , where t 0 is the age of the Universe.This means that all quantities will be measured in units of t 0 so, for instance, we will have Therefore, in this case t 0 plays, in terms of fitting parameters, the role usually ascribed to the Hubble constant H 0 in the standard approach.Since our parameterizations are explicitly expressed in terms of time, it will be more convenient to use all the standard integrals involved in the calculation of cosmological distances directly expressed as integrations over time, instead of transforming them into integrations over redshift, being the two of them related as The integrations over redshift are more convenient in the usual case because the observational data are given in terms of redshift.Thus, we will need to find the values of x that correspond to the given redshifts, i.e., we need to find the functions z = z(x) or, equivalently a = a(x).This can be easily obtained from the corresponding expression for H(x) by solving the differential equation where the prime stands for derivative with respect to x.This equation will be solved with the boundary condition a(x = 1) = 1, i.e., we normalize the scale factor to be 1 today.Thus we operationally define the time t 0 in our models by such condition.Notice that for all our parameterizations this can be done analytically.Therefore, we can obtain the values of x i corresponding to the measured values z i by numerically solving the equation A. Hubble data from early-type galaxies We use the compilation of Hubble parameter measurements estimated with the differential evolution of passively evolving early-type galaxies as cosmic chronometers, in the redshift range 0 < z < 1.97 and recently updated in [28].The corresponding χ 2 H estimator is defined as with σ H (z i ) the observational errors on the measured H obs (z i ) values, and θ is the vector of cosmological parameters, i.e., (t 0 , n, t s ) in our case.Moreover, we will add a gaussian prior, derived from the Hubble constant value given in [29], H 0 = 69.6 ± 0.7.Notice that now H 0 is a derived quantity depending on the actual fitting parameters so where the numerator H(x = 1, θ) now depends on the parameters n and x s .
as Hs > 0 Ḣs < 0 ∞ ρs 0 0 Table I.In this table we summarize the 5 parameterizations that we propose to describe the different types of future singularities that we consider throughout this work.In the first column we give the label we will use for each case, while the second column indicates the type of singularity according to the classification in [8].In the columns 3 and 4 we give the analytical expressions for H(t) and a(t) (where, as explained in the main text, the normalization a0 must be chosen so that a(x = 1) = 1).In the last columns we give the behaviour of a, H and some of its derivatives at the singularity.We also give the values of ρ, p and w eff for the theoretical interpretation discussed in Section III.It is important to keep in mind that those values depend on the underlying theoretical model and we only give them here for illustrative purposes.

B. Type Ia Supernovae
We use the SN Ia data from the Union2.1 compilation [30].The χ 2 SN in this case is generally defined as with ∆F SN = µ theo − µ obs the difference between the observed and theoretical value of the distance modulus µ, the observable quantity for Union2.1 SN Ia, defined as: with d L (z) the dimensionless luminosity distance given by where E(z) = H(z)/H 0 is the dimensionless Hubble function; µ 0 a nuisance parameter combining the Hubble constant H 0 (or t 0 in our case) and the absolute magnitude of a fiducial SN Ia.As usual, we marginalize the χ 2 SN over µ 0 .Finally, C is the covariance matrix.In terms of integration over time, the dimensionless luminosity distance can be expressed as:

C. Baryon Acoustic Oscillations
We have also made use of baryon acoustic oscillations (BAO), in particular, the data collected in [31,32] 4 .In this case the χ 2 BAO is defined as where, as before, ∆F BAO = F theo − F obs is the difference between the observed and theoretical value of the Alcock-Paczynski distortion parameter measured in a BAO survey, and defined as: with c the speed of light, H(z) the Hubble function, and D A the angular diameter given by: Even in a standard scenario, the quantity F (z) is independent of the parameter H 0 and can be written which in our notation translates into which is independent of the parameter t 0 .One important remark about BAO data concerns the possibility to use other than the Alcock-Paczynski variables like the angular diameter distance or other conveniently defined quantities (see, for example, [31,32]).Those quantities involve the calculation of the sound horizon at dragging epoch, which in turn requires knowledge about the density parameters of baryons and radiation.However, as explained in previous sections, our models only provide phenomenological parameterizations of the Hubble function in terms of cosmic time.For that we introduce some parameters whose relation to the more physical density id.Table II.In this table we present the obtained results for the best fit of each parameterization.In column 1 we give the label identifying each parameterization in Table I.In columns 2-5 we give the 1σ confidence levels for our primary model parameters (notice that they are different for the different cases and, in particular, Ωm is not within the fitting parameters of our parameterizations).In column 6 we show the age of the Universe.We also show the effective equation of state parameter (as defined in ( 14)) for each parameterization evaluated at the present.Finally, in columns 8 and 9 we give the Bayesian evidence and ratio with respect to ΛCDM for Jeffreys' interpretation.
parameters can only be inferred after assuming a particular theory of gravity.Thus, as a very conservative choice, we have decided to perform our analysis using only the Alcock-Paczynski method.For sake of correctness, we have also to stress that even the Alcock-Paczynski quantities are derived making some assumptions related to the background cosmology (at least, at a very early stage of raw data analysis), as it is pointed out in [33].But, in the same reference, it is claimed that their final results are not very sensitive to the fiducial model.Finally, the total χ 2 to be minimized will be χ 2 = χ 2 H +χ 2 H0 +χ 2 SN +χ 2 BAO .We minimize the total χ 2 using the Markov Chain Monte Carlo (MCMC) method and we check its convergence with the method developed in [34].In order to compare the models in the best statistical way possible, we have calculated the Bayesian evidence for each of them.The Bayesian evidence is defined as the probability of the data D given the model M with a set of parameters θ, E(M ) = dθL(D|θ, M )π(θ|M ): π(θ|M ) is the prior on the set of parameters, normalized to unity, and L(D|θ, M ) is the likelihood function.
We have been very careful in imposing priors; our parameters are, basically, t 0 , x s and n.For numerical convenience, we have used the parameter α instead of x s defined as in order to compactify the range x s ∈ (1, ∞) into α ∈ (0, 1).Since we are compactifying the range of our parameter x s from an infinite range to a finite one, we face the question of imposing an appropriate prior.For the parameter α, we have imposed two different priors on this range: a flat uniform prior, and a logarithmic one.This is due to the relation between α and the variable we are really interested in, x s , so that we want to have a full sampling of the very low range for α, which maps into x s → ∞ and, thus, it corresponds to nearly non-singular scenarios (i.e., the singularity is pushed to the far infinity).For n (and t 0 ) we assume a flat prior for only positive values, n > 0, given that this is the condition to ensure present accelerated expansion for all the models.The only exception is for the singularity D, where acceleration is guaranteed for n > n(α) > 0, with n(α) numerically found imposing the condition q(t) = 0, with q(t) being the deceleration parameter.Thus, the parameters span sufficiently wide and general ranges in order to have the same weight for each model when calculating the Bayesian evidence.The evidence is estimated using the algorithm in [35]; in order to reduce the statistical noise we run the algorithm many times obtaining a distribution of ∼ 100 values from which we extract the best value of the evidence as the mean of such distribution.
In order to compare the goodness of the different parameterizations, we further calculate the Bayes Factor, defined as the ratio of evidences of two models, M i and given the data.We will use the ΛCDM model as the reference model j (we have performed a further analysis with this model using the same data sets we have described above).The Bayesian evidence may be interpreted using Jeffreys' Scale [36], which tries to quantify the preference of a model against another based on the value of the evidence.In particular, if the ln B ij < 1, the evidence in favor of the highest-evidence model is not significant; if 1 < ln B ij < 2.5, the evidence is substantial; if 2.5 < ln B ij < 5, the evidence is strong; and if ln B ij > 5, the evidence is decisive.In [37], it is shown that the Jeffreys' scale is not a fully-reliable tool for model comparison, but at the same time the statistical validity of the Bayes factor as an efficient model-comparison tool is not questioned: a Bayes factor B ij > 1 unequivocally states that the model i is more likely than model j.We present results in both contexts for reader's interpretation.

V. RESULTS
After the datasets introduced in the previous section and the discussed considerations, we have proceeded to run the MCMC chains in order to obtain the confidence regions of each parameterization and, therefore, achieving the main goal of this work, namely, obtaining a lower bound for the time of a future singularity.The results corresponding to our different cases are shown in Fig. 1, where we display the marginalized contours of the parameters for each model, and in Table II.In order to have a further criterium to judge the statistical validity of our models, we have also analyzed, using the same data sets we have described in the previous section, the ΛCDM model and the Chevallier-Polarski-Linder (CPL) parametrization [38], which is widely used as the most basic generalization of a constant dark energy to a dynamical fluid.
When considering the combination of all the datasets, as can be visually checked from Fig. 1, we obtain that the lowest value for the singularity time is achieved for model B and turns out to be t s,min ≃ 1.2 (at the 2σ level), which corresponds to 2.8 Gyrs from today.Remarkably, this lower bound is prior-independent, i.e., it is for both the flat and the logarithmic priors on α that we have used; this helps us to state that such limit is not statisticallybiased, but physically compelling.In that regard, it is advantageous to report that the minimum in the χ 2 for models A and B is located, respectively, at x s ∼ 2.18 and ∼ 1.55.This does not happen for all other models, whose minima are located at α 0, and are only limited by numerical resolution at ∼ 10 −4 .Interestingly, the obtained lower limit for the occurrence of the singularity is shorter than the expected time for the Sun to burn all its fuel (estimated to be 5-7 Gyrs).
An interesting feature of models A and B is that having the singularity at infinity is excluded at the 1σ level, when a uniform prior is assumed; when using a logarithmic prior, model A only lies within the 1σ region, while model B still excludes the singularity at infinity at the 1σ level.We should remember that t s = ∞ corresponds to having a de Sitter universe in the asymptotic future, so for those two models, such a scenario seems to be disfavoured.This highlights that having an asymptotically de Sitter universe in our parameterizations does not necessarily implies being close to a ΛCDM model.And we also have to point out that when we use a logarithmic prior, which is going to give a better sampling than the uniform prior in the range of very small α, such higher bound disappears for model A, but not for model B. Such feature makes this latter model the most interesting among all those we have considered.Remarkably, these two models present a Bayesian evidence which make them equivalent to ΛCDM from a statistical point of view, i.e., they provide fits as good as those of ΛCDM, and they are also even better than the widely used CPL parameterization.Notice moreover that the effective equation of state parameter today is close to the one of ΛCDM.For a more direct comparison of our models with ΛCDM, in Fig. 2 we plot the effective equation of state for a total fluid as introduced in Eq. ( 14) which influences the background dynamics of the universe and the expansion history H(a).However, we need to note that this only happens for the restricted dataset considered in our analysis, while ΛCDM give a good fit to a much wider variety of cosmological observations, while it is not clear whether the models with singularities will fit all those observations as well as ΛCDM.
For model C, the possibility of having the singularity at infinity is within the 1σ region.The Bayesian evidence in this case is worse than that for models A and B, but still not strongly disfavoured with respect to ΛCDM.Interestingly the effective equation of state parameter today for this case is substantially lower than for ΛCDM.Models D and E are strongly disfavoured with respect to the baseline ΛCDM.Again, these models allow to have t s = ∞ at the 1σ level.In these cases, we find that w eff,0   is higher than in the ΛCDM case.Finally, in Fig. 3, we plot the interaction term given by Eq. ( 16), assuming dark energy equation of state equal to −1 and a dynamical one given by the CPL best fit we have found in our analysis, and reported in .In this plot we show the interaction Q (normalized to units of H 3 0 /(8πG) ) derived from our singularity models best fit and assuming a constant dark energy equation of state, w = −1.Singularity A: black -Singularity B: blue -Singularity C: magenta -Singularity D: red -Singularity E: green.Left: all models for all times -Right: zoom of the best models in the actual Universe time life range.We used solid (dashed) for positive (negative) Q. by using SN Ia, BAO and H(z) data.We have briefly reviewed the cosmological singularities emphasizing the fact that a divergence in a given cosmological parameter does not necessarily implies a singular spacetime.We have then discussed under which conditions a given cosmological singularity actually corresponds to a singular spacetime so that we can discern the severity of the different cosmological singularities.Our discussion focused on the geodesic completeness of the spacetime as well as the presence of divergent tidal forces when approaching the singularity.
After this brief theoretical review, we have constructed a set of parameterizations comprising different types of singularities.These parameterizations have been designed so that we recover an early time matter dominated phase that transits to a phase with a future singularity where a given time-derivative of the scale factor diverges, but not the lower ones.We have then run a series of MCMC chains to confront our parameterizations to SN Ia, BAO and H(z) data.The obtained results are then summarized in Table.II.Our main conclusion is that within our family of parameterizations, a potential future singularity cannot be closer to the present time than ∼ 0.2t 0 , that roughly corresponds to 2.8 Gyr.We found that the proximity of the singularity to the present time has a mild dependence on the type of singularity for our parameterizations, but we can conclude that in all cases there is a consistent lower bound around 1.2 − 1.5t 0 .
Another interesting conclusion that we have found is that, following results from the Bayesian evidence, our parameterizations A and B provide fits which are not significantly worse than ΛCDM for the considered datasets.This was not obvious a priori since none of our parameterizations contain ΛCDM in its parameter space.Hence, as shown in previous references [17], a singular scenario can be discarded right away from tests of the background evolution and the time remaining for the occurrence of a future singularity may be shorter than expected.However, we need to stress that ΛCDM has become the standard model of cosmology because of its outstanding performance in fitting most cosmological observations, not only the ones considered in our analysis, so that in order to be able to establish a compelling scenario with a future singularity on equal footing as ΛCDM, we would need to show its ability to fit the rest of cosmological observations, including those sensitive to the perturbations.

Figure 2 .
Figure 2. (Top Panel:) Effective EoS for the combined effect of matter and "singularity-fluid" when a uniform prior is applied.(Center Panel:) Effective EoS for the combined effect of matter and "singularity-fluid" when a logarithmic prior is applied.(Bottom Panel:) Rate expansion in terms of redshift, H(z), in the approximate range covered by data.ΛCDM from Table: solid light grey -CPL from Table: dashed light grey -Singularity A: black -Singularity B: blue -Singularity C: magenta -Singularity D: red -Singularity E: green.(Left: all models for all times -Right: zoom of the best models in the approximate time range covered by data).
Table: solid light grey -CPL from Table: dashed light grey -Singularity A: black -Singularity B: blue -Singularity C: magenta -Singularity D: red -Singularity E: green.(Left: all models for all times -Right: zoom of the best models in the approximate time range covered by data).
Table.II.In this work we have reconsidered the subject of future cosmological singularities occurring at a finite time.The aim of the work has been to establish a general lower bound for the time of a potential future singularity