Recent Advances in Cosmological Singularities

Abstract

The discovery of the Universe’s late-time acceleration and dark energy has led to a great deal of research into cosmological singularities, and in this brief review, we discuss all the prominent developments in this field for the best part of the last two decades. We discuss the fundamentals of spacetime singularities, after which we discuss in detail all the different forms of cosmological singularities that have been discovered in recent times. We then address methods and techniques to avoid or moderate these singularities in various theories and discuss how these singularities can also occur in non-conventional cosmologies. We then discuss a useful dynamical systems approach to deal with these singularities and finish up with some outlooks for the field. We hope that this work serves as a good resource to anyone who wants to update themselves with the developments in this very exciting area.

1. Introduction

Observations of the late-time acceleration of the Universe came as a huge surprise to the cosmological community [1], and ever since then a lot of work has been carried out in order to explain this expansion. The cosmological expansion problem has been addressed from multiple facets, which include the standard approaches of the cosmological constant [2,3,4] alongside more exotic scenarios like modified gravity theories [5,6,7] and scalar-field-driven late-time cosmic acceleration scenarios [8,9,10,11,12,13]. Several approaches to quantum gravity have also weighed in on the cosmic acceleration puzzle, ranging from the braneworld cosmology of string theory to the likes of loop quantum cosmology and asymptotically safe cosmology [14,15,16,17,18,19,20,21,22,23,24]. This, however, has also revealed some discrepancies that seem to be pointing towards the limits of our current understanding of the Universe, most famous of which is arguably the Hubble tension, referring to the disagreements between the values of the Hubble constant measured from detailed Cosmic Microwave Background (CMB) radiation maps combined with Baryon Acoustic Oscillation data and those from SNeIa data [25,26,27]. Hence, the current epoch of the Universe has certainly provided us with a wide range of questions and looks set to become an avenue whereby advanced gravitational physics will lead the way towards a better understanding of cosmology.

An expansive body of literature has also been published in recent times that has been devoted to the study of various types of singularities that could occur in the present and far future of the Universe, with the observation of late-time acceleration having given a significant boost to such works [28,29,30,31,32,33,34,35,36,37,38]. Even the term singularity comprises many different definitions. With regards to cosmological cases, until the end of the 20th century, the only popular possibilities for singularity formation were the initial Big Bang singularity and, in the case of spatially closed cosmological models, the final Big Crunch singularity. The definition of a singular point in cosmology was given by Hawking and Penrose, and most of the theorems they proved make use of the null energy condition and the fact that at a singular point in spacetime, geodesic incompleteness occurs and the curvature scalars diverge. Although in modified gravity the null energy condition may be different in general compared to the Einstein–Hilbert case (see, for example, [6,7]), it is generally accepted that the geodesic incompleteness and the divergence of the curvature invariants strongly indicate the presence of a crushing singularity. The singularities in cosmology vary in their effects, and a complete classification of these singularities was performed in [31,39]. While one can treat singularities as points at which a cosmological theory somewhat fails, one might also consider them as windows to new physics, thus providing a different kind of appeal. In particular, finite-time singularities (those which happen in a finite time) could be viewed as either flaws in the classical theory or, alternatively, a doorway towards a quantum description of general relativity. This is due to the fact that these cannot be addressed in a similar way to the spacelike singularities of black holes, for instance, and so one is left to ponder the accuracy of the predictions of classical gravitational theories. Hence, studying singularities in cosmological contexts and how they could (possibly) be removed provides a route towards a deeper understanding of the relationship between quantum descriptions of cosmology and classical ones.

These cosmological singularities that have been discussed in recent times can be classified broadly into two types: strong and weak (such a classification was initially put forward by [40]). Strong singularities are those singularities that can distort finite objects and can mark either the beginning or the end of the Universe, with the Big Bang being the one for the start of the Universe and the so-called “Big Rip” signaling the end of the Universe. Weak singularities, as the name might suggest, are those that do not have such far-reaching implications and do not distort finite objects in the same sense as their strong counterparts. We discuss these various singularities in more detail as follows, in accordance with the classification provided in [31,39]:

• Type $-1$ (“Grand Bang/Grand Rip”): In this case, the scale factor becomes null (bang) or diverges (rip) for $w=-1$ [41].
• Type 0 (“Big Bang”): In this case, the scale factor becomes null for $w\ne -1$.
• Type I (“Big Rip”): In this case, the scale factor effective energy density and effective pressure density diverge for $w\ne -1$. This results in a scenario of universal death, where everything within the Universe undergoes progressive disintegration [42].
• Type II (“sudden/quiescent singularity”): In this case, the pressure density diverges and so do the derivatives of the scalar factor from the second derivative onwards [43]. It is also known as a quiescent singularity, but this name originally appeared in contexts related to non-oscillatory singularities [44]. A special case of this is the Big Brake singularity [45].
• Type III (“Big Freeze”): In this case, the derivative of the scale factor from the first derivative onwards diverges. This was detected in generalized Chaplygin gas models [46].
• Type IV (“generalized sudden singularities”): These are finite-time singularities with finite density and pressure instead of diverging pressure. In this case, the derivative of the scale factor diverges from a derivative higher than the second [31,47].
• Type V (“w-singularities”): In this case, the scale factor and the energy and pressure densities are all finite, but the barotropic index $w=\frac{p}{\rho }$ becomes singular [48].
• Type ∞ (“directional singularities”): Curvature scalars vanish at the singularity, but there are causal geodesics along which the curvature components diverge [49] and, in this sense, the singularity is encountered for just some observers.
• Inaccessible singularities: These singularities appear in cosmological models with toral spatial sections, due to the infinite winding of trajectories around the tori—for instance, compactifying spatial sections of the de Sitter model to cubic tori. However, these singularities cannot be reached by physically well-defined observers, which prompts the name inaccessible singularities [50].

Figure 1 summarizes the classification of these singularities in a pedagogical way. All of the singularities discussed above have been studied in a variety of different contexts, and in this review, we would like to summarize works primarily of the past two decades on these topics and discuss the current status of such singularities. In Section 2, we give an overview of the subtleties of spacetime singularities, while in Section 3, we will discuss all the cosmological singularities mentioned above in detail. In Section 4, we discuss various methods that have been presented to remove such singularities (in some cases). In Section 5, we discuss a particular dynamical system analysis method (known as the Goriely–Hyde method) that has been shown to be very useful for cosmological singularity discussions. Finally, in Section 6, we summarize our brief review and discuss the future outlook for cosmology with regards to singularities.

2. An Overview of Spacetime Singularities

After Einstein proposed the general theory of relativity, which describes gravity in terms of spacetime curvature, the field equations were introduced to relate the geometry of spacetime to the matter content of the Universe. Early solutions included the Schwarzschild metric and the Friedmann models, which described the gravitational field around isolated objects and the overall geometry of the Universe, respectively. These models exhibited spacetime singularities where curvatures and energy densities became infinitely high, leading to a breakdown of the physical description. The Schwarzschild singularity at the center of symmetry could be eliminated by a coordinate transformation, but the genuine curvature singularity at $r=0$ remained. It was initially believed that these singularities were a result of the high symmetry in the models.

However, further research by Hawking, Penrose, Geroch, and others demonstrated that spacetime could have singularities under more general conditions. Singularities are an inherent feature of the theory of relativity and also apply to other gravitational theories based on spacetime manifolds. These singularities indicate super-ultra-dense regions in the Universe where physical quantities become infinitely large. In classical theories of gravity, singularities are an unavoidable aspect of describing physical reality. The behavior of these regions is beyond the scope of classical theory, and a quantum theory of gravity is needed to understand them. The field of gravitational physics saw significant developments in the 1960s due to observations of high-energy astrophysical phenomena and advancements in the study of spacetime structure and singularities. These advancements led to progress in black hole physics, relativistic astrophysics, and cosmology.

Singular behavior is observed in spacetime models described by general relativity. Examples include the Friedmann–Robertson–Walker (FRW) cosmological models and Schwarzschild spacetime. These models exhibit singularities where energy density and curvatures become infinitely large, leading to a breakdown of the conventional description of spacetime. Schwarzschild spacetime displays an essential curvature singularity at $r=0$, where the Kretschmann scalar $\alpha =R^{ijkl}{R}_{ijkl}$ diverges along any non-spacelike trajectory approaching the singularity. Similarly, for FRW models with $\rho +3p>0$ at all times (where $\rho$ is total energy density and p is pressure), a singularity arises at $t=0$, representing the origin of the Universe. Along past-directed trajectories approaching this singularity, both $\rho$ and the curvature scalar $R=R_{ij}{R}^{ij}$ become infinite. In both cases, past-directed non-spacelike geodesics are incomplete, and these essential singularities cannot be eliminated through coordinate transformations.

These singularities represent profound anomalies in spacetime, where the usual laws of physics fail. Geodesic incompleteness implies that a timelike observer will cease to exist in the spacetime after a finite amount of proper time. While singular behavior can occur without extreme curvature, such cases are considered artificial. An example is Minkowski spacetime with a removed point, where timelike geodesics encounter the hole and become future-incomplete. However, it is desirable to exclude such situations by requiring the spacetime to be “inextendible”, meaning that it cannot be isometrically embedded into a larger spacetime as a proper subset.

Nevertheless, non-trivial examples of singular behavior exist, such as conical singularities. These singularities do not involve diverging curvature components but are characterized by a Weyl-type solution. An example is the metric given by $d{s}^2=-d{t}^2+d{r}^2+r^2(d{\theta }^2+{sin}^2\theta d{\varphi }^2)$ with the identification $\varphi =0$ and $\varphi =a$ (with $a\ne 2\pi$), creating a conical singularity at $r=0$. The fundamental question is whether such singularities persist in general models and under what conditions they arise. Precisely defining a singularity in a general spacetime reveals that singularities likely exist in a broad range of spacetimes, subject to reasonable conditions. These singularities can emerge as the endpoint of gravitational collapse or in cosmological scenarios, such as the origin of the Universe.

The initial observation to make here is that, by its very definition, the metric tensor must possess a well-established meaning at every typical point within the spacetime. However, this principle ceases to hold at a spacetime singularity, like those previously discussed. Such a singularity cannot be considered a standard point within the spacetime; instead, it is a boundary point connected to the manifold. Consequently, difficulties arise when attempting to characterize a singularity based on the requirement that curvatures become infinite in proximity to it. The issue stems from the fact that, since the singularity lies outside the spacetime domain, it is not feasible to define its vicinity in the usual sense, which is essential for discussing the behavior of curvature quantities in that specific region.

An alternative approach might involve defining a singularity in relation to the divergence of elements within the Riemann curvature tensor along trajectories that do not follow spacelike directions. However, a challenge arises here as well: the behavior of these elements can change depending on the reference frames employed, rendering this approach less useful. One might consider utilizing curvature scalars or scalar polynomials involving the metric and Riemann tensor, demanding that they reach exceedingly large values. Instances of such divergence are encountered in models such as those of Schwarzschild and Friedmann. However, it remains possible that such a divergence only occurs at infinity for a given non-spacelike path. In a broader sense, it seems reasonable to expect some form of curvature divergence to occur along non-spacelike trajectories that intersect a singularity. Nevertheless, attempting to universally characterize singularities through curvature divergence encounters various complications.

Taking into account these scenarios and analogous ones, the presence of non-spacelike geodesic incompleteness is widely accepted as a criterion indicating the existence of a singularity within a spacetime. Although this criterion may not encompass all potential forms of singular behavior, it is evident that the occurrence of incomplete non-spacelike geodesics within a spacetime manifold signifies definite singular behavior. This manifests when a timelike observer or a photon abruptly vanishes from the spacetime after a finite interval of proper time or a finite value of the affine parameter. The singularity theorems, which emerge from an analysis of gravitational focusing and the global attributes of a spacetime, establish this incomplete nature for a broad array of spacetimes under a set of relatively general conditions.

From a physical standpoint, a singularity in any physics theory typically indicates that the theory becomes invalid either in the vicinity of the singularity or directly at the singularity. This implies a need for a broader and more comprehensive theory, necessitating a revision of the existing framework. Similar reasoning applies to spacetime singularities, suggesting that a description involving quantum gravity is warranted within these regions of the Universe, rather than relying solely on a classical framework.

