Recent advances in cosmological singularities

The discovery of universe's late-time acceleration and dark energy has overseen 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 2 decades. We discuss the fundamentals of space-time singularities after which we discuss about all the different forms of cosmological singularities which have been discovered in recent times in detail. We then talk about methods and techniques to avoid or moderate these singularities in various theories and discuss how these singularities can occur in non-conventional cosmologies too. 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.


Introduction
Observations of 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 done in order to explain this expansion. The cosmological expansion problem has been addressed from multiple facets till now, which include the standard approaches of the Cosmological constant [2][3][4] alongside more exotic scenarios like Modified gravity theories [5][6][7], scalar field driven late-time cosmic acceleration scenarios [8][9][10][11][12][13][14]. 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 [15][16][17][18][19][20][21][22][23][24][25]. This, however, has also sprung up some discrepancies which seem to be pointing towards the limits of our current understanding of the universe, most famous of which is arguably the Hubble tension which refers to the disagreements between the values of the Hubble constant measured from. Detailed Cosmic Microwave Background Radiation(CMB) maps, combined with Baryon Acoustic Oscillations data, and those from the SNeIa data [26][27][28]. Hence, the current epoch of the universe has certainly provided us with a wide range of questions and looks set to become an avenue where advanced gravitational physics will lead the way towards better understanding of cosmology.
There has also been an expansive literature in recent times which has been devoted to the study of various types of singularities that could occur during the current and far future of the Universe, with the observation of late-time acceleration having given a significant boost to such works .Even the term singularity comprises many different definitions and with regards to cosmological cases, until the end of the 20th century, the only popular possibilites of singularty 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 proven by them make use of the null energy condition and also of the facts that at a singular point of the spacetime, geodesics incompleteness occurs and also the curvature scalars diverge. Although in modified gravity the null energy condition may be different in general in comparison to the Einstein-Hilbert case (see for example), it is generally accepted that the geodesic incompleteness and also that 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 was performed in . While one can treat singularities as points at which a cosmological theory somewhat fails, one might also consider them as windows to new physics and thus have a different kind of appealing interest with it. Particularly finite-time singularities (those which happen in a finite time) could be viewed as either flaws the classical theory, or alternatively as a doorway towards a quantum description of general relativity. This is due to the fact that these cannot be dressed in a similar way to the spacelike singularities of black holes for instance, and so one is left to ponder about the accuracy of the predictions of classical gravitational theories. Hence studying singularities in cosmological contexts and how they could be (possibly) removed, provides a way towards a deeper understanding of relating quantum descriptions of cosmology with its classical ones.
These cosmological singularities which have been discussed in recent times can be classified broadly into two types ; strong and weak (such a classification was initially put forward by [52]). Strong singularities are those singularities which 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" signalling the end of the universe. Weak singularities, as the name might suggest, are those which do not have such far reaching implications and do not distort finite objects in the same sense as their strong counterparts. We can discuss these various singularities in more detail as follows, in accordance to the classification provided in [32,53] : • Type -1 ("Grand Bang/ Grand rip") : In this scale the scale factor becomes null (Bang) or diverges (rip ) for w = −1 [54] • Type 0 ("Big Bang") : In this case, the scale factor becomes null for w ̸ = −1 • Type I ("Big Rip") : In this case, the scale factor , effective energy density and effective pressure density diverges for w ̸ = −1. This is a scenario of universal death, wherein everything which resides in the universe is progressively torn apart [55].
• Type II ("Sudden/Quiescent singularity") : In this case, the pressure density diverges and so does the derivatives of the scalar factor from the second derivative onwards [56]. The weak and strong energy conditions hold for this singularity. Also known as quiescent singularities, but this name originally appeared in contexts related to non-oscillatory singularities [57]. A special case of this is the big brake singularity [58].
• Type III ("Big Freeze") : In this case, the derivative of the scale factor from the first derivative on wards diverges. These were detected in generalized Chaplygin gas models [59].
• 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 [32,60].
• Type V ("w-singularities") : In this case, the scale factor, the energy and pressure densities are all finite but the barotropic index w = p ρ becomes singular [61].
• Type ∞(" Directional singularities "): Curvature scalars vanish at the singularity but there are causal geodesics along which the curvature components diverge [62] and in this sense, the singularity is encountered just for some observers.
• Inaccessible singularities: These singularities appear in cosmological models with toral spatial sections, due to 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 and hence this prompts the name inaccessible singularities [63].
• Type -1 ("Grand Bang/Rip") : In this case, the scale becomes null or diverges for w = −1 [64]. All of these 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 quo of such singularities. In section II, we would discuss in detail all the singularities we have listed above and how these have been shown to form in different cosmological scenarios. In Section III, we will discuss some lesser known singularities which are more special cases of the previously listed singularities while in section IV, we would discuss various ways which have been shown to remove such singularities (in some cases). In Section V we would discuss a particular dynamical system analysis method (known as the Goriely-Hyde method ) which has been shown to be very useful for cosmological singulartity discussions. Finally in setion VI, we would summarize our brief review and discuss the future outlooks for cosmology with regards to singularities.