The existence of an incomplete non-spacelike geodesic or an inextendible non-spacelike curve with a finite length, as measured by a generalized affine parameter, implies the presence of a spacetime singularity. The concept of “generalized affine length” for such a curve is defined as

$$L\left(\lambda \right)={\int }_0^a{\left[\sum _{i=0}^3{\left(X^i\right)}^2\right]}^{1/2}ds,$$

which remains finite. The component $X^i$ represents the tangent to the curve in a tetrad frame propagated in parallel along the curve. Each incomplete curve defines a boundary point of the spacetime, which is singular. To be considered a genuine physical singularity, it is expected that such a singularity is associated with the unbounded growth of spacetime curvatures. If all curvature components and scalar polynomials involving the metric and Riemann curvature tensor remain finite and well behaved as the singularity is approached along an incomplete non-spacelike curve, the singularity might be removable by extending the spacetime with relaxed differentiability requirements [51].

Different formalizations are possible for this requirement. A “parallely propagated curvature singularity” is one where the components of the Riemann curvature tensor are unbounded in a parallely propagated frame, forming the endpoint of at least one non-spacelike curve. Conversely, a “scalar polynomial singularity” occurs when a scalar polynomial involving the metric and Riemann tensor takes on infinitely large values along a non-spacelike curve ending at the singularity. This includes cases like the Schwarzschild singularity, where the Kretschmann scalar $(R^{ijkl}{R}_{ijkl}$ becomes infinite as r approaches 0. Curvature singularities, as further elucidated, also arise in various spacetime scenarios involving gravitational collapse. The strength of singularities and their potential to cause tidal forces on extended bodies can be assessed, and various criteria are available to determine this aspect [40]. These criteria all involve representing a finite object at each point along a causal geodesic as a volume defined by three independent Jacobi fields in the hyperspace, with the velocity of the curve as the normal vector. Tipler’s criterion [52] deems a singularity as strong if this volume tends to zero as the singularity is approached along the geodesic. On the other hand, Krolak’s criterion [53] stipulates that the derivative of this volume with respect to the normal parameter must be negative. Consequently, some singularities can be strong according to Krolak’s criterion while being weak according to Tipler’s, such as type III or Big Freeze singularities. Another criterion is outlined in [54].

Working with Jacobi fields involves solving the Jacobi equation along geodesics, a demanding task. However, conditions for lightlike and timelike geodesics, satisfying specific criteria, have been established [51]. These conditions involve integrals of the Ricci and Riemann curvatures of the spacetime metric along these curves.

• Lightlike geodesics:

  –

  According to Tipler’s criterion, a singularity is strong along a lightlike geodesic if and only if the integral

  $${\int }_0^{\tau }d{\tau }^{\prime }{\int }_0^{{\tau }^{\prime }}d{\tau }^{″}{R}_{ij}{u}^i{u}^j$$

  diverges as the proper time parameter $\tau$ approaches the singularity. Here, $u^i$ denotes the components of the velocity vector along the geodesic.

  –

  Krolak’s criterion states that the singularity is strong if and only if the integral

  $${\int }_0^{\tau }d{\tau }^{\prime }{R}_{ij}{u}^i{u}^j$$

  diverges as $\tau$ approaches the singularity.

• Timelike geodesics:

  –

  For timelike geodesics, Ref. [51] presents various necessary and sufficient conditions, though not a single characterization.

  –

  According to Tipler’s criterion for timelike geodesics, a singularity is strong if the integral

  $${\int }_0^{\tau }d{\tau }^{\prime }{\int }_0^{{\tau }^{\prime }}d{\tau }^{″}{R}_{ij}{u}^i{u}^j$$

  diverges as the proper time parameter approaches the singularity.

  –

  Krolak’s criterion for timelike geodesics specifies that the singularity is strong if the integral

  $${\int }_0^{\tau }d{\tau }^{\prime }{R}_{ij}{u}^i{u}^j$$

  diverges on approaching the singularity.

In passing, it is also of interest to talk of the cosmic censorship conjecture [55], which is the idea that all singularities arising from gravitational collapse will always be hidden by an event horizon. There are actually two versions of this conjecture: the weak version is that dynamical singularities in general relativity are generically not visible to observers at infinity, while the strong version is that dynamical singularities in general relativity are generically not visible to any observer. Singularities in violation of the weak version are dubbed globally naked, while those in violation of the strong version are dubbed locally naked. The conjectures have not yet been proven and have been a topic of recurring debates, motivating a lot of work on the topic of naked singularities. Several examples of spacetimes containing naked singularities have been found in recent times [56,57,58,59,60,61,62,63,64,65,66,67,68]. When such singularities develop in gravitational collapse, they give rise again to extremely intriguing physical possibilities and problems. The opportunity offered in this case is that we may observe the possible ultra-high energy physical processes occurring in such a region of the Universe, including the quantum gravity effects. In fact, loop quantum gravity in particular has a very amicable view of naked singularities and has been shown to be in favor of their existence [69,70]. Such observations of ultra-high-energy events in the Universe could provide observational tests and guide our efforts for a possible quantum theory of gravity, and so naked singularities could also be a good avenue for testing out the predictions of such theories. Recently, another very interesting study related to naked singularities was performed in [71]. The authors performed general relativistic ray tracing and radiative transfer simulations to generate synchrotron emission images utilizing the thermal distribution function for emissivity and absorptivity. They investigated the effects in the images of JMN-1 naked singularity (the Joshi–Malafarina–Narayan singularity, which has become quite a famous case of a naked singularity [65,68]) and a Schwarzschild black hole by varying the inclination angle, disk width, and frequency. Their results provided further support for naked singularities being a realistic scenario.

3. Types of Singularities

3.1. Strong Singularities

As mentioned before, strong singularities are those singularities that can distort finite objects in spacetime, and now we would like to discuss the prominent singularities of this category.

3.1.1. Big Bang Singularity (Type 0)

Classical models of the Universe generically feature an initial or "Big Bang" singularity. When we consider progressively earlier and earlier stages of the Universe, observable quantities stop behaving in a physically reasonable way. A more precise mathematical characterization of the cosmic Big Bang singularity can be achieved in terms of both a global notion of the incompleteness of inextendible causal (i.e., non-spacelike) past-directed curves and a local notion of the existence of a curvature pathology. Models of inflation also feature massive moving particles observing a singularity in a finite proper time. Hence, the Big Bang, as is widely known, is considered the singularity at the very beginning of the Universe. It is important to note a very interesting study on the completeness of inflationary spacetimes in the past direction that was undertaken by Borde, Guth, and Vilenkin [72], which highlighted the issue of incompleteness in certain inflationary scenarios. The applicability of their theorem varies across different cosmological frameworks. Notably, the emergent Universe paradigm, as elucidated by Ellis and Maartens [73], presents an intriguing alternative by proposing an inflationary cosmology without a singularity. The work of Guendelman et al. [74] delved into emergent cosmology, inflation, and dark energy, offering further insights into scenarios that evade the singularities discussed by Borde, Guth, and Vilenkin.).

3.1.2. Big Rip Singularity (Type 1)

In the case of a Big Rip singularity, the scale factor of the Universe becomes infinite at a finite time in the future, and the energy density and pressure density of the Universe also become infinite. In a Big Rip singularity, the dark energy density becomes so large that it causes the expansion of the Universe to accelerate at an ever-increasing rate. As a result, the scale factor of the Universe increases without bound, and the Universe becomes infinitely large at the time of the Big Rip. The energy density and pressure density of the Universe also become infinite at the time of the Big Rip. The thing to note here is that, interestingly, the Big Rip was proposed as a possible phantom scenario for the Universe [42], which means that the equation of state

$$w=\frac{p}{\rho }<-1$$

The phantom conclusion is interesting from the point of view that this presents some peculiar properties, like the energy density of phantom energy increasing with time or the fact that a phantom scenario violates the dominant energy condition [75]. Despite the fact that sound waves in quintessence travel at the speed of light, it should not be automatically assumed that disturbances in phantom energy must propagate faster than the speed of light. Indeed, there exist several scalar-field models for phantom energy where the sound speed is actually subluminal [76,77,78,79]. Phantom constructions have also been discussed in the context of quantum gravitational theories, for example in various string-theoretic realizations of dark energy [14,80,81]. Thus, it seems in principle interesting to look for a late-time Universe scenario with phantom dominance, and this is where the Big Rip comes in.

It is also worth discussing the subtleties of the Big Rip and how it would unfold. In a Universe resembling one with a cosmological constant, the scale factor’s expansion is faster than the Hubble distance, leading galaxies to gradually vanish beyond our observable horizon. If we introduce the concept of phantom energy, the expansion rate (Hubble constant) increases over time, causing the Hubble distance to shrink. Consequently, galaxies disappear at an accelerated pace as the cosmic horizon approaches. What is even more intriguing is the potential of the enhanced dark energy density to eventually tear apart objects held together by gravity. In the framework of general relativity, the gravitational potential’s source stems from the volume integral of the sum of energy density ($\rho$) and three times the pressure (p), denoted as $\rho +3p$.

For instance, a planet in orbit around a star with mass M and radius R becomes unbound approximately when the condition $-(4\pi /3)(\rho +3p)R^3\approx M$ is satisfied. In cases where the equation $-(\rho +3p)$ decreases over time due to a parameter w greater than or equal to −1, if $-(4\pi /3)(\rho +3p)R^3$ is smaller than M at present, it will continue to remain smaller indefinitely. This implies that any currently gravitationally bound system, such as the solar system, the Milky Way, the Local Group, and galaxy clusters, will remain bound in the future.

However, when dealing with phantom energy, the quantity $-(\rho +3p)$ increases with time. Consequently, at a certain point in time, every gravitationally bound system will eventually disintegrate. Analyzing the time evolution of the scale factor and the dependence of phantom-energy density on time, we deduce that a gravitationally bound system with mass M and radius R will undergo disintegration around a time $t\approx P\sqrt{2|1+3w|}/\left[6\pi \right|1+w\left|\right]$. Here, P represents the period of a circular orbit at radius R within the system. This process occurs prior to the Big Rip, with the earliest estimate of the Big Rip’s occurrence being 35 billion years. The Big Rip rips apart molecules and atoms, and even nuclei are dissociated (which makes it fitting that the name of the singularity is the Big Rip). However, it is not all gloomy, as various works have also explored ways to avoid the Big Rip (for example, [82]), and we will discuss these later on. Many more works have been published on various aspects of Big Rip singularities over the years—see [39,83,84,85,86,87,88,89,90,91,92,93,94]. A comparison of the Big Rip from dark energy and modified gravity was carried out for the first time in [95]. Less drastic variants of the Big Rip have also been found in recent years [96,97,98], which are discussed in detail in Appendix B.

3.1.3. Grand Bang and Grand Rip Singularities (Type $-1$)

The Grand Bang, although apparently different in name, is almost the same as the Big Bang singularity with null scale factors and a diverging pressure and energy density, but the one difference is that the singularity occurs with the equation of state parameter being equal to −1 [41]. This type of singularity was found initially by using a series ansatz for the scale factor. The Grand Bang and Grand Rip singularities are quite intricately linked with each other and so we shall discuss them now.

To discuss these grand singularities, we note that the equation for the parameter w is given by

$$w=\frac{p}{\rho }=-\frac{1}{3}-\frac{2}{3}\frac{a\ddot{a}}{{\dot{a}}^2}.$$

This expression holds true specifically for flat models. When considering curvature, additional terms need to be included.

The equation of state (EOS) parameter w has a close connection with the deceleration parameter q:

$$q=-\frac{a\ddot{a}}{{\dot{a}}^2}=\frac{1+3w}{2},$$

assuming flat models. Otherwise, the relationship between these parameters becomes more intricate, involving the Hubble parameter $H=\dot{a}/a$. This enables a direct translation of results from the EOS parameter to the deceleration parameter.

Alternatively, one can view this equation as the differential equation governing the evolution of the scale factor for a given time-dependent barotropic index $w\left(t\right)$. It is advantageous to introduce the variable $x=lna$:

$$\frac{\ddot{x}}{{\dot{x}}^2}=-\frac{3}{2}(w+1)=-(q+1)$$

This allows us to define