An overview on space-time 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. These models 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 spacetimes 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 space-time models described by general relativity. Examples include the Friedmann-Robertson-Walker (FRW) cosmological models and the Schwarzschild space-time. These models exhibit singularities where energy density and curvatures become infinitely large, leading to a breakdown of the conventional description of space-time.The Schwarzschild spacetime displays an essential curvature singularity at r = 0, where the Kretschmann scalar α = R ijkl R ijkl diverges along any non-spacelike trajectory approaching the singularity. Similarly, for FRW models with ρ + 3p > 0 at all times (where ρ 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 ρ 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 space-time, where the usual laws of physics fail. Geodesic incompleteness implies that a timelike observer will cease to exist in the space-time after a finite amount of proper time.While singular behavior can occur without extreme curvature, such cases are considered artificial. An example is the Minkowski space-time 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 space-time to be "inextendible," meaning it cannot be isometrically embedded into a larger space-time 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 ds 2 = −dt 2 + dr 2 + r 2 (dθ 2 + sin 2 θ dϕ 2 ) with the identification ϕ = 0 and ϕ = a (with a ̸ = 2π), 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 space-time reveals that singularities likely exist in a broad range of space-times, 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 Schwarzschild and Friedmann. However, it remains possible that such a divergence only occurs at infinity for a given nonspacelike path. In a broader sense, it seems reasonable to expect some form of curvature divergence to occur along nonspacelike 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 nonspacelike 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 nonspacelike 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 nonspacelike geodesic or an inextendible nonspacelike 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 : ds which remains finite. The components X i represent 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 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 nonspacelike curve, the singularity might be removable by extending the spacetime with relaxed differentiability requirements [65].
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 nonspacelike curve. Conversely, a "scalar polynomial singularity" occurs when a scalar polynomial involving the metric and Riemann tensor takes on infinitely large values along a nonspacelike 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 [52]. 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 [66] 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 [67] 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 [68]. Working with Jacobi fields can be demanding as it involves solving the Jacobi equation along geodesics. Nevertheless, conditions for lightlike and timelike geodesics, satisfying both criteria, have been established [65]. These conditions are expressed in terms of 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 diverges as the normal parameter τ approaches the singularity. Krolak's criterion states that the singularity is strong if and only if the integral τ 0 dτ ′ R ij u i u j diverges as τ approaches the singularity.
• Timelike geodesics: For timelike geodesics, [65] presents various necessary and sufficient conditions, but not a single characterization.
Adhering to Tipler's criterion, a singularity is strong along a timelike geodesic if the integral diverges on approaching the singularity.
Conforming to Krolak's criterion, the singularity is strong if the integral τ 0 dτ ′ R ij u i u j diverges on approaching the singularity. Additional necessary conditions exist, although they are not utilized for our purposes.
In passing, it is also of interest to talk of the cosmic censorship conjecture [69], 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 and has sprung up a lot of work the topic of naked singularities. In principle, one could think of cosmological singularities as naked singularities as well given that there is no need of an event horizon in such cases and there is no development of such horizon in cases in which such singularities develop too. Several examples of spacetimes containing naked singularities have been found in recent times [70][71][72][73][74][75][76][77][78][79][80][81][82]. When such singularities develop in gravitational collapse, they give rise again to extremely intriguing physical possibilities and problems. The opportunity offered in that case is, we may have the possibility to observe the possible ultra-high energy physical processes occurring in such a region of 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 favour of their existence [83,84]. 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 can also be a good avenue for testing out the predictions of such theories. Very recently another very interesting work has been done with regards to naked singularities was performed in [85]. The authors there a performed general relativistic ray-tracing and radiative transfer simulations to generate synchrotron emission images utilising thermal distribution function for emissivity and absorptivity. They investigate effects in the images of JMN-1 naked singularity and Schwarzschild black hole by varying inclination angle, disk width and frequency. Their results give further motivation of naked singularities being a realistic scenario.

Strong singularities
As mentioned before, strong singularities are those singularities which can distort finite objects in space-time and now we would like to discuss of the prominent singularities in this regard.

Big Bang singularity(Type 0)
Classical models of the universe generically feature an initial or 'big bang' singularity. This is when we consider progressively earlier and earlier stages of the universe, observable quantities stop behaving in a physically reasonable way. A more precise mathematical characterisation of the cosmic big bang singularity can be made in terms of both a global notion of 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 seeing a singularity in a finite proper time. The idea of the big bang hence, as is popularly known is, the singularity at the very beginning of the universe.

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, the energy density and the pressure density of the universe also becomes 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 [55], which means that the equation of state 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 [86]. 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 [87][88][89][90]. Phantom constructions have also been discussed in the context of quantum gravitational theories, for example in various string theoretic realizations of dark energy. [15,91,92]. So it seems in principle interesting to look for a late time universe scenario with phantom dominance and that is where 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's even more intriguing is the potential of 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 (ρ) and three times the pressure (p), denoted as ρ + 3p. For instance, a planet in orbit around a star with mass M and radius R becomes unbound approximately when the condition −(4π/3)(ρ + 3p)R 3 ≈ M is satisfied. In cases where the equation −(ρ + 3p) decreases over time due to a parameter w greater than or equal to -1, if −(4π/3)(ρ + 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 −(ρ + 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 ≈ P 2|1 + 3w|/[6π|1 + w|]. 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 time being 35 billion years. Big rip even rips apart molecules, 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 explored way to avoid the big rip too (for example [93]) and we will be discussing those later on here. A lot more works have been done on various aspects of big rip singularities over the years, see [94][95][96][97][98][99][100][101][102][103][104][105][106]. The comparison of Big Rip from Dark energy and modified gravity was done for the first time in [107]. Less drastic variants of the big rip have also been found in recent years [108][109][110] which have been discussed in detail in Appendix B.

Grand bang and Grand rip singularities (Type −1)
The grand bang, although apparently different in name, is quite the same as the big bang singularity with null scale factors and a diverging pressure and energy density with the one difference being that the singularity occurs with the equation of state parameter being equal to -1 [64]. This type of singularity was found initially by using a series ansatz for the scale factor. An understanding of the grand bang and grand rip singularities is quite intricately linked with each other and so we shall discuss that now.
To discuss these grand singularities, we note that the equation for the parameter w is given by: 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: assuming flat models. Otherwise, the relationship between these parameters becomes more intricate, involving the Hubble parameter H =ȧ/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(t). It is advantageous to introduce the variable x = ln a: allowing us to define as a correction around the case of a pure cosmological constant: This change of variables assists in reducing the order of the differential equation: In terms of this, one finds the scale factor to be If one then assumes a power series form of h(t) ( which has become quite a well motivated ansatz for the scale factor in various cosmological studies, see ) It can then be found out that the energy and pressure densities can be written as and the pressure, This presents us with intriguing possibilities, where our specific focus will be on the case where η 0 > 0. In this scenario, we observe that at t = 0, ρ and p exhibit divergences following t −2(η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 [111] and [112]. This is due to the behavior of the scale factor (which is an exponential of rational functions); it doesn't lend itself to convergent power expansions, whether generalized or not, with a finite number of terms featuring negative powers. However, the function x(t) 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(t) as Based on this, we make the following observations: • For h 0 > 0: In this scenario, 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, in this case, 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.

Directional singularities (Type ∞)
We follow the discussion of [62] in order to understand how directional singularities were initially found in cosmological models. If we start with the flat FLRW metric in the form indicated by Eq. (1): it is evident that the equations governing the trajectories of geodesics, followed by observers not subject to acceleration (δ = 1) and light-like particles (δ = 0) possessing a specific linear momentum P , can be simplified to: assuming constant θ and ϕ due to the symmetry inherent in these models.
Here, τ represents the intrinsic or proper time as measured by the observer.In the context of null geodesics, we find that: Consequently, to ensure that the initial event t = −∞ corresponds to a finite proper time interval ∆τ from an event at t, the requirement is: Therefore, the emergence of singular behavior exclusively at t = −∞ is possible if the scale factor can be expressed as an integrable function of coordinate time. This condition necessitates that a(t) tends towards zero as t approaches −∞, although this alone is not sufficient. Similarly, for timelike geodesics with nonzero P : indicating that the proper time interval to t = −∞ is finite provided the time interval for light-like geodesics is also finite. Consequently, t = −∞ is reachable for these observers. As a result, condition (7) implies that both light-like and timelike geodesics with non-zero P experience t = −∞ within a finite proper time interval in their past. Conversely, comoving observers tracing timelike geodesics with P = 0 exhibit dτ = dt, which leads to t = −∞ corresponding to an infinite proper time interval in their past, and 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: These differ from the "little rip" model in the sign of h(t), and are termed "little bang" if they denote an initial singularity, or "little crunch" if they represent a final singularity [113]. Instances of this case encompass models with a scale factor a(t) ∝ e −α(−t) p where p > 1 and α > 0. • Changing the sign of h(t) gives rise to a variant of the "little rip" scenario, featuring an asymptotically vanishing energy density and pressure. Models with a scale factor a(t) ∝ e −α(−t) p where p ∈ (0, 1) and α > 0 exemplify this case.
This case applies to models like a(t) ∝ t −p with p > 1, as explored in [62].
While they have recently been discussed in the context of inflationary models [113],not much work has been done in the regard of Type ∞ singularities since their discovery with regards to their avoidance or their emergence in more exotic cosmological models.

Sudden singularities (Type II)
In the case of such type II singularities, one has the pressure density diverging or equivalently, the derivatives of the scale factor diverging from the second derivative onwards. Let's start by examining informally whether there is a 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(t) approaches a non-zero or infinite value a(t s ) and the Hubble parameter H(t) 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: 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, manifested as: as t → t s , consistent with the conditions outlined in Equation (10). In such instances, the pressure singularity is concomitant with an infinite acceleration.
To illustrate an example for this, we take the most primitive example of such singularities which was put forward by Barrow in [56]. In this regard, assume that it is physically reasonable to expect that the scale factor can be written in the form of this ansatz 1 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(0) = 0, leading to A = −Ct n s > 0. Consequently, we find the expression for the Hubble parameter H s : For simplicity, we use the freedom to scale the Friedmann metric by dividing by A and set A ≡ 1 and C ≡ −t n s . This yields the simplified form of a(t): where a s ≡ a(t s ). As t approaches t s from below, the behavior of the second derivative of a can be described: whenever 1 < n < 2 and 0 < q ≤ 1. This solution is valid for 0 < t < t s . Consequently, as t approaches t s , a approaches a s , H s and ρ s > 0 (as long as 3q 2 (a s − 1) 2 t −2 s > −k) remain finite, while p s → ∞. When 2 < n < 3,ä remains finite but ... a → ∞ as t approaches t s . Here, p s remains finite, butṗ s → ∞.In contrast, there exists an initial strong-curvature singularity, where both ρ and p tend to infinity as t approaches 0. Importantly, in this scenario, both ρ and ρ + 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 [56], such singularities have been discussed in loads of various different settings, both from modified gravity perspectives and other phenomenological considerations. Work has also been done on ways to escape such singularities and we will be discussing them later in our discussions.

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 firstly shown in a phantom generalized Chaplygin gas cosmology (PGCG) in [59] 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: where the symbol A represents a positive constant and α signifies a parameter. In the scenario where α = 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:ρ + 3H(p + ρ) = 0, from which emerges the expression for energy density ρ: where B stands as a constant parameter. In a noteworthy observation made in Ref. [98], it was discerned that a negative B renders the perfect fluid, with the equation of state p = − A ρ α , unable to uphold the null energy condition, that is, p + ρ < 0. Intriguingly, under these conditions, the energy density rather escalates as the Universe expands, contrary to redshift behavior, thus earning the term "phantom generalized Chaplygin gas" (PGCG).
Further insights from the works of [98,114] unveil that for PGCG with α > −1, a FLRW Universe hosting this fluid can evade the impending big rip singularity. As the scale factors attain magnitudes far beyond, the Universe eventually approximates an asymptotically de Sitter state. In stark contrast, during the Big Freeze scenario, the PGCG energy density responds by amplifying as the scale factor matures. Especially, as the scale factor approaches minuscule values (a → 0), ρ tends towards A 1 1+α , while it experiences a surge at a finite scale factor a max : .
As a consequence, a 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: Remarkably, this singularity not only emerges at a finite scale factor but also at a distinct future cosmic time. Conversely, the history of a FLRW Universe permeated with this fluid traces back to an asymptotically de Sitter state in the past. Expressing this temporal journey succinctly: where a 0 signifies a minute scale factor. Additionally, the Universe embarks on its odyssey from a bygone infinity of cosmic time, as a → 0 and p + ρ → 0 − . Remarkably, the homogeneous and isotropic nature of the Universe propels it into a phase of super-accelerated expansion, denoted by: up until it culminates at the singularity a = a max . It's imperative to recall that the PGCG eludes satisfaction of the null energy condition [98], as embodied in p + ρ < 0. A great amount of work has been done on Big Freeze singularities since the intial one by in [59] and it has been shown that one can encounter such singularities in a lot of exotic cosmological settings and also there have been works probing how one can avoid such singularities too [115][116][117][118][119][120][121][122].

Generalized sudden singularities (Type IV)
These singularities were firstly discussed in [31] and have since then been found in a diverse variety of cosmological settings. So here we will briefly discuss the primary cases in which type IV singularities were shown. In fact, this example will illustrate all the prominent singularities we have discussed so far. We start with an equation of state of the form This sort of equation of state with f (ρ) = Aρ α for α being an arbitrary constant was first proposed in [60] and was investigated in detail in [94] and there can be diverse physical motivations behind such an equation of state. This form of EOS can also be equivalent to bulk viscosity [123]. This type of an equation of state can also come about due to modified gravity effects [30]. We now consider the following ansatz for the scale factor where n is a positive constant and 0 < t < t s . The scale factor diverges within a finite time (t → t s ), resembling the phenomenon of the Big Rip singularity. Consequently, t s represents the universe's lifetime. When t ≪ t s , the evolution of a(t) follows t n , leading to an effective EOS given by w = −1 + 2/(3n) > −1.
The Hubble rate in this case can be expressed as Utilizing Equation (19) one can deduce the relation As a result, both H and ρ exhibit minima at t = t s /2, characterized by the values Next, we examine a specific form for f (ρ) given by where A, B, α, and β are constants. As we shall see, this dark energy scenario harbors a complex structure with respect to singularities. In scenarios where α surpasses β, we observe that For non-unit values of α and β, we obtain The realm of possibilities in this cosmology is extensive. If 1 > α > β and A, B > 0 (A, B < 0), the scale factor has a minimum (maximum) at ρ = 0, extending to infinity (vanishing) as ρ → ∞. When α > 1 > β and A < 0 while B > 0 (A > 0 and B < 0), the scale factor features a minimum (maximum) at a non-trivial (non-vanishing) ρ value, reaching infinity (zero) as ρ approaches zero or a positive infinity. For α > 1 > β and A, B > 0 (A, B < 0), the scale factor becomes infinite (vanishes) as ρ → ∞ (ρ → 0), and it vanishes (increases) as ρ → 0 (ρ → ∞). When α > β > 1, the scale factor approaches a 0 as ρ → ∞. Additionally, if A > 0 (A < 0), the scale factor tends to 0 (∞) as ρ → 0. With A, B > 0 (A, B < 0), the scale factor demonstrates monotonic growth (decrease) concerning ρ. In the case of A > 0 and B < 0 (A < 0 and B > 0), the scale factor attains a nontrivial maximum (minimum) at a finite ρ value.
To summarize, the possibilities for singularity formation in this cosmological model are remarkably diverse. It's worth noting that some of the identified singularities may violate one or more energy conditions. These energy conditions encompass: ρ + p ≥ 0 "null energy condition" ρ ≥ 0 ρ + p ≥ 0 "weak energy condition" (27) ρ + 3p ≥ 0 ρ + p ≥ 0 "strong energy condition" With these considerations, we can succinctly summarize the findings for the cosmological model defined by the f (ρ) function as follows: • For A/B < 0, a type II singularity is inevitable, irrespective of the values of β.
• Regardless of the sign of A/B, the nature of singularities varies according to the values of β.
6. β < 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 (−∞) for B < 0 (B > 0).
So with this example (as was discussed in [31] ) 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 comes out to be qualitative differences when one considers singularities in Jordan and Einstein frames, something which was discussed in detail and discovered in [124,125]. It is also worth noting that when one considers viscous fluids, as in [126], then it may give rise to different types of singularities. The occurrence of singularities in an oscillating universe has been also been discussed, firstly in [127]. Singularities have also been considered in detail bounce cosmologies [128][129][130]. The realisation of all known 4 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 [131], while also in Teleparallel constructions like the one considered in [132].
A crucial point to that we should note here in passing with regards to all the singularities we have discussed so far is that the tidal forces which manifest for these singularities as the (infinite) impulse which reverses (or stops) the increase of separation of geodesics and the geodesics themselves can evolve further ; the universe can continue its evolution through a singularity. Moreover, it's 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 [133]: (with T denoting string tension, and τ, σ representing the string's worldsheet coordinates, η µν corresponding to the worldsheet metric, µ, ν = 0, 1, and g ab standing for the spacetime metric), the scenario involves the string interacting with a non-BB singularity [134]. The crux of the matter is that a measurable property of the string, its invariant size S(τ ) = 2πa(η(τ ))R(τ ) (assuming a circular assumption with radius R), reveals certain characteristics. Specifically, at a Big-Rip singularity, the string undergoes infinite stretching (S → ∞), resulting in its destruction. In contrast, at a type II singularity, the scale factor remains finite at the η-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. 2

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 firstly introduced in [139] and then expanded upon in later works [61,140]. The authors in [139] arrived at w-singularities by firstly choosing the scale factor ansatz as It contains seven arbitrary constants: A, B, C, D, γ, n, and t s . The last of the constants t s is the time when we expect the singularity. Having the scale factor (30), they imposed the following conditions The first of the conditions (31) is 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(0) = ∞, which is equivalent to taking γ < 0). One can see that after introducing (31), 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Ḧ(t s ) ̸ = 0 inḦ, and so it is not of the type IV singularity according to the classification of Ref. [31]. On the other hand, even though bothä(t s ) andȧ(t s ) vanish in the limit t → t s , the deceleration parameter blows-up to infinity, i.e., and consequently as one can find the EOS parameter to be related with the deceleration parameter as one finds that w(t s ) → ∞. Then, we face a very strange singularity. It has vanishing pressure and energy density, a constant scale factor, but the deceleration parameter and, in particular, a time-dependent barotropic index w(t) are singular. Another ansatz for the scale factor which can give w-singularities was proposed by Dabrowski and Marosek in [141], which has an exponential form. That ansatz is given by where a s has the units of length and is a constant while m and n are also constants 3 . 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 aṡ where the overdots now denote differentiation with respect to time. From this, one can see that for 1 < n < 2ȧ(0) → ∞ andȧ(t s ) = mas ts = const. , while a(t s ) = a s ,ä(0) → ∞ andä(t s ) → −∞ and we have sudden future singularities. Furthermore, it was shown in [141] that for the simplified case of the scale factor (20) with m = 0, one can get w-singularities for n > 0 and n ̸ = 1. Finally, yet another ansatz to get w-singularities was provided in [61] and is of a power series form given by 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 amount of works which have considered how these singularities can occur in non-standard cosmologies and how they can be avoided too. But in passing, a discussion is in order over the cosmological significance of w-singularities. While Type I-Type IV singularities deal with more direct cosmological parameters like the scale factor, 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 [142] discussed how w-singularities can occur in interacting dark energy models(while the background cosmology in this case was still general relativistic and the continuity equation had its usual form), while in [143] it was 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 [103,[144][145][146][147]. Hence while type V singularities deal primarily with a more indirect cosmological parameter, it by no means diminishes its cosmological importance and it does appear in a variety of cosmological expansion scenarios.

Singularity removal/avoidance methods
With huge influx of 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, look at [53] for a detailed account of avoiding singularities in both Jordan and Einstein frames). What we would like to do here is to discuss some of the prominent works in this regard, focusing on the use of quantum effects and modified gravity effects to deal with singularities.