$$h\left(t\right):=\frac{3}{2}(w\left(t\right)+1)=q\left(t\right)+1$$

as a correction around the case of a pure cosmological constant:

$$w\left(t\right)=-1+\frac{2}{3}h\left(t\right),q\left(t\right)=-1+h\left(t\right).$$

This change of variables assists in reducing the order of the differential equation:

$$h=-\frac{\ddot{x}}{{\dot{x}}^2}={\left(\frac{1}{\dot{x}}\right)}^{·}\Rightarrow \dot{x}={\left(\int hdt+K_1\right)}^{-1}$$

In these terms, one finds the scale factor to be

$$a\left(t\right)=exp\left(\int \frac{dt}{\int h\left(t\right)dt}\right)$$

(1)

If one then assumes a power series form of $h\left(t\right)$ (which has become a quite well supported ansatz for the scale factor in various cosmological studies),

$$h\left(t\right)=h_0{t}^{{\eta }_0}+h_1{t}^{{\eta }_1}+\dots ,{\eta }_0<{\eta }_1<\dots ,$$

(2)

the energy and pressure densities can then be written as

$$\rho \left(t\right)=\left\{\begin{array}{cc}{\displaystyle 3{\left(\frac{{\eta }_0+1}{h_0}\right)}^2{t}^{-2({\eta }_0+1)}+\dots }\hfill & \mathrm{if}-1\ne {\eta }_0\ne 0\hfill \\ \\ {\displaystyle \frac{3}{h_0^2}\frac{1}{{ln}^2\left|t\right|}+\dots }\hfill & \mathrm{if}{\eta }_0=-1\hfill \\ \\ {\displaystyle \frac{3t^{-2}}{h_0^2}+\dots }\hfill & \mathrm{if}{\eta }_0=0,\hfill \end{array}\right.$$

and the pressure as

$$p\left(t\right)=\left\{\begin{array}{cc}{\displaystyle \frac{2{({\eta }_0+1)}^2}{h_0}{t}^{-{\eta }_0-2}+\dots }\hfill & \mathrm{if}-1\ne {\eta }_0<0\hfill \\ \\ {\displaystyle \frac{2}{h_0}\frac{1}{t{ln}^2\left|t\right|}+\dots }\hfill & \mathrm{if}{\eta }_0=-1\hfill \\ \\ {\displaystyle \frac{2h_0-3}{h_0^2}{t}^{-2}+\dots }\hfill & \mathrm{if}{\eta }_0=0\hfill \\ \\ {\displaystyle -3{\left(\frac{{\eta }_0+1}{h_0}\right)}^2{t}^{-2({\eta }_0+1)}+\dots }\hfill & \mathrm{if}{\eta }_0>0\hfill \end{array}\right.$$

This presents us with intriguing possibilities, but our focus will be on the case where ${\eta }_0>0$. In this scenario, we observe that at $t=0$, $\rho$ and p exhibit divergences following $t^{-2({\eta }_0+1)}$, and the parameter w converges to the value of $-1$. The consideration of such a singularity has not been explored within previous frameworks. The reason behind this omission is rooted in its incompatibility with the classifications established in [99,100]. This is due to the behavior of the scale factor (which is an exponential of rational functions); it does not lend itself to convergent power expansions, whether generalized or not, with a finite number of terms featuring negative powers. However, the function $x\left(t\right)$ does exhibit such behavior. The nature of the singularity is governed by the sign of the coefficient $h_0$. This is evident in the approximation of $a\left(t\right)$ as

$$a\left(t\right)\approx e^{-\mathrm{sgn}\left(h_0\right)\alpha /t^{{\eta }_0}},\alpha =\frac{{\eta }_0+1}{{\eta }_0\left|h_0\right|}>0,t>0$$

Based on this, we make the following observations:

• For $h_0>0$, the exponential term in Equation (1) decreases as t increases, and the scale factor a approaches zero as t approaches 0. This resembles an exponential-type Big Bang singularity or, if we swap t for $-t$, a Big Crunch. Given that $h_0$ is positive, the barotropic index w consistently remains below the phantom divide near $t=0$. Specifically, the value $w=-1$ is approached from values below it. These types of singularities are known as Grand Bang singularities.
• For $h_0<0$, conversely, the exponential term increases as t increases, causing the scale factor a to diverge to infinity as t approaches 0. This resembles an exponential-type Big Rip singularity at $t=0$, which, when considering the future, can be located by substituting t with $-t$. In this instance, the barotropic index w consistently remains above the phantom divide, and the value $w=-1$ is approached from values above it. This scenario is termed the Grand Rip singularity.

3.1.4. Directional Singularities (Type ∞)

The FLRW cosmological models are described by the metric

$$d{s}^2=-d{t}^2+a^2\left(t\right)\left\{f^2\left(r\right)d{r}^2+r^2\left(d{\theta }^2+{sin}^2\theta d{\varphi }^2\right)\right\},$$

(3)

where $f\left(r\right)=\frac{1}{\sqrt{1-k{r}^2}}$, and $k=0,\pm 1$. The spherical coordinates are $r,\theta$, and $\varphi$, and the coordinate time has a range depending on the cosmological model.

Free-falling observers follow timelike geodesics parametrized by proper time $\tau$. The velocity vector u is written as a unitary vector:

$$\delta =-g_{ij}{\dot{x}}^i{\dot{x}}^j,$$

(4)

where ${\dot{x}}^i$ denotes the derivative with respect to $\tau$. Three types of geodesics exist: timelike ($\delta =1$), spacelike ($\delta =-1$), and lightlike ($\delta =0$). We focus on causal geodesics, $\delta =0,1$.

For simplicity, we consider homogeneous and isotropic spacetimes, allowing us to set $\dot{\theta }=0=\dot{\varphi }$. The vector ${\partial }_R={\partial }_r/f\left(r\right)$, with R defined based on k, generates an isometry along straight lines. This leads to a conserved quantity P, the specific linear momentum of the observer:

$$\pm P=u·{\partial }_R=a^2\left(t\right)f\left(r\right)\dot{r},$$

where · represents the inner product defined by the metric (3), and the double sign ensures that P is positive.

To complete the set of equations, we use the unitarity condition (4):

$$\delta ={\dot{t}}^2-a^2\left(t\right)f^2\left(r\right){\dot{r}}^2.$$

It is evident that the equations governing the trajectories of geodesics, followed by observers not subject to acceleration ($\delta =1$) and lightlike particles ($\delta =0$) possessing a specific linear momentum P, can be simplified to

$$\frac{dt}{d\tau }=\sqrt{\delta +\frac{P^2}{a^2\left(t\right)}}$$

(5)

$$\frac{dr}{d\tau }=\pm \frac{P}{a^2\left(t\right)}$$

(6)

assuming a constant $\theta$ and $\varphi$ due to the symmetry inherent in these models. Here, $\tau$ represents the intrinsic or proper time as measured by the observer. In the context of null geodesics, we find that

$$\mathsf{\Delta }\tau =\frac{1}{P}{\int }_{-\infty }^{t}a\left(t\right)dt$$

(7)

Consequently, to ensure that the initial event $t=-\infty$ corresponds to a finite proper time interval $\mathsf{\Delta }\tau$ from an event at t, the requirement is

$${\int }_{-\infty }^{t}a\left(t\right)dt<\infty .$$

(8)

Therefore, the emergence of singular behavior exclusively at $t=-\infty$ is possible if the scale factor can be expressed as an integrable function of coordinate time. This condition necessitates that $a\left(t\right)$ tends towards zero as t approaches $-\infty$, although this alone is not sufficient. Similarly, for timelike geodesics with a non-zero P,

$$\mathsf{\Delta }\tau ={\int }_{-\infty }^t\frac{dt}{\sqrt{1+\frac{P^2}{a^2\left(t\right)}}}<\frac{1}{P}{\int }_{-\infty }^{t}a\left(t\right)dt,$$

(9)

indicating that the proper time interval to $t=-\infty$ is finite provided that the time interval for lightlike geodesics is also finite. Consequently, $t=-\infty$ is reachable for these observers. As a result, condition (8) implies that both lightlike and timelike geodesics with a non-zero P experience $t=-\infty$ within a finite proper time interval in their past. Conversely, comoving observers tracing timelike geodesics with $P=0$ exhibit $d\tau =dt$, which leads to $t=-\infty$ corresponding to an infinite proper time interval in their past; thus, they cannot encounter the singularity. This dichotomy is responsible for the directional nature of Type ∞ singularities, as they are accessible to causal geodesics, except those with $P=0$. Ultimately, it can be concluded that Type ∞ singularities can manifest in three scenarios:

• For a finite ${\int }_{-\infty }hdt$ with $h\left(t\right)>0$: $a_{-\infty }=0$, ${\rho }_{-\infty }=\infty$, $p_{-\infty }=-\infty$, $w_{-\infty }=-1$. These differ from the “Little Rip” model in the sign of $h\left(t\right)$, and are termed a “Little Bang” if they denote an initial singularity, or a “Little Crunch” if they represent a final singularity [101]. Instances of this case encompass models with a scale factor $a\left(t\right)\propto e^{-\alpha {(-t)}^p}$ where $p>1$ and $\alpha >0$.
• When $h_{-\infty }=0$ and $\left|h\left(t\right)\right|\gtrsim {\left|t\right|}^{-1}$ with $h\left(t\right)<0$: $a_{-\infty }=0$, ${\rho }_{-\infty }=0$, $p_{-\infty }=0$, $w_{-\infty }=-1$. Changing the sign of $h\left(t\right)$ gives rise to a variant of the “Little Rip” scenario, featuring an asymptotically vanishing energy density and pressure. Models with a scale factor $a\left(t\right)\propto e^{-\alpha {(-t)}^p}$ where $p\in (0,1)$ and $\alpha >0$ exemplify this case.
• For a finite $h_{-\infty }\in (-1,0)$: $a_{-\infty }=0$, ${\rho }_{-\infty }=0$, $p_{-\infty }=0$, and a finite $w_{-\infty }\ne -1$. This case applies to models like $a\left(t\right)\propto t^{-p}$ with $p>1$, as explored in [49].

While they have recently been discussed in the context of inflationary models [101], not much work has been carried out on Type ∞ singularities since their discovery with regard to their avoidance or their emergence in more exotic cosmological models.

3.2. Weak Singularities

3.2.1. Sudden Singularities (Type II)

In the case of such Type II singularities, the pressure density diverges, or, equivalently, the derivatives of the scale factor diverge from the second derivative onwards. Let us start by examining informally whether there is potential for the emergence of singularities in which a physical scalar quantity becomes unbounded at a finite future comoving proper time $t_s$. This might occur when the scale factor $a\left(t\right)$ approaches a non-zero or infinite value $a\left(t_s\right)$ and the Hubble parameter $H\left(t\right)$ approaches a finite value $H_s$ (where $H_s$ is positive and not infinite). If such scenarios are feasible, the following conditions need to be satisfied:

$$\rho \to 3H_s^2+\frac{k}{a_s^2}={\rho }_s<\infty$$

(10)

and

$$\frac{\ddot{a}}{a}\to \frac{p}{2}-\frac{H_s^2}{2}-\frac{k}{6a_s^2}$$

(11)

$$\dot{\rho }\to -3H_s({\rho }_s+p)$$

(12)

Hence, it becomes apparent that the density must inevitably remain finite at $t_s$. However, there is still a possibility for a singularity in pressure to arise, manifesting as

$$p\left(t\right)\to \infty$$

(13)

as $t\to t_s$, consistent with the conditions outlined in Equation (11). In such instances, the pressure singularity is concomitant with an infinite acceleration.

To illustrate this, we take the most primitive example of such a singularity, which was put forward by Barrow in [43]. In this regard, assume that it is physically reasonable to expect that the scale factor can be written in the form of the following ansatz (Appendix A provides a detailed overview of the motivations that allow such a consideration for the scale factor):

$$a\left(t\right)=A+B{t}^q+C{(t_s-t)}^n,$$

(14)

where $A>0$, $B>0$, $q>0$, C, and $n>0$ are constants that we will determine. We set the origin of time such that $a\left(0\right)=0$, leading to $A=-C{t}_s^n>0$. Consequently, we find the expression for the Hubble parameter $H_s$:

$$H_s=\frac{qB{t}_s^{q-1}}{A+B{t}_s}.$$

(15)

For simplicity, we use the freedom to scale the Friedmann metric by dividing by A and set $A\equiv 1$ and $C\equiv -t_s^n$. This yields the simplified form of $a\left(t\right)$:

$$a\left(t\right)={\left(\frac{t}{t_s}\right)}^q(a_s-1)+1-{\left(1-\frac{t}{t_s}\right)}^n,$$

(16)

where $a_s\equiv a\left(t_s\right)$. As t approaches $t_s$ from below, the behavior of the second derivative of a can be described:

$$\ddot{a}\to q(q-1)B{t}^{q-2}-\frac{n(n-1)}{t_s^2{(1-\frac{t}{t_s})}^{2-n}}\to -\infty ,$$

(17)

whenever $1<n<2$ and $0<q\le 1$. This solution is valid for $0<t<t_s$. Consequently, as t approaches $t_s$, a approaches $a_s$; $H_s$ and ${\rho }_s>0$ (as long as $3q^2{(a_s-1)}^2{t}_s^{-2}>-k$) remain finite; and $p_s\to \infty$.

When $2<n<3$, $\ddot{a}$ remains finite but $\stackrel{⃛}{a}\to \infty$ as t approaches $t_s$. Here, $p_s$ remains finite, but ${\dot{p}}_s\to \infty$. In contrast, there exists an initial strong-curvature singularity, where both $\rho$ and p tend to infinity as t approaches 0. Importantly, in this scenario, both $\rho$ and $\rho +3p$ remain positive. Such behavior can even arise in a closed universe ($k=+1$), where the pressure singularity prevents expansion from reaching a maximum. This is the most primitive example of a pressure singularity, but ever since the work in [43] was presented, such singularities have been discussed in many different settings, both from modified gravity perspectives and based on other phenomenological considerations. Work has also been carried out on ways to escape such singularities, which we will discuss later in this paper.

3.2.2. Big Freeze Singularity (Type III)

The Big Freeze singularity is similar to the Big Rip but is still quite different from it. This singularity was first shown in a phantom generalized Chaplygin gas (PGCG) cosmology in [46], and we shall quickly see how it unfolds in such a scenario.

The equation of state governing PGCG closely resembles that of the conventional generalized Chaplygin gas. It can be succinctly expressed as

$$p=-\frac{A}{{\rho }^{\alpha }},$$

where the symbol A represents a positive constant, and $\alpha$ signifies a parameter. In the scenario where $\alpha =1$, the equation assumes the form of a simple Chaplygin gas equation of state. This relationship is crucially connected to the continuity equation, given by

$$\begin{array}{c}\hfill \dot{\rho }+3H(p+\rho )=0\end{array}$$

(18)

from which emerges the expression for energy density $\rho$:

$$\rho ={\left(A+\frac{B}{a^{3(1+\alpha )}}\right)}^{\frac{1}{1+\alpha }},$$

where B stands as a constant parameter. In a noteworthy observation made in Ref. [87], it was discerned that a negative B renders the perfect fluid, with the equation of state $p=-\frac{A}{{\rho }^{\alpha }}$, unable to uphold the null energy condition, that is, $p+\rho <0$. Intriguingly, under these conditions, the energy density escalates as the Universe expands, contrary to redshift behavior, thus earning the label “phantom generalized Chaplygin gas” (PGCG).

Further insights from the works of [87,102] revealed that for a PGCG with $\alpha >-1$, an FLRW universe hosting this fluid can evade the impending Big Rip singularity. As the scale factors attain far greater magnitudes, the Universe eventually approximates an asymptotically de Sitter state. In stark contrast, during the Big Freeze scenario, the PGCG energy density responds by increasing as the scale factor matures. Specifically, as the scale factor approaches minuscule values ($a\to 0$), $\rho$ tends towards $A^{\frac{1}{1+\alpha }}$, while it experiences a surge at a finite scale factor $a_{\max }$:

$$a_{\max }={\left|\frac{B}{A}\right|}^{\frac{1}{3(1+\alpha )}}.$$

As a consequence, an FLRW universe saturated with PGCG is destined to confront a finite-radius future singularity. Notably, the vicinity of this singularity lends itself to a cosmological evolution described by the relation

$$a\simeq a_{\max }\left\{1-{\left[\frac{1+2\alpha }{2(1+\alpha )}\right]}^{\frac{2(1+\alpha )}{1+2\alpha }}{A}^{\frac{1}{1+2\alpha }}{\left|3(1+\alpha )\right|}^{\frac{1}{1+2\alpha }}{(t_{\max }-t)}^{\frac{2(1+\alpha )}{1+2\alpha }}\right\}.$$

Remarkably, this singularity emerges at not only a finite scale factor but also a distinct future cosmic time. Conversely, the history of an FLRW universe permeated with this fluid traces back to an asymptotically de Sitter state in the past. This temporal journey is expressed succinctly as follows:

$$a\simeq a_{0}exp\left(A^{\frac{1}{2(1+\alpha )}}t\right),$$

where $a_0$ signifies a minute scale factor. Additionally, the Universe embarks on its odyssey from a bygone infinity of cosmic time as $a\to 0$ and $p+\rho \to 0^{-}$. Remarkably, the homogeneous and isotropic nature of the Universe propels it into a phase of super-accelerated expansion, denoted by

$$\dot{H}=-\frac{3}{2}(p+\rho )>0,$$

until it culminates at the singularity $a=a_{\max }$. It is imperative to recall that the PGCG eludes the satisfaction of the null energy condition [87], as embodied in $p+\rho <0$. A great amount of work has been carried out on Big Freeze singularities since the initial study presented in [46], and it has been shown that one can encounter such singularities in many exotic cosmological settings. Furthermore, there have also been studies probing how one can avoid such singularities [103,104,105,106,107,108,109,110].

3.2.3. Generalized Sudden Singularities (Type IV)

These singularities were first discussed in [30] and have since been found in a diverse variety of cosmological settings. Thus, here we will briefly discuss the primary cases in which Type IV singularities have been shown. In fact, the following example will illustrate all the prominent singularities we have discussed so far. We start with an equation of state of the form

$$p=-\rho -f\left(\rho \right)$$

(19)

This sort of equation of state with $f\left(\rho \right)=A{\rho }^{\alpha }$, where $\alpha$ is an arbitrary constant, was first proposed in [47] and was investigated in detail in [83]. There can be diverse physical motivations behind such an equation of state. This form of an EOS can also be equivalent to bulk viscosity [111] or come about due to modified gravity effects [29]. We now consider the following ansatz for the scale factor:

$$a\left(t\right)=a_0{\left(\frac{t}{t_s-t}\right)}^n,$$

(20)

where n is a positive constant, and $0<t<t_s$. The scale factor diverges within a finite time ($t\to t_s$), resembling the phenomenon of the Big Rip singularity. Consequently, $t_s$ represents the Universe’s lifetime. When $t\ll t_s$, the evolution of $a\left(t\right)$ follows $t^n$, leading to an effective EOS given by $w=-1+2/\left(3n\right)>-1$. Conversely, when $t\sim t_s$, the effective EOS assumes $w=-1-2/\left(3n\right)<-1$. The Hubble rate in this case can be expressed as

$$H=n\left(\frac{1}{t}+\frac{1}{t_s-t}\right).$$

(21)

Utilizing Equation (21), one can deduce the relation

$$\rho =\frac{3n^2}{{\kappa }^2}{\left(\frac{1}{t}+\frac{1}{t_s-t}\right)}^2.$$

(22)

As a result, both H and $\rho$ exhibit minima at $t=t_s/2$, characterized by the values

$$H_{\min }=\frac{4n}{t_s},{\rho }_{\min }=\frac{48n^2}{{\kappa }^2{t}_s^2}.$$

(23)

Next, we examine a specific form for $f\left(\rho \right)$ given by

$$f\left(\rho \right)=\frac{AB{\rho }^{\alpha +\beta }}{A{\rho }^{\alpha }+B{\rho }^{\beta }},$$

(24)

where A, B, $\alpha$, and $\beta$ are constants. As we shall see, this dark energy scenario harbors a complex structure with respect to singularities.

In scenarios where $\alpha$ surpasses $\beta$, we observe that

$$f\left(\rho \right)\to \left\{\begin{array}{cc}A{\rho }^{\alpha }\hfill & \mathrm{as}\rho \to 0\hfill \\ B{\rho }^{\beta }\hfill & \mathrm{as}\rho \to \infty \hfill \end{array}\right..$$

(25)

For non-unit values of $\alpha$ and $\beta$, we obtain

$$a=a_{0}exp\left\{-\frac{1}{3}\left[\frac{{\rho }^{-\alpha +1}}{(\alpha -1)A}+\frac{{\rho }^{-\beta +1}}{(\beta -1)B}\right]\right\}.$$

(26)

The realm of possibilities in this cosmology is extensive. If $1>\alpha >\beta$ and $A,B>0$ ($A,B<0$), the scale factor has a minimum (maximum) at $\rho =0$, extending to infinity (vanishing) as $\rho \to \infty$. When $\alpha >1>\beta$ and $A<0$ while $B>0$ ($A>0$ and $B<0$), the scale factor features a minimum (maximum) at a non-trivial (non-vanishing) $\rho$ value, reaching infinity (zero) as $\rho$ approaches zero or a positive infinity. For $\alpha >1>\beta$ and $A,B>0$ ($A,B<0$), the scale factor becomes infinite (vanishes) as $\rho \to \infty$ ($\rho \to 0$), and it vanishes (increases) as $\rho \to 0$ ($\rho \to \infty$). When $\alpha >\beta >1$, the scale factor approaches $a_0$ as $\rho \to \infty$. Additionally, if $A>0$ ($A<0$), the scale factor tends to 0 (∞) as $\rho \to 0$. With $A,B>0$ ($A,B<0$), the scale factor demonstrates a monotonic increase (decrease) regarding $\rho$. In the case of $A>0$ and $B<0$ ($A<0$ and $B>0$), the scale factor attains a non-trivial maximum (minimum) at a finite $\rho$ value.

To summarize, the possibilities for singularity formation in this cosmological model are remarkably diverse. It is worth noting that some of the identified singularities may violate one or more energy conditions. These energy conditions encompass the following:

$$\begin{array}{cc}\hfill & \rho \ge 0\rho \pm p\ge 0“\mathrm{dominant}\mathrm{energy}\mathrm{condition}”\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & \rho +p\ge 0“\mathrm{null}\mathrm{energy}\mathrm{condition}”\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & \rho \ge 0\rho +p\ge 0“\mathrm{weak}\mathrm{energy}\mathrm{condition}”\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & \rho +3p\ge 0\rho +p\ge 0“\mathrm{strong}\mathrm{energy}\mathrm{condition}”\hfill \end{array}$$

(30)

With these considerations, we can succinctly summarize the findings for the cosmological model defined by the $f\left(\rho \right)$ function as follows:

• For $A/B<0$, a Type II singularity is inevitable, irrespective of the values of $\beta$.
• Regardless of the sign of $A/B$, the nature of singularities varies according to the values of $\beta$:
  • $0<\beta <1/2$: A Type IV future singularity is evident. The parameter w approaches infinity ($-\infty$) for $B<0$ ($B>0$).
  • $\beta >1$: A Type III future singularity emerges, accompanied by a breach of the dominant energy condition. The parameter w approaches infinity ($-\infty$) for $B<0$ ($B>0$).
  • $3/4<\beta <1$: A Type I future singularity emerges if $A>0$. The dominant energy condition is violated for $A>0$, and w approaches $-1+0$ ($-1-0$) for $A<0$ ($A>0$).
  • $1/2\le \beta \le 3/4$: No finite future singularity is present.
  • $\beta =0$: A finite future singularity is absent, yet as $\rho \to 0$, w approaches infinity ($-\infty$) for $B<0$ ($B>0$).
  • $\beta <0$: A Type II future singularity emerges. The dominant energy condition is broken, though the strong energy condition remains intact for $B<0$. The parameter w approaches infinity ($-\infty$) for $B<0$ ($B>0$).

Thus, this example (as was discussed in [30]) shows us how one can find not only Type IV singularities but also the other singularities we have discussed so far. Another interesting thing to note is that there turn out to be qualitative differences when one considers singularities in Jordan and Einstein frames, something which was discussed in detail and discovered in [112,113]. It is also worth noting that when one considers viscous fluids, as in [114], then different types of singularities may arise. The occurrence of singularities in an oscillating universe has also been discussed, first in [115]. Singularities have also been considered in detail for bounce cosmologies [116,117,118]. The realization of all four known types of future singularities (Type I-Type IV) has also been found in very exotic modified gravity theories, for example, in an f(R) version of Horava–Lifschitz gravity [119], while also in teleparallel constructions like the one considered in [120].

A crucial point that we should note here in passing with regards to all the singularities we have discussed so far is that the tidal forces manifest for these singularities as the (infinite) impulse that reverses (or stops) the increase in the separation of geodesics, and the geodesics themselves can evolve further; the Universe can then continue its evolution through a singularity. Moreover, it is intriguing to consider the potential consequences of these singularities on the constructs of quantum gravity. Although there exists a considerable body of literature exploring the emergence of cosmological singularities in quantum gravitational scenarios like braneworlds, for instance, a more profound inquiry pertains to the influence of such singularities on fundamental entities like strings.

If we contemplate an elongated structure such as a classical string, modeled using the Polyakov formalism [121]

$$S=-\frac{T}{2}\int d\tau d\sigma {\eta }^{\mu \nu }{g}_{ab}{\partial }_{\mu }{X}^a{\partial }_{\nu }{X}^b$$

(31)

(with T denoting string tension; $\tau ,\sigma$ representing the string’s worldsheet coordinates; ${\eta }^{\mu \nu }$ corresponding to the worldsheet metric; $\mu ,\nu =0,1$; and $g_{ab}$ standing for the spacetime metric), the scenario involves the string interacting with a non-BB singularity [122]. The crux of the matter is that a measurable property of the string, its invariant size $S\left(\tau \right)=2\pi a\left(\eta \right(\tau \left)\right)R\left(\tau \right)$ (using a circular assumption with radius R), reveals certain characteristics. Specifically, at a Big Rip singularity, the string undergoes infinite stretching ($S\to \infty$), resulting in its destruction. In contrast, at a Type II singularity, the scale factor remains finite at the $\eta$-time, consequently maintaining a finite invariant string size. Analogously, the same holds true for Type III and Type IV singularities. This implies that strings remain intact when encountering such singularities. This also underscores the “weakness” of these singularities in the sense that they do not display geodesic incompleteness. As a result, particles [123] and even more extensive entities like extended objects [122] can traverse them without obstruction. Hence, they lack a “dangerous” quality, which explains their potential emergence in the relatively proximate future (for instance, in around 10 million years for Type II, or the idea that a pressure singularity has happened in the recent past) [38,124,125,126].

3.2.4. w-Singularities (Type V)

As the name suggests, w-singularities occur when the equation of state parameter (w) blows up in some cosmological models. The singularities were first introduced in [127] and then expanded upon in later works [48,128]. The authors in [127] arrived at w-singularities by first choosing the scale factor ansatz as follows:

$$a\left(t\right)=A+B{\left(\frac{t}{t_s}\right)}^{\frac{2}{3\gamma }}+C{\left(D-\frac{t}{t_s}\right)}^n.$$

(32)

This contains seven arbitrary constants: A, B, C, D, $\gamma$, n, and $t_s$. The last of the constants $t_s$ is the time when we expect the singularity. Using the scale factor (32), the authors imposed the following conditions:

$$a\left(0\right)=0,a\left(t_s\right)=const.\equiv a_s,\dot{a}\left(t_s\right)=0,\ddot{a}\left(t_s\right)=0.$$

(33)

The first of the conditions (33) was chosen in order for the evolution to begin with a standard Big Bang singularity at $t=0$ (note that in order to have a Big Rip, one would have to impose $a\left(0\right)=\infty$, which is equivalent to taking $\gamma <0$). One can see that after introducing (33), the energy density and the pressure vanish at $t=t_s$. The model does not admit a singularity of the higher derivatives of the Hubble parameter since $\ddot{H}\left(t_s\right)\ne 0$ in $\ddot{H}$, and so it is not a Type IV singularity according to the classification of Ref. [30]. On the other hand, even though both $\ddot{a}\left(t_s\right)$ and $\dot{a}\left(t_s\right)$ vanish in the limit $t\to t_s$, the deceleration parameter blows up to infinity, i.e.,

$$q\left(t_s\right)=-\frac{\ddot{a}\left(t_s\right)a_s}{{\dot{a}}^2\left(t_s\right)}\to \infty$$

(34)

Consequently, one can find that the EOS parameter is related to the deceleration parameter as follows:

$$w\left(t\right)=\frac{c^2}{3}\left[2q\left(t\right)-1\right].$$

(35)

Thus, one finds that $w\left(t_s\right)\to \infty$. Then, we face a very strange singularity. It has vanishing pressure and energy density and a constant scale factor, but the deceleration parameter and, in particular, the time-dependent barotropic index $w\left(t\right)$ are singular. Another ansatz for the scale factor that can give w-singularities was proposed by Dabrowski and Marosek in [129] and has an exponential form. This ansatz is given by

$$a\left(t\right)=a_s{\left(\frac{t}{t_s}\right)}^{m}exp{\left(1-\frac{t}{t_s}\right)}^n$$

(36)

where $a_s$ has the units of length and is a constant, and m and n are also constants. While the ansatz on the surface looks quite different from a power series ansatz, which we will consider later on, it can be a subcase of a series ansatz within certain limits as well. The scale factor is zero (a = 0) at t = 0, thus signifying the Big Bang singularity. One can write the first and second derivatives of the scale factor as

$$\dot{a}\left(t\right)=a\left(t\right)\left[\frac{m}{t}-\frac{n}{t_s}{\left(1-\frac{t}{t_s}\right)}^{n-1}\right]$$

(37)

$$\ddot{a}\left(t\right)=\dot{a}\left(t\right)\left[\frac{m}{t}-\frac{n}{t_s}{\left(1-\frac{t}{t_s}\right)}^{n-1}\right]+a\left(t\right)\left[-\frac{m}{t^2}+\frac{n(n-1)}{t^2}{\left(1-\frac{t}{t_s}\right)}^{n-2}\right]$$

(38)

where the overdots denote differentiation with respect to time. From this, one can see that for $1<n<2$, $\dot{a}\left(0\right)\to \infty$ and $\dot{a}\left(t_s\right)=\frac{m{a}_s}{t_s}=$ const., while $a\left(t_s\right)=a_s$, $\ddot{a}\left(0\right)\to \infty$, and $\ddot{a}\left(t_s\right)\to -\infty$ and we have sudden future singularities. Furthermore, it was shown in [129] that for the simplified case of the scale factor (20) with $m=0$, one can obtain w-singularities for $n>0$ and $n\ne 1$. Finally, yet another ansatz to obtain w-singularities was provided in [48] and is of a power series form, given by

$$a\left(t\right)=c_0+c_1{(t_s-t)}^{n1}+c_2{(t_s-t)}^{n2}\dots$$

(39)

where $t_s$ is the time of the singularity. In order for pressure to be finite, $n_1>1$. There have of course been a significant number of works that have considered how these singularities can occur in non-standard cosmologies and how they can be avoided. However, in passing, a discussion of the cosmological significance of w-singularities is in order. While Type I–Type IV singularities deal with more direct cosmological parameters like the scale factor and Hubble parameter alongside energy and pressure densities, Type V singularities deal with a somewhat indirect parameter in the form of w. This is not to say, however, that these singularities cannot occur in cosmological and, in particular, dark energy models. For example, Ref. [130] discussed how w-singularities can occur in interacting dark energy models (the background cosmology in this case was still general relativistic, and the continuity equation had its usual form), while [131] showed how varying Chaplygin gas models can also have w-singularities. The occurrence of w-singularities in various other contexts has also been discussed in [91,132,133,134,135]. Hence, while Type V singularities deal primarily with a more indirect cosmological parameter, they by no means diminish their cosmological importance and they do appear in a variety of cosmological expansion scenarios.

4. Singularity Removal/Avoidance Methods

With the increasing interest in finding singularities in cosmological models, a natural interest also grew in investigating ways in which such singularities could either be completely removed or at least mildly alleviated/avoided in some cases. This has also resulted in an impressive amount of literature (for example, refer to [39] for a detailed account of avoiding singularities in both Jordan and Einstein frames). We would like to discuss some of the prominent works in this area, focusing on the use of quantum effects and modified gravity effects to deal with singularities.

4.1. Conformal Anomaly Effects near Singularities

The effect of the quantum backreaction of conformal matter around Type I, Type II, and Type III singularities was taken into consideration in the works of Nojiri and Odintsov [28,30,136]. In these cases, the curvature of the Universe becomes large around the singularity time $t=t_s$, although the scale factor a is finite for Type II and III singularities. Since quantum corrections usually contain the powers of the curvature or higher derivative terms, such correction terms are important near the singularity. At this point, it becomes important to add some context regarding what conformal anomalies are and how they are usually perceived in high-energy physics. It is fair to assume that there were many matter fields during inflation in the early Universe because the standard model of particle physics has almost 100 fields, and this number may increase by two if the standard model is contained in a supersymmetric theory. Although the behavior of these (massless) matter fields—scalars, the Dirac spinors, and vectors in curved spacetime—is conformally invariant, some divergences are observed because of the presence of the one-loop vacuum contributions. In the renormalized action, some counterterms are required to break the matter action’s conformal invariance in order to cancel the poles of the divergence component. From the classical point of view, the trace of the energy momentum tensor in a conformally invariant theory is null. However, renormalization procedures can lead to the trace of an anomalous energy momentum tensor, which is the so-called quantum anomaly or conformal anomaly (we would recommend the reader refer to [137,138,139,140] for more details on conformal anomaly effects). The conformal anomaly we described can be considered to have the following form [30]:

$$T_A=b\left(F+\frac{2}{3}\square R\right)+b^{\prime }G+b^{″}\square R$$

(40)

where $T_A$ is the trace of the stress energy tensor, F is the square of the 4D Weyl tensor, and G is a Gauss–Bonet curvature invariant, which are given by

$$F=(1/3)R^2-2R_{ij}{R}^{ij}+R_{ijkl}{R}^{ijkl}$$

(41)

$$G=R^2-4R_{ij}{R}^{ij}+R_{ijkl}{R}^{ijkl}$$

(42)

b and $b^{\prime }$, on the other hand, are given by

$$b=\frac{N+6N_{1/2}+12N_1+611N_2-8N_{HD}}{120{\left(4\pi \right)}^2}$$

(43)

$$b^{\prime }=-\frac{N+11N_{1/2}+62N_1+1411N_2-28N_{HD}}{360{\left(4\pi \right)}^2}$$

(44)

with N being a scalar, $N_{1/2}$ a spinor, $N_1$ vector fields, $N_2$ (= 0 or 1) gravitons, and $N_{HD}$ higher derivative conformal scalars. For usual matter, $b>0$ and $b^{\prime }<0$, except for higher derivative conformal scalars, while $b^{″}$ can be arbitrary. Quantum effects due to the conformal anomaly act as a fluid with energy density ${\rho }_A$ and pressure $p_A$. The total energy density is ${\rho }_{tot}=\rho +{\rho }_A$. The conformal anomaly, also known as the trace anomaly, can be given by the trace of the fluid stress energy tensor

$$T_A=-{\rho }_A+3p_A$$

(45)

The conformal-anomaly-corrected pressure and energy densities still obey the continuity equation (18).Using this, we can write

$$T_A=-4{\rho }_A-\frac{\dot{{\rho }_A}}{H}\left(\frac{1}{2{\rho }_A}-1\right)$$

(46)

The conformal-anomaly-corrected pressure and energy densities still obey the continuity equation. Using this, we can write [30]

$$T_A=-4{\rho }_A-\frac{\dot{{\rho }_A}}{H}$$

(47)

One can then express ${\rho }_A$ as an integral in terms of $T_A$ as

$${\rho }_A=-\frac{1}{a^4}\int a^{4}H{T}_{A}dt$$

(48)

Furthermore, $T_A$ can be expressed in terms of the Hubble parameter as

$$T_A=-12b{\dot{H}}^2+24b^{\prime }(-{\dot{H}}^2+H^2\dot{H}+H^4)-(4b+6b^{″})(H^{(3)}+7H\ddot{H}+4{\dot{H}}^2+12H^2\dot{H})$$

(49)

Using this, one can obtain an expression for ${\rho }_A$ taking into account conformal anomaly effects near the singularity:

$${\rho }_A=-\frac{1}{a^4}\int dt{a}^{4}H{T}_A=-\frac{1}{a^4}\int dt{a}^{4}H\left[-12b{\dot{H}}^2+24b^{\prime }(-{\dot{H}}^2+H^2\dot{H}+H^4)-(4b+6b^{″})\left(\stackrel{⃛}{H}+7H\ddot{H}+4{\dot{H}}^2+12H^2\dot{H}\right)\right]$$

(50)

The quantum-corrected Friedmann equation is

$$\frac{3}{{\kappa }^2}{H}^2=\rho +{\rho }_A.$$

(51)

Note that to maintain consistency with the notation used in [30], we consider the Friedmann equation to be of the form $H^2=\frac{{\kappa }^2}{3}(\rho +{\rho }_m).$ Since the curvature is expected to be large near the time of the singularity, one can assume that $(3/{\kappa }^2)H^2\ll \left|{\rho }_A\right|$. Then, $\rho \sim -{\rho }_A$ from (50), which gives

$$\dot{\rho }+4H\rho =H\left[-12b{\dot{H}}^2+24b^{\prime }(-{\dot{H}}^2+H^2\dot{H}+H^4)-(4b+6b^{″})\left(\stackrel{⃛}{H}+7H\ddot{H}+4{\dot{H}}^2+12H^2\dot{H}\right)\right]$$

(52)

Finally, the continuity equation $\dot{\rho }+3H\left(\rho +p\right)=0$ for $p=-\rho -f\left(\rho \right)$ gives

$$H=\frac{\dot{\rho }}{3f\left(\rho \right)}.$$

(53)

Now, we can appreciate the implications of these effects on both strong and weak singularities. First, we consider the Big Rip. The first attempt to address the issue of the Big Rip with conformal anomalies was presented in [141,142]. For this, we consider the model given by

$$f\left(\rho \right)\sim B{\rho }^{\beta }$$

(54)

with $1/2<\beta <1$ when $\rho$ is large. In this case, there exists the Big Rip singularity, as we discussed previously in Section 4. We note that the classical evolution is characterized by $\rho \propto {(t_s-t)}^{\frac{2}{1-2\beta }}$ and $H\propto {(t_s-t)}^{\frac{1}{1-2\beta }}$, both of which exhibit divergence for $\beta >1/2$. When quantum corrections are taken into account, it is natural to assume that near the singularity $\rho$ behaves as follows:

$$\rho ={\rho }_0{\left(t_s-t\right)}^{\tilde{\gamma }}$$

(55)

As $\rho$ may diverge at $t=t_s$, we consider negative values of $\tilde{\gamma }$. Since $\tilde{\gamma }\left(1-\beta \right)<0$ in this case, we might expect that (52) would give the following approximate relation around $t=t_s$:

$$\rho \sim 6b^{\prime }{H}^4$$

(56)

The term on the right hand side grows as $H^4\propto {(t_s-t)}^{-4+4\tilde{\gamma }(1-\beta )}$, but this does not give a consistent result, since $\rho$ becomes negative for $b^{\prime }<0$. This tells us that our assumptions are wrong, and $\rho$ does not become infinite. If $\rho$ has an extremum, (53) tells us that H vanishes there since $\dot{\rho }=0$. Furthermore, the authors of [30] showed numerically that in this scenario the Hubble rate approaches zero in finite time, thus coming to the conclusion that conformal anomaly effects can alleviate the Big Rip in this case.

Let us again consider the model in (54), but now for the range $\beta >1$, in which case we see that a Type III singularity develops with $\rho \propto {\left(t_s-t\right)}^{\frac{2}{1-2\beta }}$. Again, we consider that near the singularity $\rho$ behaves as in (55). Using (53), one finds that

$$H=-\frac{\tilde{\gamma }{\rho }_0^{1-\beta }}{3B}{\left(t_s-t\right)}^{-1+\tilde{\gamma }\left(1-\beta \right)}$$

(57)

Since we are considering the case $\beta >1$ and $\tilde{\gamma }<0$, we know that $\tilde{\gamma }\left(1-\beta \right)>0$. By picking up the most singular term in the right hand side of (52), it follows that

$$\dot{\rho }\sim -6\left(\frac{2}{3}b+b^{″}\right)H\stackrel{⃛}{H}$$

(58)

Then, substituting (55) and (57) for (58), we obtain

$$\tilde{\gamma }=\frac{4}{1-2\beta }$$

(59)

This means that $\rho$ and H evolve around $t=t_s$ as follows:

$$\rho \propto {(t_s-t)}^{\frac{4}{1-2\beta }},H\propto {(t_s-t)}^{\frac{3-2\beta }{1-2\beta }}$$

(60)

Numerically solving the background equations shows that in the presence of quantum corrections, one has $H\propto {(t_s-t)}^{1/3}$ around $t=t_s$, which means that H approaches zero. Meanwhile, in the absence of quantum corrections, we have $H\propto {(t_s-t)}^{-1/3}$, thereby showing the divergence of H at $t=t_s$. From (57), we obtain

$$a\sim a_{0}exp\left[\frac{{\rho }_0^{1-\beta }}{3B\left(1-\beta \right)}{\left(t_s-t\right)}^{\tilde{\gamma }\left(1-\beta \right)}\right]$$

(61)

where $a_0$ is a constant. Comparing the classical case ($\tilde{\gamma }=2/(1-2\beta )$) with the quantum-corrected one ($\tilde{\gamma }=4/(1-2\beta )$), we find that the power of $\left(t_s-t\right)$ is larger in the presence of quantum corrections. Then, the scale factor approaches a constant $a_0$ more rapidly if we account for the quantum effect, implying that the spacetime tends to be smooth, although the divergence of $\rho$ is stronger. Thus, quantum effects moderate the classical singularity.

However, conformal anomaly effects may not always be of huge help in order to alleviate singularities. Take, for example, the case of an asymptotically safe cosmology that was considered in [143]. The capacity to build gravitational RG flow approximations outside of perturbation theory is necessary for conceptually testing asymptotic safety. A very strong framework for performing these calculations is the functional renormalization group equation (FRGE) for the gravitational effective average action ${\mathsf{\Gamma }}_k$:

$${\partial }_k{\mathsf{\Gamma }}_k[g,\overline{g}]=\frac{1}{2}Tr\left[{({\mathsf{\Gamma }}_k^{\left(2\right)}+{\mathcal{R}}_k)}^{-1}{\partial }_k{\mathcal{R}}_k\right]$$

(62)

The construction of the FRGE uses the background field formalism, where the metric $g_{\mu \nu }$ is split into a fixed background ${\overline{g}}_{\mu \nu }$ and fluctuations $h_{\mu \nu }$ (see [144] for more details on asymptotically safe cosmologies). The authors of [143] considered the simplest approximation of the gravitational RG flow, which could be obtained from projecting the FRGE onto the Einstein–Hilbert action approximating ${\mathsf{\Gamma }}_k$ by [144]:

$${\mathsf{\Gamma }}_k=\frac{1}{16\pi G_k}\int d^{4}x\sqrt{-g}\left[-R+2{\mathsf{\Lambda }}_k\right]+\mathrm{gauge-fixing}\mathrm{and}\mathrm{ghost}\mathrm{terms}$$

(63)

where R, ${\mathsf{\Lambda }}_k$, and $G_k$ are the Ricci Scalar, the running cosmological constant, and the running Newton’s gravitational constant, respectively. The scale dependence of these couplings can be written in terms of their dimensionless counterparts as follows:

$${\mathsf{\Lambda }}_k=k^2{\lambda }_{\ast }$$

(64)

$$G_k=g_{\ast }/k^2$$

(65)

where $g_{\ast }=0.707$ and ${\lambda }_{\ast }=0.193$. Considering a background FLRW metric and a perfect fluid for the stress energy tensor $T_{\mu }^{\nu }=\mathrm{diag}[-\rho ,p,p,p]$, one can obtain the Friedmann equation and the continuity equation in this scenario as follows:

$$H^2=\frac{8\pi G_k}{3}+\frac{{\mathsf{\Lambda }}_k}{3}$$

(66)

$$\dot{\rho }+3H(\rho +p)=-\frac{\dot{{\mathsf{\Lambda }}_k}+8\pi \dot{G_k}}{8\pi G}$$

(67)

where the continuity equation comes about from the Bianchi identity satisfied by Einstein’s equations $D^{\mu }[\lambda \left(t\right)g_{\mu \nu }-8\pi G\left(t\right)T_{\mu \nu }]=0$, which usually means that the divergence $D^{\mu }$ of the Einstein tensor vanishes. The extra terms of the right hand side in (67) can be interpreted as an illustration of the energy transfer between gravitational degrees of freedom and matter. Using this new continuity equation, we can write the conformal anomaly term in this case as

$$T_A=-4{\rho }_A-\frac{\dot{{\rho }_A}}{H}\left(\frac{1}{2{\rho }_A}-1\right)$$

(68)

We note that in the conventional cosmology, one could represent the conformal anomaly corrections to the pressure in the form of an integral, but it is clear that this could not be the case for the asymptotically safe cosmology. However, obtaining a corresponding integral for ${\rho }_A$ in Equation (68) is not possible in the same way. Hence, it is not feasible to address a possible removal of Type I–Type III singularities using conformal anomaly effects in this asymptotically safe cosmology.

4.2. Varying Constants Approach

Cosmologies with varying physical constants, like the speed of light or the gravitational constant [145], have been shown to regularize cosmological singularities in certain scenarios [129,146,147,148]. Here, we shall discuss briefly the fundamentals of such theories and how they can be helpful in alleviating both strong and weak cosmological singularities. Examining the generalized Einstein–Friedmann equations within the context of the theories involving a varying speed of light $c\left(t\right)$ (VSL) and varying gravitational constant $G\left(t\right)$ (VG) as presented by Barrow in [145], one can deduce the following expressions for mass density $\varrho \left(t\right)$ and pressure $p\left(t\right)$:

$$\varrho \left(t\right)=\frac{3}{8\pi G\left(t\right)}\left(\frac{{\dot{a}}^2}{a^2}+\frac{k{c}^2\left(t\right)}{a^2}\right)$$

(69)

$$p\left(t\right)=-\frac{c^2\left(t\right)}{8\pi G\left(t\right)}\left(2\frac{\ddot{a}}{a}+\frac{{\dot{a}}^2}{a^2}+\frac{k{c}^2\left(t\right)}{a^2}\right)$$

(70)

These equations highlight the influence of varying c and G on mass density and pressure. For instance, if $\dot{a}$ approaches infinity while $G\left(t\right)$ increases more rapidly than $\dot{a}$, the singularity in $p\left(t\right)$ can be eliminated. In the case of flat models, a direct relationship between pressure p and mass/energy density $\rho /\epsilon$ can be established, albeit with a time-dependent equation of state parameter, expressed as

$$p\left(t\right)=w\left(t\right)\epsilon \left(t\right)=w\left(t\right)\varrho \left(t\right)c^2\left(t\right)$$

(71)

Here, the parameter $w\left(t\right)$ is defined as $w\left(t\right)=\frac{1}{3}\left[2q\left(t\right)-1\right]$, with $q\left(t\right)=-\ddot{a}a/{\dot{a}}^2$ being the dimensionless deceleration parameter. Notably, the variation in the speed of light $c\left(t\right)$ brings about a key distinction between mass density $\varrho$ and energy density $\epsilon =\varrho c^2$, impacting the Einstein mass–energy relationship $E=m{c}^2$, which is transformed here into the mass density–pressure formula $p=\varrho c^2$ after division by volume. The variability of physical constants can be explored through the scale factor, allowing for the examination of scenarios like Big Bang, Big Rip, sudden future, finite scale factor, and w-singularities, as expressed by the scale factor equation

$$a\left(t\right)=a_s{\left(\frac{t}{t_s}\right)}^{m}exp{\left(1-\frac{t}{t_s}\right)}^n$$

(72)

The constants $t_s,a_s,m,$ and n are determined accordingly [129]. This approach illustrates how the varying constant concept aids in regularizing singularities. By inspecting Equations (69) and (70), it becomes evident that a time-dependent gravitational constant variation of the form $G\left(t\right)\propto \frac{1}{t^2}$ eliminates a Type 0 Big Bang singularity in the Friedmann cosmology, addressing both p and $\varrho$ singularities. In Dirac’s scenario [149], where $G\left(t\right)\propto 1/t$, only the $\varrho$ singularity is removed. Moreover, the time dependence of $G=1/t^2$ is less constrained by the geophysical limitations on Earth’s temperature [150].

Another proposal suggests that if the scale factor (72) does not tend to zero as $t\to 0$, it could be rescaled by a “regularizing” factor $a_{rg}=(1+1/t^m)$ ($m\ge 0$), resulting in

$$a_{sm}=\left(1+\frac{1}{t^m}\right){\left(\frac{t}{t_s}\right)}^m={\left(\frac{t}{t_s}\right)}^m+\frac{1}{t_s^m}$$

(73)

Consequently, a varying constant approach (in this case, related to the gravitational constant) can effectively eliminate a strong singularity, such as the Big Bang singularity. A scenario where the varying speed of light contributes to singularity regularization begins by considering a form for the ansatz of $c\left(t\right)$. One common assumption regarding the speed of light’s variation is that it follows the evolution of the scale factor [145]:

$$c\left(t\right)=c_0{a}^s\left(t\right)$$

(74)

With $c_0$ and s as constants, the field equations (69) and (70) can be expressed as

$$\varrho \left(t\right)=\frac{3}{8\pi G\left(t\right)}\left(\frac{{\dot{a}}^2}{a^2}+k{c}_0^2{a}^{2(s-1)}\right)$$

(75)

$$p\left(t\right)=-\frac{c_0^2{a}^{2s}}{8\pi G\left(t\right)}\left(2\frac{\ddot{a}}{a}+\frac{{\dot{a}}^2}{a^2}+k{c}_0^2{a}^{2(s-1)}\right)$$

(76)

In the presence of the time dependence of $c\left(t\right)$ as given by (74), and for the choice of $a\left(t\right)=t^m$, it is possible to eliminate a pressure singularity (Type II) if certain conditions are met: $s>1/m$ for $k=0$, $m>0$, and $s>1/2$ or $m<0$, $s<1/2$ for $k\ne 0$.

4.3. Modified Gravity Effects/Quantum Gravitational Cosmologies

In recent times, there has been wide interest in dark energy models based in exotic non-general relativistic regimes, particularly because such theories display properties that are not evident in conventional cosmological models. For example, many works have considered the possibility of viable scalar-field-based dark energy regimes in quantum-gravity-corrected cosmologies like the RS-II braneworld and loop quantum cosmology [14,15,16,17,18]. There has been substantial work on new dark energy models based on thermodynamic modifications like modified area–entropy relations [151,152,153,154,155,156], even more exotic possibilities like generalized uncertainty principles [157,158,159], or non-canonical approaches like DBI [160,161,162,163,164,165,166]. This vast dark energy literature has prompted the study of cosmological singularities with a wide range of cosmological backgrounds, as there have been multiple works that have discussed Type I–IV singularities in various cosmologies [29,30,32,33,93,143,167,168,169,170,171,172,173,174,175,176,177,178]. In this vast array of literature, one can find quite a few examples where cosmologies affected by these modified gravity theories or quantum gravitational paradigms (like the braneworld or LQC) have alleviated certain singularities. Here, we would like to consider an example of how such effects can help in alleviating Type V singularities, as we have not yet discussed ways to remove or moderate these.

We would like to consider the treatment in [178] for our example here. We would like to again consider a model with an inhomogeneous EOS of the form $p=-\rho -f\left(\rho \right)$. It was shown in [129] that for the simplified case of the scale factor (36) with $m=0$, one can obtain w-singularities for $n>0$ and $n\ne 1$. The scale factor for the case $m=0$ takes the form

$$a\left(t\right)=a_{s}exp{\left(1-\frac{t}{t_s}\right)}^n$$

(77)

We will be using this form of the scale factor for this example. The modified gravity theory we are interested in is an f(R) gravity model with the action [179]

$$S=\frac{m_p^2}{2}\int d^{4}x\sqrt{-g}\left(R-\frac{{\alpha }^2}{R}\right)+\int d^{4}x\sqrt{-g}{\mathcal{L}}_m$$

(78)

where $\alpha$ is a constant that has the units of mass, ${\mathcal{L}}_m$ is the Lagrangian density for matter, and $m_p$ is the reduced Planck’s constant. The field equation for this action is

$$\left(1+\frac{{\alpha }^2}{R^2}\right)R_{\mu \nu }-\frac{1}{2}\left(1-\frac{{\alpha }^2}{R^2}\right)R{g}_{\mu \nu }+{\alpha }^2\left[g_{\mu \nu }{\nabla }_a{\nabla }^a-\nabla {(}_{\mu }{\nabla }_{\nu })\right]R^{-2}=\frac{T_{\mu \nu }^M}{m_p^2}$$

(79)

The Friedmann equation in this case can take the form

$$\frac{6H^2-\frac{\alpha }{2}}{11/8-\frac{8H^2}{4\alpha }}=\frac{\rho }{3}$$

(80)

where $\rho$ is the total energy density. This $F\left(R\right)$ gravity regime was used to explain late-time cosmic acceleration as an alternative to dark energy in [179]. The use of f(R) gravity regimes for avoiding cosmological singularities by adding an $R^2$ term was considered in detail in [180], with the same scenario later being extended in [31,169]. Moreover, based on the properties of the $R^2$ term, non-singular modified gravity was proposed in [181]. The action (79) guides one towards the notion that very tiny corrections to the usual Einstein–Hilbert action in the form of $R^n$ with $n<0$ can produce cosmic acceleration. As corrections of the form $R^n$ with $n>0$ can lead to inflation in the early Universe [182], the authors of [179] proposed a purely gravitational paradigm through (78) to explain both the early and late-time accelerations of the Universe.

Now, we consider $f\left(\rho \right)={\rho }^{\alpha }$ and first examine the status quo of w-singularities for such a model in the standard cosmology given by $H^2=\frac{\rho }{3}$ (written here in natural units for simplicity). We can write the w-parameter for this cosmology as

$$w=-3^{\alpha -1}{\left(\frac{n^2{\left(1-\frac{t}{t_s}\right)}^{2(n-1)}}{t_s^2}\right)}^{\alpha -1}-1$$

(81)

From this, we can make the following observations:

• For n = 1, no w-singularities occur, as is the case in the usual scenario with a conventional equation of state.
• For $\alpha <0$, w-singularities occur for all positive values of n besides unity, but w-singularities do not occur for any negative values of n.
• For $\alpha >0$, we see a very interesting behavior. In this case, completely in contrast to what happens in the usual case, no w-singularities occur for positive values of n ($n>0$), but they occur only when n has negative values ($n<0$). Hence, here we see the first sign of departure in the occurrence conditions of w-singularities when one considers inhomogeneous equations of state.

Thus, we see here that incorporating an inhomogeneous EOS can be of use in moderating w-singularities, but this still does not remove them per say as it only changes the conditions under which they occur with regards to what happens in the conventional cosmology. Now, the w-parameter for the case in (80) is given by

$$w=-{12}^{\alpha -1}{\left(\frac{\eta \left(12n^2{\left(1-\frac{t}{t_s}\right)}^{2n}-\eta {(t_s-t)}^2\right)}{11\eta {(t_s-t)}^2-18n^2{\left(1-\frac{t}{t_s}\right)}^{2n}}\right)}^{\alpha -1}-1$$

(82)

For the w-parameter expressed above, we have the following observations:

• For n = 1, contrary to the other cases we have considered, a w-singularity can occur, but this is possible only in the extreme case that $\alpha \to \infty$. This cannot realistically be expected, but in principle singularities can appear in this case.
• The most interesting detail that arises when one considers this scenario is that w-singularities do not occur for any value of n and $\alpha$. For both positive and negative values of $\alpha$ and n, the w-parameter remains regular and does not diverge.

Thus, we see that just by incorporating the effects of a modified gravity theory, in this case a particular form of $f\left(R\right)$ gravity, one can also alleviate singularities. Furthermore, $f\left(R\right)$ gravity theories have been of great use in alleviating various other singularities, which we have discussed quite extensively; hence, it seems appropriate to discuss this example to illustrate how Type V singularities could be moderated too.

5. Dynamical Systems Approach and the Goriely–Hyde Method

While it is seems quite natural to study singularities and their avoidance methods in various cosmological settings, as we have discussed so far, often it is very difficult to classify and study the cosmological singularities that may occur in extremely non-conventional cosmologies motivated by quantum gravitational/phenomenological considerations (for example, see the classification of singularities in asymptotically safe cosmology [143]), and this may not even be possible in an orthodox fashion. Hence, it becomes essential to look for non-conventional ways to find cosmological singularities in exotic cosmologies, and, in this regard, a particular dynamical systems method can be of huge help. From a dynamical standpoint, one of the most intriguing aspects of studying various dynamical systems lies in understanding their singularity structure, which becomes particularly relevant when these systems describe physically significant phenomena. While numerous approaches have been proposed to explore the singularity structure of autonomous dynamical systems, one particularly interesting method is the Goriely–Hyde procedure [183]. As cosmology presents a multitude of captivating dynamical systems [184], the investigation of singularity structure in such systems has gained considerable attention, with the Goriely–Hyde method proving particularly useful for cosmological explorations [185,186,187,188,189,190]. This method has previously been applied to study finite- and non-finite-time singularities in certain classes of quintessence models as well [32,170,191].The Goriely–Hyde method provides an elegant approach to determining the presence of singularities in dynamical systems, and the procedure can be outlined as follows:

• We begin by considering a dynamical system described by n differential equations of the form

  $${\dot{x}}_i=f_i\left(x\right),$$

  (83)

  where $i=1,2,\dots ,n$, and the overdot represents differentiation with respect to time t, which in the case of quintessence models can be better represented by the number of e-foldings N. We identify the parts of the equation $f_i$ that become significant as the system approaches the singularity. These significant parts are referred to as “dominant parts” [183]. Each dominant part constitutes a mathematically consistent truncation of the system, denoted as ${\widehat{f}}_i$. The system can then be written as

  $${\dot{x}}_i={\widehat{f}}_i\left(x\right).$$

  (84)
• Without loss of generality, the variables $x_i$ near the singularity can be expressed as

  $$x_i=a_i{\tau }^{p_i},$$

  (85)

  where $\tau =t-t_c$, and $t_c$ is an integration constant. Substituting Equation (4) into Equation (3) and equating the exponents, we can determine the values of $p_i$ for different i values, which form the vector $\mathbf{p}=(p_1,p_2,\dots ,p_n)$. Similarly, we calculate the values of $a_i$ to form the vector $\overrightarrow{a}=(a_1,a_2,\dots ,a_n)$. It is important to note that if $\overrightarrow{a}$ contains only real entries, it corresponds to finite-time singularities. Conversely, if $\overrightarrow{a}$ contains at least one complex entry, it may lead to non-finite-time singularities. Each set $(a_i,p_i)$ is known as a dominant balance of the system.
• Next, we calculate the Kovalevskaya matrix given by

  $$R=\left(\begin{array}{ccccc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_1}{\partial x_2}& .& .& \frac{\partial f_1}{\partial x_n}\\ \frac{\partial f_2}{\partial x_1}& \frac{\partial f_2}{\partial x_2}& .& .& \frac{\partial f_2}{\partial x_n}\\ .& .& .& .& .\\ .& .& .& .& .\\ \frac{\partial f_n}{\partial x_1}& \frac{\partial f_n}{\partial x_2}& .& .& \frac{\partial f_n}{\partial x_n}\end{array}\right)-\left(\begin{array}{ccccc}{p}_1& 0& .& .& 0\\ 0& p_2& .& .& 0\\ .& .& .& .& .\\ .& .& .& .& .\\ 0& 0& .& .& p_n\end{array}\right).$$

  (86)
  After obtaining the Kovalevskaya matrix, we evaluate it for different dominant balances and determine the eigenvalues. If the eigenvalues are of the form $(-1,r_2,r_3,\dots ,r_n)$, with $r_2,r_3,\cdots >0$, then the singularity is considered general and will occur regardless of the initial conditions of the system. Conversely, if any of the eigenvalues $r_2,r_3,...$ are negative, the singularity is considered local and will only occur for certain sets of initial conditions.

After applying the method, one can then classify singularities using well-supported ansatzes for the scale factor or the Hubble parameter. The most general form of the Hubble parameter for investigating singularities within the aforementioned classified types is expressed as [191]

$$H\left(t\right)=f_1\left(t\right)+f_2\left(t\right){(t-t_s)}^{\alpha }$$

(87)

Here, $f_1\left(t\right)$ and $f_2\left(t\right)$ are assumed to be non-zero regular functions at the time of the singularity, and similar conditions apply to their derivatives up to the second order. Additionally, $\alpha$ is a real number. It is not mandatory for the Hubble parameter (34) to be a solution to the field equations; however, we will consider this case and explore the implications of this assumption on the singularity structure based on our dynamic analysis. First, we observe that none of the variables x, y, or z as defined in (10) can ever become singular for any cosmic time value. The singularities that can occur considering the Hubble parameter as defined in (34) are as follows:

• For $\alpha <-1$, a Big Rip singularity occurs.
• For $-1<\alpha <0$, a Type III singularity occurs.
• For $0<\alpha <1$, a Type II singularity occurs.
• For $\alpha >1$, a Type IV singularity occurs.

Another ansatz useful for classifying singularities was introduced in [38], where the scale factor was written as

$$a\left(t\right)=g\left(t\right){(t-t_s)}^{\alpha }+f\left(t\right)$$

(88)

Here, g(t) and f(t) and all their higher-order derivatives with respect to the cosmic time are smooth functions of the cosmic time. For this ansatz, according to the values of the exponent $\alpha$, one can have the following singularities:

• For $\alpha <0$, a Type I singularity occurs.
• For $0<\alpha <1$, a Type III singularity develops.
• For $a<\alpha <2$, a Type II singularity occurs.
• For $\alpha >2$, a Type IV singularity occurs.

Again, it is not mandatory that the scale factor in (88) be a solution to the field equations, but we would like to consider this and (87) in order to obtain a well-motivated estimation of the type of cosmological singularities we can deal with in the various models we have discussed so far.

As an example of this method, let us consider singularities in an RS-II braneworld cosmology where dark energy can be described by a scalar field paradigm, following the treatment of [32]. The action for the inclusion of both the scalar and the background fluid term can be written as

$$\begin{array}{c}S=S_{RS}+S_B+S_{\varphi }=\int d^{5}x\sqrt{-g^{\left(5\right)}}\left({\mathsf{\Lambda }}^{\left(5\right)}+2R^{\left(5\right)}\right)+\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\int d^{4}x\sqrt{-g}\left(\sigma -\frac{1}{2}\mu \left(\varphi \right){(\nabla \varphi )}^2-V\left(\varphi \right)+{\mathcal{L}}_B\right)\end{array}$$

(89)

where $R^{\left(5\right)}$, $g_{\mu \nu }^{\left(5\right)}$, and ${\mathsf{\Lambda }}^{\left(5\right)}$ are the bulk Ricci Scalar, metric, and cosmological constant, respectively, with $\sigma$ being the brane tension on the 3-brane, $g_{\mu \nu }$ being the 3-brane metric, and $\mu \left(\varphi \right)$ being a scalar coupling function. Note that here we are working in Planck units with ${(m_p^{(5)})}^2=1$ and $m_p^{\left(5\right)}$ being the five-dimensional Planck mass. Assuming that the brane metric has the usual FLRW form, we obtain the Friedmann equation as follows: [192]

$$H^2=\rho \left(1+\frac{\rho }{2\sigma }\right)$$

(90)

where $\rho ={\rho }_{\varphi }+{\rho }_B$ is the total cosmological energy density taking into account contributions from both the scalar field and the background fluid term, and the bulk cosmological constant has been set to zero for simplicity. One can similarly find that

$$2\dot{H}=-\left(1+\frac{\rho }{\sigma }\right)\left(\mu \left(\varphi \right){\dot{\varphi }}^2+{\rho }_B\right)$$

(91)

The equation for the motion of the scalar is given by

$$\mu \left(\varphi \right)\ddot{\varphi }+\frac{1}{2}\frac{d\mu }{d\varphi }{\dot{\varphi }}^2+3H\mu \left(\varphi \right)\dot{\varphi }+\frac{dV}{d\varphi }=0$$

(92)

Finally, using the variables introduced in [193],

$$x=\frac{\dot{\varphi }}{\sqrt{6}H}y=\frac{\sqrt{V}}{\sqrt{3}H}z=\frac{\rho }{3H^2}$$

(93)

and setting the background fluid to have the form of pressureless dark matter, such that $w_B=0$, we obtain the dynamical system for this model as follows:

$$x^{\prime }=-\sqrt{\frac{3}{2\mu }}\lambda y^2-3x+\frac{3x}{2}\left(z+x^2-y^2\right)\left(\frac{2}{z}-1\right)$$

(94)

$$y^{\prime }=\sqrt{\frac{3}{2\mu }}\lambda xy+\frac{3y}{2}\left(z+x^2-y^2\right)\left(\frac{2}{z}-1\right)$$

(95)

$$z^{\prime }=3(1-z)(z+x^2-y^2)$$

(96)

where the primes denote differentiation with respect to the e-folding number N, and $\lambda =\frac{V^{\prime }}{V}$. We can finally start with the analysis as we have a proper autonomous dynamical system; the first truncation that we consider is

$$\widehat{f}=\left(\begin{array}{c}-k\lambda y^2\\ -3y^3{z}^{-1}\\ 3x^2\end{array}\right)$$

(97)

where $k=\sqrt{\frac{3}{2\mu }}$. Using the ansatz of the Goriely–Hyde method, we obtain $\mathbf{p}=(-1,-1,-1)$, and, using these, we obtain

$$\begin{array}{c}\hfill a_1=\left(-\frac{1}{k\lambda },\frac{i}{k\lambda },-\frac{3}{k^2{\lambda }^2}\right)\\ \hfill a_2=\left(-\frac{1}{k\lambda },-\frac{i}{k\lambda },-\frac{3}{k^2{\lambda }^2}\right)\end{array}$$

(98)

As both $a_1$ and $a_2$ have complex entries, only non-finite-time singularities will be possible with regards to this truncation. The Kovalevskaya matrix then takes the form

$$R=\left(\begin{array}{ccc}1& -2k\lambda y& 0\\ 0& 1-\frac{9y^2}{z}& \frac{3y^3}{z^2}\\ 6x& 0& 1\end{array}\right)$$

(99)

We then finally find the eigenvalues of the matrix, which are given by

$$r=(-1,-1,2)$$

(100)

Hence, the singularities in this case will only be local singularities, which will only form for a limited set of initial conditions. In [32], it was determined that there are two more possible truncations. The balances and corresponding eigenvalues for the first of these are, respectively,

$$\begin{array}{c}\hfill a_1=\left(\frac{1}{\sqrt{3}},\frac{i}{\sqrt{3}},\frac{1}{3}\right)\\ \hfill a_2=\left(\frac{1}{\sqrt{3}},-\frac{i}{\sqrt{3}},\frac{1}{3}\right)\\ \hfill a_3=\left(-\frac{1}{\sqrt{3}},\frac{i}{\sqrt{3}},\frac{1}{3}\right)\\ \hfill a_4=\left(-\frac{1}{\sqrt{3}},-\frac{i}{\sqrt{3}},\frac{1}{3}\right)\end{array}$$

(101)

and

$$r=(-1,\sqrt{\frac{3}{2}},-\sqrt{\frac{3}{2}})$$

(102)

The other truncation has the balances

$$\begin{array}{c}\hfill a_1=\left(\frac{1}{\sqrt{3}},\sqrt{\frac{2}{3}},\frac{2}{3}\right)\\ \hfill a_2=\left(\frac{1}{\sqrt{3}},-\sqrt{\frac{2}{3}},\frac{2}{3}\right)\\ \hfill a_3=\left(-\frac{1}{\sqrt{3}},\sqrt{\frac{2}{3}},\frac{2}{3}\right)\\ \hfill a_4=\left(-\frac{1}{\sqrt{3}},-\sqrt{\frac{2}{3}},\frac{2}{3}\right)\end{array}$$

(103)

with

$$r=(-1,-1,1)$$

(104)

We see from (101) and (102) that the truncations to which they belong seem to still suggest that only non-finite-time local singularities are possible in the system, but we note from (103) that the other truncation will allow for finite-time singularities, albeit they would still be local as (104) has $r_2=-1$. To proceed further and now classify the singularities physically, we use the ansatz for the Hubble parameter (87), and we need to express $\dot{\varphi }$ and $V\left(\varphi \right)$ in terms of the Hubble parameter. For simplicity, we will consider that the coupling constant $\mu =1$ and $\dot{{\rho }_B}=0$. With these considerations, we can write

$$-2\dot{H}={\dot{\varphi }}^2\left(1+\frac{\rho }{\sigma }\right)$$

(105)

One can then write

$${\dot{\varphi }}^2=-2\left[\left(\sigma +V+\sigma {\rho }_B\right)+\sqrt{{\left(\sigma +V+\sigma {\rho }_B\right)}^2-2\dot{H}}\right]$$

(106)

Furthermore, one can now write $V\left(\varphi \right)$ in terms of the dark energy equation of state as

$$V\left(\varphi \right)=\frac{{\dot{\varphi }}^2}{2}\frac{(1-w)}{(1+w)}$$

(107)

Note that here we are only considering the dark energy equation of state with no background contributions; hence, here we will only consider scalar field contributions. Then, we can write the potential as

$$V=\frac{2b(1+k)+\sqrt{{\left(2b(1+k)\right)}^2-2\dot{H}(k^2-1)}}{2(k^2-1)}$$

(108)

where $k=\frac{2w}{1-w}$ and $b=\sigma (1+{\rho }_B)$ (note that both k and b will always be positive for a positive brane tension). Notice that V is now completely in terms of the Hubble parameter (for constant values of $\sigma$, w, and ${\rho }_B$), and so one can use this form of V to find $\dot{\varphi }$ in terms of the Hubble parameter as well. It is necessary to express these quantities in terms of $H\left(t\right)$ as now we can find out which type of singularities are possible in this scenario, in view of the fact that x, y, and z have to remain regular. By studying the expressions for these variables, one can determine that Type I, Type III, and Type IV singularities are allowed in this scenario, while Type II singularities are not. This also makes us realize that even if the cosmology is heavily motivated by quantum gravitational considerations (like the braneworld in this scenario), it can still have quite a few cosmological singularities.

6. Future Outlook and Conclusions

In this brief review paper, we discussed (almost) all the prominent developments in the field of cosmological singularities that have taken place in the past 25 years or so. We firstly provided a detailed outlook on what spacetime singularities are and discussed their various nuances, like their various strength criteria. After this, we discussed in detail the prominent strong and weak cosmological singularities in accordance with the classification scheme provided by Odintsov and Nojiri. We detailed under which conditions these singularities can occur in various scenarios and under which cosmological settings they were initially discovered. We then saw how one can moderate or even remove these singularities using various techniques with quantum or modified gravity origins, and we also discussed the Goriely–Hyde method and its usefulness in singularity works. As a whole, one general point that we can safely make is that such singularities provide a revealing arena for the interface of cosmology and quantum gravitational theories. The scales at which such events could take place lie on the horizon for testing quantum gravity ideas, and, with the constant increase in the precision of various observational setups like the WST, CTA, SKA, Euclid, LSST, and Roman Space Telescope [194,195,196,197,198,199], one would not be wrong to think that investigating such singularities in detail could possibly shed light on problems in both cosmology and quantum gravity. Besides this, we also saw that cosmological singularities may have subtle connections to current cosmological tensions like the H0 crisis.