Conformal anomaly effects near singularities
The effect of quantum backreaction of conformal matter around Type I, Type II and Type III singularities were taken into consideration in the works of Nojiri and Odintsov [29,31,148]. 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 a bit of context about what conformal anomalies are and how they are usually perceived in high energy physics. It is fair to assume that there are 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 behaviour of these (massless) matter fields-scalars, the Dirac spinors, and vectors in curved space-time-is conformal 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. But renormalization procedures can lead to the trace of an anomalous energy momentum tensor, which is the so-called quantum anomaly or the conformal anomaly(we would recommend the reader [149][150][151][152] for more details on conformal anomaly effects). The conformal anomaly we have described be considered to have the following form [31] T 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 b and b ′ on the other hand are given by with N scalar, N 1/2 spinor , N 1 vector fields , N 2 (= 0 or 1 ) gravitons and N HD being higher derivative conformal scalars. For usual matter, b > 0 and b ′ < 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 ρ A and pressure p A . The total energy density is ρ tot = ρ + ρ A . The conformal anomaly, also known as the trace anomaly, can be given by the trace of the fluid stress-energy tensor The conformal anomaly corrected pressure and energy densities still obey the continuity equation (2.10) and using that, we can write The conformal anomaly corrected pressure and energy densities still obey the continuity equation and using that, we can write [31] T One can then express ρ A as an integral in terms of T A as Furthermore T A can be expressed in terms of the Hubble parameter as And using this, one can have an expression for ρ A taking into account conformal anomaly effects near the singularity While the quantum corrected Friedmann equation 4 is 3 Since the curvature is expected to be large near the time of the singularity, one can be warranted to think that (3/κ 2 )H 2 ≪ |ρ A |. Then ρ ∼ −ρ A from (48), which giveṡ (50) Finally then, the continuity equationρ + 3H (ρ + p) = 0 for p = −ρ − f (ρ), gives Now we can appreciate the implications of these effects on both strong and weak singularities, where firstly we consider the big rip. The first attempt to address the issue of big rip with conformal anomalies was done in [153,154] For this, we consider the model given by with 1/2 < β < 1 when ρ is large and in this case there exists the Big Rip singularity, as we had discussed previously too. in Sec. IV. We note that the classical evolution is characterized by ρ ∝ (t s − t) 2 1−2β and H ∝ (t s − t) 1 1−2β , both of which exhibit divergence for β > 1/2. When quantum corrections are taken into account, it is natural to assume that near the singularity ρ behaves as 4 Note that to maintain consistency with the notation used in [31] we are considering the Friedmann equation to be of the form As ρ may diverge at t = t s , we consider negative values ofγ. Sinceγ (1 − β) < 0 in this case, we might expect that (50) would give the following approximate relation around t = t s : The term on the r.h.s. grows as H 4 ∝ (t s − t) −4+4γ(1−β) , but this does not give a consistent result, since ρ becomes negative for b ′ < 0. This tells that our assumptions should be wrong and ρ does not become infinite. If ρ has an extremum, (51) tells that H vanishes there sinceρ = 0. Furthermore, the authors in [31] 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 (52) but now for the range β > 1, in which case we see that a type III singularity develops with ρ ∝ (t s − t) 2 1−2β . Again, we consider that near the singularity ρ behaves as (53). Using (51), one finds Since we are considering the case β > 1 andγ < 0, we have thatγ (1 − β) > 0. By picking up most singular term in the r.h.s of (50), it followṡ Then substituting (53) and (55) for (56), we obtaiñ γ = 4 1 − 2β (57) This means that ρ and H evolve as around t = t s . Numerically solving the background equations shows that in the presence of quantum corrections one has H ∝ (t s − t) 1/3 around t = t s , which means that H approaches zero. Meanwhile in the absence of quantum where a 0 is a constant. Comparing the classical case [γ = 2/(1 − 2β)] with the quantum corrected one [γ = 4/(1 − 2β)], we find that the power of (t s − t) 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 ρ is stronger. Thus quantum effects moderate the classical singularity.
But conformal anomaly effects may not always be of a huge help in order to alleviate singularities, take for example the case of an asymptotically safe cosmology which was considered in [155]. The capacity to build gravitational RG flow approximations outside of perturbation theory is necessary for conceptually testing asymptotic safety. a very strong framework for doing these calculations is the functional renormalization group equation (FRGE) for the gravitational effective average action Γ k The construction of the FRGE uses the background field formalism, where the metric g µν is split into a fixed background g µν and fluctuations h µν ( see [156] for a more detailed on asmyptotitcally safe cosmologies ). The authors in [155] considered the simplest approximation of the gravitational RG flow, which could be obtained from projecting the FRGE onto the Einstein-Hilbert action approximating Γ k by [156] gauge-fixing and ghost terms (61) where R, Λ k and G k are the Ricci Scalar, the running cosmological constant and the running Newton's gravitational constant. The scale-dependence of these couplings can be written in terms of their dimensionless counterparts as where g * = 0.707 and λ * = 0.193 . Considering a background FLRW metric and a perfect fluid for the stress energy tensor T ν µ = diag[−ρ, p, p, p] , one can get the Friedmann equation and the continuity equation in this scenario to be Where the continuity equation comes about from the Bianchi identity which is satisfied by Einstein's equations D µ [λ(t)g µν − 8πG(t)T µν ] = 0, which has the usual meaning that the divergence D µ of the Einstein tensor vanishes . The extra terms of the right hand side in (65) 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 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 as But obtaining a corresponding integral for ρ A for the equation (66) is not possible in the same way . Hence its not feasible to address a possible removal of Type I-Type III singularities using conformal anomaly effects in this asymptotically safe cosmology.

Varying constants approach
Cosmologies with varying physical constants, like speed of light or the gravitational constant [157] have been shown to regularize cosmological singularities in certain scenarios [141,[158][159][160]. Here we shall discuss briefly about 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 the varying speed of light c(t) (VSL) and varying gravitational constant G(t) (VG) as presented by Barrow in their work [157], one can deduce the following expressions for mass density ϱ(t) and pressure p(t): These equations highlight the influence of varying c and G on mass density and pressure. For instance, ifȧ approaches infinity while G(t) increases more rapidly thanȧ, the singularity in p(t) can be eliminated. In the case of flat models, a direct relationship between pressure p and mass/energy density ρ/ε can be established, albeit with a time-dependent equation of state parameter, expressed as: Here, the parameter w(t) is defined as w(t) = 1 3 [2q(t) − 1], with q(t) = −äa/ȧ 2 being the dimensionless deceleration parameter. Notably, the variation of the speed of light c(t) brings about a key distinction between mass density ϱ and energy density ε = ϱc 2 , impacting the Einstein mass-energy relationship E = mc 2 , which is transformed here into the mass density-pressure formula p = ϱ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: The constants t s , a s , m, n are determined accordingly [141]. This approach illustrates how the varying constant concept aids in regularizing singularities. By inspecting Equations (67) and (68), it becomes evident that a time-dependent gravitational constant variation of the form G(t) ∝ 1 t 2 eliminates a Type 0 Big Bang singularity in Friedmann cosmology, addressing both p and ϱ singularities. In Dirac's scenario [161], where G(t) ∝ 1/t, only the ϱ singularity is removed. Moreover, the time dependence of G = 1/t 2 is less constrained by geophysical limitations on Earth's temperature [162].
Another proposal suggests that if the scale factor (70) doesn't tend to zero as t → 0, it could be rescaled by a "regularizing" factor a rg = (1 + 1/t m ) (m ≥ 0), resulting in: 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(t). One common assumption regarding the speed of light's variation is that it follows the evolution of the scale factor [157]: With c 0 and s as constants, the field equations (67) and (68) can be expressed as: In the presence of the time dependence of c(t) as given by (72), and for the choice of a(t) = 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 ̸ = 0.

Modified gravity effects/Quantum gravitational cosmologies
In recent times, there has been a wide interest in dark energy models based in exotic non general relativistic regimes particularly because such theories display properties which are not evident in conventional cosmological models. For example,a lot of 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 [15][16][17][18][19]. There has been substantial work on new dark energy models based on thermodynamic modifications like modified area-entropy relations [163][164][165][166][167][168] or even more exotic possibilities like generalized uncertainty principles [169][170][171] or non-canonical approaches like DBI etc. as well [172][173][174][175][176][177][178].This vast dark energy literature has prompted the study of cosmological singularities in a wide range of cosmological backgrounds as well, as there have been multiple works which have discussed Type I-IV singularities in various cosmologies [30,31,33,34,105,155,[179][180][181][182][183][184][185][186][187][188][189][190]. In this vast array of works, one has seen quite a few examples in which cosmologies affected by the effects of these modified gravity theories or quantum gravitational paradigms ( like the Braneworld or LQC for example ) 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 discussed ways to remove or moderate this till now.
We would like to consider the treatment in [190] for our example here. We would like to again consider a model with an inhomogeneous EOS of the form p = −ρ − f (ρ). It was shown in [141] that for the simplified case of the scale factor (34) with m = 0, one can get w-singularities for n > 0 and n ̸ = 1. The scale factor for the case m = 0 takes the form and we will be using this form of the scale factor for this example. The modified gravity theory we would be interested in is an f(R) gravity model with the action [191] S = m 2 where α is a constant which has the units of mass, L m is the Lagrangian density for matter and m p is the reduced Planck's constant. The field equation for this action is The Friedmann equation in this case can take the form where ρ is the total energy density. This F (R) gravity regime was used to explain late time cosmic acceleration as an alternative to dark energy in [191].
The use of f(R) gravity regimes for avoiding cosmological singularities by adding an R 2 term was considered in detail in [192], with the same scenario later being extended in [32] and [181]. Moreover, based on properties of R 2 term, nonsingular modified gravity were proposed in [193]. The action (77) prompts one towards the notion that very tiny corrections to the usual Einstein Hilbert 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 [194], the authors in [191] proposed a purely gravitational paradigm through (76) to explain both the early and late time accelerations of the universe. Now we consider f (ρ) = ρ α and firstly see the status quo of w-singularities for such a model in the standard cosmology given by ( written here in natural units just for simplicity) H 2 = ρ 3 . We can write the w-parameter for this cosmology using as From this we can make the following observations : • For n = 1, no w-singularities occur as is the case in the usual scenario with conventional equation of state.
• For α < 0 , w-singularities occur for all values of positive values of n besides unity but w-singularities do not occur for any negative values of n • For α > 0 we see a very interesting behaviour. 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.
So 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 (78) case is given by For the w-parameter as expressed above, we have the following observations : • For n = 1, contrary to the other cases we have considered till now, one can have a w-singularity but that is possible only in the extreme case that α → ∞ which is not pretty realistic to expect but in principle singularities can appear in this case.
• The most interesting thing that comes out when one considers this scenario is that w-singularities do not occur for any value of n and α ! For both positive and negative values of α and n, the w-parameter remains regular and does not diverges.
And so we see that just by incorporating the effects of a modified gravity theory, in this case a particular form of f (R) gravity, one can alleviate singularities too. Furthermore, f (R) gravity theories which been of great use in alleviating various other singularities too but we have discussed about other singularities quite extensively before and so it seems appropriate to discuss this example to illustrate how type V singularities could be moderated too.

Dynamical systems approach and the Goriely-Hyde method
While it is seems pretty natural to study singularities and their avoidance methods in various cosmological settings like we have discussed so far,often it is very difficult to classify and study the cosmological singularities which may occur in extremely non-conventional cosmologies which are motivated by quantum gravitational/phenomenological considerations (for example, see the classification of singularities in asymptotically safe cosmology [155] ) and often it may not even be possible to do so in an orthodox fashion. Hence it becomes essential to look for non-conventional ways to find out 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 [195]. As cosmology presents a multitude of captivating dynamical systems [196], the investigation of singularity structure in such systems has gained considerable attention, with the Goriely-Hyde method proving particularly useful for cosmological explorations [197][198][199][200][201][202]. This method has previously been applied to study finite and non-finite time singularities in certain classes of quintessence models as well [33,182,203].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:ẋ where i = 1, 2, ..., 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" [195]. Each dominant part constitutes a mathematically consistent truncation of the system, denoted asf i . The system can then be written as: • Without loss of generality, the variables x i near the singularity can be expressed as: where τ = 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, which form the vector p = (p 1 , p 2 , ..., p n ). Similarly, we calculate the values of a i to form the vector ⃗ a = (a 1 , a 2 , ..., a n ). It is important to note that if ⃗ a contains only real entries, it corresponds to finite-time singularities. Conversely, if ⃗ 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: 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 , ..., r n ), with r 2 , r 3 , ... > 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.
• Without loss of generality, the variables x i near the singularity can be expressed as: where τ = 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, which form the vector p = (p 1 , p 2 , ..., p n ). Similarly, we calculate the values of a i to form the vector ⃗ a = (a 1 , a 2 , ..., a n ). It is important to note that if ⃗ a contains only real entries, it corresponds to finite-time singularities. Conversely, if ⃗ 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: 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 , ..., r n ), with r 2 , r 3 , ... > 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 motivated ansatz' 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 [203]: Here, f 1 (t) and f 2 (t) are assumed to be nonzero regular functions at the time of the singularity, and similar conditions apply to their derivatives up to the second-order. Additionally, α 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 α < −1, a big rip singularity occurs.
Another ansatz useful for classifying singularities was introduced in [39] whereby the scale factor was written as (88) where 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 α one can have the following singularities • For α < 0, a type I singularity occurs • For 0 < α < 1, a type III singularity develops • For a < α < 2, a type II singularity occurs • For α > 2, a type IV singularity occurs Again, it is not mandatory that the scale factor in (88) will necessarily be a solution to the field equations but we would like to consider this and (87) in order to get a well-motivated feel for the type of cosmological singularities we can deal with in the various models we have discussed so far. As an example for this method, let's consider singularities in an RS-II braneworld cosmology where dark energy can be described by a scalar field paradigm, where we shall follow the treatment of [33]. The action for inclusive of both the scalar and the background fluid term can be written as where R (5) , g µν and Λ (5) are the bulk Ricci Scalar, metric and the cosmological constant respectively with σ being the brane tension on the 3-brane, g µν being the 3-brane metric and µ(ϕ) being a scalar coupling function. Note that here we are working in Planck units with (m (5) p ) 2 = 1 with m (5) p being the 5-dimensional Planck mass. Assuming that the brane metric has the usual FLRW form, we get the Friedmann equation to be [204] where ρ = ρ ϕ + ρ 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 And the equation motion of the scalar is given by Finally, using the following variables introduced in [205] x Choosing the background fluid to be of the form of pressurelees dark matter, in a way that w B = 0 , we get the dynamical system for this model to be where the primes denote differentiation with respect to the e-folding number N and λ = V ′ V . We can finally start with the analysis as we have proper autonomous dynamical system, with the first truncation that we consider beinĝ where k = 3 2µ . Using the ansatz of the Goriely-Hyde method, we get p = (−1, −1, −1) and using these, we get 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 We then finally find the eigenvalues of the matrix, which are given by Hence the singularities in this case will only be local singularities which will only form for a limited set of initial conditions. In [33], it was worked out that there are two more possible truncations, with the balances and corresponding eigenvalues in them being with eigenvalues being And another truncation with the balance with We see from (101) and (102) that the truncations for which they belong to seem to still tell the story that only non-finite time local singularities could be 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φ and V (ϕ) in terms of the Hubble parameter. For simplicity, we will consider that the coupling constant µ = 1 andρ B = 0. Making these considerations, we can write −2Ḣ =φ 2 1 + ρ σ One can then writė Furthermore, one can now write V (ϕ) in terms of the dark energy equation of state 5 as Then we can write the potential as where k = 2w 1−w and b = σ(1 + ρ 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 σ, w and ρ B ) and so one can use this form of V in to findφ in terms of the Hubble parameter as well. It is necessary to express these quantities in terms of H(t) as now we can find out which type of singularities are possible in this scenario, in the view of the fact that x,y and z as described in have to remain regular. By studying the expressions for these variables, one can make out that Type I, Type III and Type IV singularities are allowed in this scenario while Type II is 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.

Future outlooks and conclusions
In this brief review paper, we have discussed (almost) all the prominent developments in the field of cosmological singularities which have taken place in the recent 25 years or so. We firstly provided a detailed outlook on what space-time singularities are and discussed the various nuances regarding them, like various strength criterion etc. After that, we discussed in detail about 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 found in. We then saw how one can moderate or even remove these singularities using various techniques having 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 could safely make is that such singularities provide a revealing arena on the interface of cosmology and quantum gravitational theories. The scales at which such events could take place lie in the horizons for testing quantum gravity ideas and with the constant increase in the precision of various observational setups, one would not be wrong to think that investigating such singularities in detail can possibly shed light on problems in both cosmology and quantum gravity.

Acknowledgments
The author would like to thank Sergei Odintsov for the invitation to write this review article and for the numerous discussions with him on various aspects related to singularities. I would also like to thank Maxim Khlopov, Pankaj Joshi, Robert Scherrer, Alexander Timoshkin, Vasilis Oikonomou, Jackson Levi Said and Rafael Nunes for various discussions on singularities. I would also like to thank Parth Bhambhaniya for discussions on some particular aspects of singularities. Finally, I would also like to thank Sunny Vagnozzi for discussions on cosmology and dark energy in general which have been very helpful for this work.

Appendix A
The use of power series expansions for the scale factor and related quantities in cosmology has gathered significant pace in recent times(for a detailed overview, see [112]) and especially in the context of studies of cosmological singularities. Hence it is fitting to discuss such expansions in a bit of detail here. Generalized Frobenius series find frequent application in the expansion of solutions to differential equations around their singular points. With this characteristic in mind, we will assume that in the vicinity of the key point, the scale factor exhibits an expansion based on a generalized power series. This concept extends the familiar notions of Taylor series, meromorphic Laurent series, Frobenius series, and Liapunov expansions, as referenced in [206]. Furthermore, this generalized power series is more encompassing than the one employed in [207].
In the current context, if the scale factor a(t) can be expressed using such a generalized power series, then the Friedmann equations dictate that both ρ(t) and p(t) can likewise be represented using such power series. Employing formal series reversion, it follows that the equation of state ρ(p), and consequently the function ρ(a), exhibit these generalized power series expansions. Conversely, when ρ(a) is described by such a generalized power series, the first Friedmann equation indicates thatȧ(t) possesses a power series of similar nature, which, upon integration, implies that a(t) itself is characterized by such a power series.
Similarly, if the equation of state p(ρ) can be expressed as a generalized power series, then integrating the conservation equation leads to the expression: This equation also adopts a generalized power series representation. The potential value of expanding the conventional notion of a Frobenius series becomes evident through the analysis presented in [207].
It is important to clarify the types of entities that fall outside this category of generalized power series. First, essential singularities, effectively infinite-order poles that emerge, for instance, in functions like exp(−1/x) near x = 0, lie beyond this classification. Secondly, certain variations on the concept of Puiseux series, specifically those containing terms like (ln x) n , (ln ln x) n , (ln ln ln x) n , and so forth, also exist beyond this classification. However, there is currently no awareness of any scenarios where these exceptional cases become pertinent in a physical context.
It has been shown to be reasonable to assume that in the vicinity of some cosmological singularity, happening at some time t 0 , the scale factor has a (possibly one-sided) generalized power series expansion of the form where the indicial exponents η i are generically real (but are often non-integer) and without loss of generality are ordered in such a way that they satisfy Finally we can also without loss of generality set c 0 > 0. (112) There are no a priori constraints on the signs of the other c i , though by definition c i ̸ = 0.
From a physical point of view, this definition is really generic and can be applied to any type of cosmological milestone. This generalized power series expansion is sufficient to encompass all the models we are aware of in the literature, and as a matter of fact, the indicial exponents η i will be used to classify the type of cosmological singularity we are dealing with. For many of the calculations in this chapter, the first term in the expansion is dominant, but even for the most subtle of the explicit calculations below it will be sufficient to keep only the first three terms of the expansion: a(t) = c 0 |t − t 0 | η0 + c 1 |t − t 0 | η1 + c 2 |t − t 0 | η2 . . . ; η 0 < η 1 < η 2 ; c 0 > 0. (113) The lowest few of the indicial exponents are sufficient to determine the relationship between these cosmological milestones, the curvature singularities and even the energy. Note also that this expansion fails if the cosmological milestone is pushed into the infinite past or infinite future. Using such an expansion, one can encounter quite a few cosmological singularities and we shall list some conditions under which some prominent singularities 6 can be found as follows : • Big bang (Type 0): If a big bang occurs at time t 0 7 , we define the behavior with indicial exponents (0 < η 0 < η 1 . . . ) when the scale factor has a generalized power series near the singularity, given by: The series is carefully constructed such that a(t 0 ) = 0.
• Big rip (Type 1): If a big rip occurs at time t 0 , the indicial exponents of the rip (η 0 < η 1 . . . ) are defined when the scale factor has a generalized power series near the rip: where η 0 < 0 and c 0 > 0. The series is constructed to satisfy a(t 0 ) = ∞. The only difference from the big bang case is the sign of the exponent η 0 .
Here, c 0 > 0 and η 1 is a non-integer. The condition a(t 0 ) = c 0 ensures finiteness, and a sufficient number of differentiations yields: a (n) (t → t 0 ) ∼ c 0 η 1 (η 1 −1)(η 1 −2) . . . (η 1 −n+1) |t−t 0 | η1−n → ∞. (117) The toy model by Barrow [56] can be expressed as: where t b is the time of the big bang. This model fits into the general classification when expanded around the sudden singularity time, t 0 , and into the classification of big bang singularities when expanded around the big bang time, t b .

Appendix B
Over the years, several alternatives to the Big rip have been found and the first one that we shall discuss in this regard is the Little rip (LR). It is characterized by a growing energy density ρ over time, but this increase follows an asymptotic pattern, necessitating an infinite amount of time to approach the singularity. This situation corresponds to an equation of state parameter w that falls below -1; however, it approaches -1 as time progresses towards infinity. The energy density's growth is gradual, preventing the emergence of the Big Rip singularity. The LR models depict transitional behaviors between a asymptotic de Sitter expansion and a BR evolution. In their work [108], the authors presented an elegant method for comprehending the implications of the little rip, distinguishing it from the Big Rip, which we will explore in the following. During the universe's expansion, the relative acceleration between two points separated by a comoving distance l can be expressed as lä/a, where a signifies the scale factor. If an observer is situated at a comoving distance l from a mass m, they will detect an inertial force acting on the mass as follows: Let's assume that the two particles are bound by a constant force F 0 . When the positive value of F iner surpasses F 0 , the two particles become unbound. This scenario corresponds to the phenomenon known as the "rip," which emerges due to the accelerating expansion. Equation (119) demonstrates that a rip always occurs when either H orḢ diverges (assumingḢ > 0). The divergence of H leads to a "big rip", while if H remains finite butḢ diverges withḢ > 0, it results in a Type II or "sudden future" singularity [31,56,197], which also causes a rip.
Nonetheless, as pointed out in [208], it's feasible for H and, consequently, F iner , to grow boundlessly without inducing a future singularity at a finite time. This phenomenon is referred to as the little rip. Both the big rip and the little rip share the characteristic of F iner → ∞; the distinction lies in the fact that for the big rip, F iner → ∞ occurs at a finite time, whereas for the little rip, it occurs as t → ∞. Two possible ansatz/models which have been shown to lead to little rip behaviour [108] are given by the following forms of the Hubble parameter where H 0 and λ are positive constants while another viable model which is similar to this is given by where H 0 , λ and C are positive constants as well. Another interesting possibility for the evolution of the universe is the so-called Pseudo-Rip [109], where the Hubble parameter, although increasing, tends to a "cosmological constant" in the remote future. That means, H(t) → H ∞ < ∞, t → +∞, where H ∞ is a constant. A possible model for this is given by the Hubble ansatz given as where H 0 , H 1 and λ are positive constants with H 0 > H 1 . Yet another possible alternative for the rip is a model in which the dark energy density ρ monotonically increases (w < −1) in the first stage, and thereafter monotonically decreases (w > −1), known as the "Quasi rip" [110]. At first, it thus tends to disintegrate bound structures in the universe, but then in the second stage the disintegration becomes reversed, implying that already disintegrated structures have the possibility to be recombined again. As an example model for this, we consider the energy density of dark energy to be a function of the scale factor and consider it's anastz to be ρ(a) = ρ 0 a α−β ln a where a is the scale factor, α and β are constants with ρ 0 being the energy density at some past time t 0 . Yet another possibility is the Little sibling of the big rip [209] wherein the Hubble rate and the scale factor blow up but the derivatives of the Hubble rate does not. This takes place at an infinite cosmic time with the scalar curvature blowing up too. An example model for this also involves taking the energy density of dark energy as a function of the scale factor, given by [209] ρ(a) = Λ + A ln a a 0 (124) Table 1 summarizes all these various scenarios as discussed till now. Further- Quasi Rip Dark energy density ρ first increases (w < −1) and then decreases (w > −1), implying disintegration and recombination of structures.
ρ(a) = ρ 0 a α−β ln a Little Sibling of the Big Rip Hubble rate and scale factor diverge, but derivatives of Hubble rate do not, with scalar curvature divergence.