Singularity-free and cosmologically viable Born-Infeld gravity with scalar matter

The early Cosmology driven by a single scalar field, both massless and massive, in the context of Eddington-inspired Born-Infeld gravity, is explored. We show the existence of nonsingular solutions of bouncing and loitering type (depending on the sign of the gravitational theory's parameter, $\epsilon$) replacing the Big Bang singularity, and discuss their properties. In addition, in the massive case we find some new features of the cosmological evolution depending on the value of the mass parameter, including asymmetries in the expansion/contraction phases, or a continuous transition between a contracting phase to an expanding one via an intermediate loitering phase. We also provide a combined analysis of cosmic chronometers, standard candles, BAO, and CMB data to constrain the model, finding that for roughly $\vert \epsilon \vert \lesssim 5\cdot 10^{-8} \text{ m}^2$ the model is compatible with the latest observations while successfully removing the Big Bang singularity. This bound is several orders of magnitude stronger than the most stringent constraints currently available in the literature


Introduction
The standard concordance cosmological ΛCDM model, framed within Einstein's General Theory of Relativity (GR), including an early phase of inflationary expansion, a cold dark matter component, and a tiny cosmological constant driving the accelerated late-time expansion of the Universe, has successfully met all observations [1,2]. Within this model, scalar fields have found new and imaginative applications. For instance, inflationary models in the early Universe involve from one to many scalar fields [3][4][5][6][7][8][9][10][11][12]. In the slow-roll approximation the exact form of the scalar field potential is unknown since many different potentials have been studied and confronted to observations. On the other hand, a way to parameterize dark energy is by using a scalar field, the so-called quintessence model (for canonical scalar fields [13,14]) or its generalizations to K-essence models (when a non-canonical scalar Lagrangian is considered [15][16][17][18]), in such a way that the cosmological constant gets replaced by a dark energy fluid with a nearly constant density today [19][20][21][22][23]. Dark matter can be also parameterized in terms of weakly-interacting massive particles, which can be scalar particles still undiscovered at colliders and other dark matter detection experiments. Models for dark matter can also be based on other kinds of scalar fields, for instance, via fuzzy dark matter [24], or by using a Lagrange multiplier that changes the behaviour of the kinetic term [25][26][27][28]. Scalar field models may also be the result of complex effective interactions of other fundamental fields in equilibrium, such as in Bose-Einstein condensates, thus allowing for an even broader range of phenomenological justifications. Therefore, the consideration of scalar fields within cosmological models is justified from the point of view of building a fully consistent history of the cosmological evolution compatible with observations. S EiBI = 1 1 This proposal can actually be framed within the tradition of considering square-root action, such as in the DBI one, see e.g. [36][37][38][39].
where κ 2 ≡ 8πG/c 4 is Newton's constant, is EiBI parameter, λ is a dimensionless constant, g is the determinant of the space-time metric g µν and q the determinant of an auxiliary metric defined as: where the (symmetric part of the) Ricci tensor R µν (Γ) ≡ R ρ µρν (Γ) is a function solely of the (torsionless) affine connection Γ ≡ Γ λ µν , assumed to be a priori independent of the space-time metric g µν (metric-affine or Palatini formalism). This symmetrizing requirement ensures that the theory is invariant under projective transformations, which avoids the presence of ghost-like instabilities [52,53]. As for the matter sector, S m = d 4 x √ −gL m (g µν , ψ m ), it is assumed to be minimally coupled to the space-time metric g µν (see e.g. [54] for a definition of minimal coupling in metric-affine theories), with ψ m denoting collectively the matter fields. It is worth pointing out that EiBI gravity recovers the GR dynamics and its solutions in the |R µν | −1 limit [51] (from now on, vertical bars indicate a determinant), inheriting an effective cosmological constant Λ e f f = λ−1 κ 2 . This fact allows EiBI gravity to naturally pass weak-field limit tests, such as those based on solar system observations.
EiBI gravity is nowadays a very well known theory in the community thanks to its many applications (for a detailed account of its properties we refer the reader to the review [51]). It is actually a member of the so-called Ricci-based family of gravitational theories [55], all of which admit an Einstein-like representation of their field equations given by where L G is the gravitational Lagrangian, T µν = 2 √ −g δS m δg µν is the stress-energy tensor of the matter fields, T its trace, and G µ ν (q) is the Einstein tensor of the auxiliary metric q µν , such that Γ is Levi-Civita of it, that is This q µν metric is related to the space-time one g µν via the relation where the deformation matrixΩ depends on-shell on the matter fields (and possibly the space-time metric g µν as well). For EiBI gravity this matrix is implicitly determined via the equation while the EiBI Lagrangian in (1) can also be expressed in terms of this matrix as Let us point out that all terms on the right-hand side of the equations (3) are functions of the matter fields and the metric g µν , thus representing a system of second-order field equations with new couplings engendered by the matter fields. In vacuum, T µ ν = 0, one recovers the GR solutions 2 , which ensures the propagation of the two polarizations of the gravitational field travelling at the speed of light [57].

EiBI cosmology with scalar fields
As the matter sector of our model let us consider a single (real) scalar field described by the action with kinetic term X ≡ g µν ∂ µ φ ∂ ν φ and scalar potential V(φ). The variation of this action with respect to the scalar field leads to the field equations (where V φ ≡ dV/dφ) and to the associated stress-energy tensor We consider next a spatially flat Friedman-Lemaitre-Robertson-Walker (FLRW) space-time 3 , with line element: where a(t) is the expansion factor and d x 2 = δ ij dx i dx j , with i, j = 1 . . . 3 for the spatial part. As from now on we shall assume φ = φ(t), the stress-energy tensor (10) reads: where a dot denotes a derivative with respect to t, as usual. From Eq.(6) we can conclude that the deformation matrix needs to have a similar algebraic (diagonal) structure as that of the stress-energy tensor of the scalar field, namely As a remark, we would like to note that it has been shown that in generic RBGs, due to the non-linearities of the field equations, the deformation matrix admits other (typically pathological) solutions besided the one that boils down to GR in vacuum. These solutions can lead to anisotropic deformation matrices even if the stress-energy tensor is isotropic. Remarkably, it was shown that for EiBI gravity, no such anisotropic solutions exist in presence of an isotropic stress-energy tensor, see [56].
where the identity matrix I 3×3 represents the spatial sector. Now, plugging this ansatz and Eq.(12) into Eq. (6) gives the components of this matrix as: where we have used the shorthand notations Therefore, we are ready to cast the field equations (3) for this problem as where the determinant of the deformation matrix reads |Ω| 1 2 = Ω 1/2 + Ω 3/2 − . Now, to compute the left-hand side of the field equations (18) one can write another line element for the auxiliary metric also of the FLRW form, that is where in the second line we have used the fundamental relation (5) together with the ansatz (13). These gravitational field equations must be supplemented with the scalar field equations (9) which, in the FLRW background (11), readφ Now, using the fact that in the q µν geometry (19) we have the well known formulas with the relations it follows that the combination can be written as and a little algebra leads to 1 a d dt aΩ which allows us to find an expression for H ≡ȧ/a as The square of this quantity is the generalized FLRW equation for EiBI gravity coupled to a scalar field with a potential V(φ). The ± signs in front of the square root yield expanding (+) or contracting (-) universes and must be chosen on physical grounds. In the limit → 0 the above expression yields which recovers the FLRW dynamics of GR coupled to a scalar field. Note that, from Eq.(12), for a scalar field the energy density ρ φ and pressure P φ are related to the field variables by which allows to write the quantities inside the functions Ω ± of Eqs. (14) and (15) in terms of ρ φ and P φ as

Massless scalar fields
To explicitly solve the gravitational and scalar field equations (18) and (9) we need to specify a form of the scalar potential V(φ). For the sake of simplicity, and to make contact with similar settings in black hole scenarios [59], let us consider first the case of a free scalar field, V = 0, for which the scalar field equation (20) has a first integralφ withφ 0 an integration constant. This equation is formally the same as in the GR case, though the scale factor implicitly contains the -corrections via the resolution of Eq. (27), smoothly reducing to their GR values in the limit → 0, as follows from Eq.(28). Since in this case P φ = ρ φ =φ 2 /2, the components of the deformation matrix (13) assume the relatively simple form where we have defined the critical density ρ = 1/(κ 2 | |) and s = ±1 denotes the sign of . The corresponding Hubble function takes the form As mentioned before, at low densities as compared to the scale ρ , this equation recovers the GR expression, which can be solved leading to the approximate solution a(t) = H 0 t 1/3 . The Hubble factor here can be conveniently set to H 2 0 = 3κ 2φ2 0 /2 to adjust the integration constantφ 0 in order to reproduce the standard cosmological evolution of GR coupled to a massless scalar field. At higher densities, however, the dynamics strongly departs from that of GR and it is evident that when the energy density of the scalar field approaches its maximum value, ρ φ = λρ , the Hubble function vanishes. This corresponds to a minimum value of the expansion factor for both branches of solutions s = ±1, but there are two different mechanisms by which the initial Big Bang singularity is avoided, which we discuss next. For simplicity, from now on we focus on asymptotically flat configurations, λ = 1.
For the s = −1 branch, at the critical density ρ φ = ρ , the Hubble factor in Eq. (35) and its derivative behave as This indicates that at the maximum achievable density, ρ φ = ρ , the Hubble factor vanishes while its derivative goes to infinity. This implies that the universe contracts to a minimum (maximum) value of the expansion factor (energy density), before re-expanding. This the typical behaviour expected in standard bouncing solutions with a transition from a contracting phase to an expanding one, where the universe could even undergo a sequence of cyclic cosmologies with both phases. In Fig. 1 we plot the form of H(ρ) (blue curves) according to Eq. (35) to show that while at late-times (i.e. low densities, ρ φ ρ ) these solutions converge to those of GR (in agreement to Eq.(28)), at high densities (ρ φ /ρ → 1) they depart form the standard Big Bang singularity of GR.
For the s = +1 branch we find that the evolution of the Hubble factor as we approach the maximum density ρ φ = ρ is given instead by which means that the Hubble factor vanishes there as well, but instead of being divergent its derivative takes a finite value (which is actually zero). This implies that the expansion factor reaches a fixed value a(t) = a m , corresponding to an asymptotically Minkowski past, and starts expanding as we move forward in time after some reference time t = t p . In Fig. 1 we depict this behaviour (red curves) of these so-called loitering solutions, where the qualitative differences with the bouncing solutions are manifest. Again, at late times (low densities) the standard GR evolution is recovered. Therefore we see that a free massless scalar field coupled to EiBI gravity is able to yield nonsingular evolutions in both s = ±1 branches according to these two mechanisms, something not possible within GR.

Massive scalar field
Let us now consider a massive scalar field with a potential of the form where µ is a constant. This choice for the potential, besides being the simplest form, corresponds to the natural mass associated to a canonical scalar field. In this case, the Hubble function (27) develops an explicit dependence on both ρ φ and P φ [to lighten the notation from now on we shall drop the label φ in these quantities] via Eqs. (31) and (32). Since the resulting expression does not provide any useful insight, we will not write it explicitly here. Instead, for the sake of comparison with the massless case of the previous section, we find it more illustrative to provide a parametric representation of H(t) versus ρ(t). These functions are obtained by numerically integrating the second-order equation forä that follows from  (22) together with its corresponding right-hand side. The result of the integration is used to construct the quantityȧ/a, which is then compared with the formula (27) to check the consistency of the numerical integration. Again we split our discussion of the corresponding results into the s = ±1 cases.    Fig. 4. Note that the blue and magenta curves will bounce at some point and begin a growing oscillatory trajectory similar to that of the red dashed curve, being this a generic behavior of all these solutions.
In Figs. 2 and 3, we show the function H(ρ(t)) from Eq. (27) with the help of (14) and (15) particularized to the potential above, for solutions corresponding to several values of the reduced mass σ ≡ µ 2 /ρ in the case s = −1. The corresponding expansion factors a(t) [obtained from integrating the previous function H] appear in Fig. 4, while in Fig. 5 we show the functionȧ(t) of those same examples and in Fig. 6 their Hubble functions. As one can see from all these plots, the case s = −1 still represents bouncing solutions for all the mass range explored, which goes from σ ≈ 10 −8 up to σ ≈ 0.4, with little variations with respect to the massless scenario for masses as high as σ ≈ 10 −3 (see Fig. 2). For higher masses, the egg-shaped Hubble function develops a fish-like structure, with asymmetric fins. This asymmetry is also manifest in the contracting branch of the expansion factor (see Fig. 4), which becomes increasingly asymmetric as the mass parameter grows from zero. A similar behaviour can be observed in the density profiles of the scalar field depicted in Fig. 7.  Focusing now on the case with s = +1, which represents loitering solutions in the massless case, Fig.  8 shows that the function H(ρ(t)) is now qualitatively different (compare to Fig. 1). Indeed, when σ = 0, H(t) never crosses the horizontal axis as one approaches the maximum allowed energy density. However, for σ = 0 one either has an expanding branch which starts from an almost constant a(t) = a m initial phase or a contracting phase that ends up in an almost constant final a(t) = a m . Now we observe that all solutions develop a bounce at some high density, effectively crossing the axis and establishing a continuous  transition from a contracting phase to an expanding one. The loitering phase, with an almost constant a(t), may last for a long period after an initial contraction but it will always end up in an expanding branch. This is illustrated in Fig. 9, where the blue dashed curve represents a phase with a(t) almost constant for a long time (very low mass). The orange curve has essentially the same future behaviour as the blue curve but its loitering phase (almost constant expansion factor) does not last so much, exhibiting a previously  contracting phase that arises abruptly. Similar behaviours arise for larger masses, and one verifies that the instabilities in the expansion factor always lead to a bounce, never finding fully collapsed solutions. This is a remarkable property of the EiBI model, since it always avoids the development of singular solutions. For completeness, the behaviour ofȧ(t) and H(t) for this case s = +1 are shown in Figs. 10 and 11. The corresponding energy density profiles appear in Figs. 12 and 13. It is amusing to see how for low mass configurations (orange solid curve, for instance) the energy density at early times is very low (large universe with very diluted energy) until it rapidly increases to reach a maximum, where it stays for some time at an almost constant value (loitering phase) with a slight decay (decompression) as we move forward in time. Then the density suddenly drops again giving rise to an expanding phase. For larger values of the scalar field mass, this process can be significantly deformed but the qualitative features remain.

Observational Constraints
In this section we shall put to observational test the EiBI cosmologies considered in the previous sections by using the latest observational data to constrain the EiBI parameter preventing the singularity in the early Universe.

BBN
In the first sections we began with the original action with the scalar field inside, and shifted the solution to other frame. This frame, as we showed, removes the singularity of the solution and for → 0 it recovers GR. Here, we put the density of the ΛCDM that produces in the "Einstein frame" a solution that still removes the singularity and recovers GR for → 0. Doing the same procedure for a theory with matter fields ρ m in the original frame gives the density: in the physical frame. For → 0 the relation givesρ = ρ, as expected.
In order to constrain the EiBI parameter we use different measurements from our Universe. The strongest one is from the Big Bang Nucleosynthesis (BBN) where z ∼ 10 9 . In principle the condition for the BBN constraint is [60,61]: In this era the matter dominance is radiation, so we take the density and the pressure to be: The density for the EiBI theory gives: From the condition Eq.(42) we get the bound 4 : 4 Recall that restoring dimensions we have → κ 2 , which represents the inverse of a matter density. Accordingly, the bound on , which represents a squared length, can be written also as | | 5.53 · 10 −8 m 2 .
where we took Ω r ∼ 10 −4 . It is important to note that this bound improves any previous existing bounds on by several orders of magnitude. In fact, the most stringent constraints so far came from elementary particles scattering experiments [62][63][64] and implied κ 2 1.86 · 10 −28 (or equivalently, 10 −2 m 2 ).
We now proceed to describe the observational data sets along with the relevant statistics in constraining the model, while we use the first constraint from the BBN above. The matter fields considered are the dark energy Ω Λ , dark matter Ω m , and radiation Ω r components. The corresponding energy density reads Inserting this equation in the expression of the Hubble factor (27), and assuming a small EiBI parameter, |R µν | −1 , yields the extended Friedman equation: where the EiBI corrections to GR solutions are apparent.

Direct measurements of the Hubble expansion
Cosmic Chronometers (CC): This data set exploits the evolution of differential ages of passive galaxies at different redshifts to directly constrain the Hubble parameter [65]. We use uncorrelated 30 CC measurements of H(z) discussed in [66][67][68][69]. The corresponding χ 2 H function reads: where H i is the observed Hubble rates at redshift z i (i = 1, ..., N) and H pred is the predicted one from the model.

Standard Candles
As Standard Candles (SC) we use measurements of the Pantheon type Ia supernova [70], that were collected in [71], as well as quasars [72] and gamma ray bursts [73]. The model parameters are fitted by comparing the observed µ obs i value to the theoretical µ th i value of the distance moduli, which are the logarithms: where m and M are the apparent and absolute magnitudes and µ 0 = 5 log H −1 0 /Mpc + 25 is the nuisance parameter that has been marginalized. The luminosity distance is defined by (here Ω k = 0, i.e., a flat space-time and E(z) is the dimensionless Hubble parameter). For the SnIa data the covariance matrix is not diagonal and the distance modulus is given as µ i = µ B,i − M, where µ B,i is the apparent magnitude at maximum in the rest frame for redshift z i and M is treated as a universal free parameter, quantifying various observational uncertainties [70]. Following standard lines, the chi-square function of the standard candles is given by where µ s = {µ 1 − µ th (z 1 , φ ν ) , ... , µ N − µ th (z N , φ ν )} and the subscript 's' denotes SnIa and QSOs.

Baryon acoustic oscillations
We use uncorrelated data points from different Baryon Acoustic Oscillations (BAO). BAO are a direct consequence of the strong coupling between photons and baryons in the pre-recombination epoch. After the decoupling of photons, the over densities in the baryon fluid evolved and attracted more matter, leaving an imprint in the two-point correlation function of matter fluctuations with a characteristic scale of around r d ≈ 147Mpc that can be used as a standard ruler and to constrain cosmological models. Studies of the BAO feature in the transverse direction provide a measurement of D H (z)/r d = c/H(z)r d , with the comoving angular diameter distance [74,75]: The angular diameter distance D A = D M /(1 + z) and the quantity D V (z)/r d with are a combination of the BAO peak coordinates above. The surveys provide the values of the measurements at some effective redshift. We employ the BAO data points collected in [76] from [77][78][79][80][81][82][83][84][85][86][87][88], in the redshift range 0.106 < z < 2.34. The BAO scale is set by the sound horizon at the drag epoch z d ≈ 1060 when photons and baryons decouple. In our analysis we used r d as independent parameter. The BAO data represent the absolute distance measurements in the Universe. From the measurements of correlation function or power spectrum of large scale structure, we can use the BAO signal to estimate the distance scales at different redshifts. In practice, the BAO data are analyzed based on a fiducial cosmology and the sound horizon at drag epoch. We calculate the redshift of the drag epoch as: where The uncorrelation of this data set yields the corresponding χ 2 : where D i is the observed distant module rates at redshift z i (i = 1, ..., N) and D pred is the predicted one from the model.

CMB distant priors
We use the CMB distant priors that based on the latest CMB measurements [2]. Its contribution in the likelihood analysis is expressed in terms of the compressed form with CMB shift parameters: where r s (z) is the comoving sound horizon at redshift z, and z * is the redshift to the photon-decoupling surface. These two CMB shift parameters together with ω b = Ω b h 2 and spectral index of the primordial power spectrum n s can give an efficient summary of CMB data for the dark energy constraints. The comoving sound horizon is given by The radiation term in the expression of E(z) for the CMB data analysis should not be ignored. It can be determined by the matter-radiation equality relation Ω r = Ω m /(1 + z eq ), and z eq = 2.5 × We assume the CMB temperature T CMB = 2.7255K. The redshift z * can be calculated by the fitting formula: where (62)

Direct Detection of the Hubble parameter
We include the latest measurement of the Hubble parameter: reported by [89]. The measurement presents an expanded sample of 75 Milky Way Cepheids with Hubble Space Telescope (HST) photometry and Gaia EDR3 parallaxes which uses the extragalactic distance ladder in order to recalibrate and refine the determination of the Hubble constant. The combination is related via the relation The χ 2 Hub estimates the deviation from the latest measurement of the Hubble constant.

Joint analysis and model selection
In order to perform a joint statistical analysis of 4 cosmological probes we need to use the total likelihood function, consequently the χ 2 tot expression is given by In order to perform a joint statistical analysis of these four cosmological probes we need to use the total likelihood function. Regarding the problem of likelihood maximization, we use an affine-invariant Markov Chain Monte Carlo sampler [90], as it is implemented within the open-source packaged Polychord [91] with the GetDist package [92] to present the results. The prior we choose is with a uniform distribution, where Ω m ∈ [0.; 1.], Ω Λ ∈ [0.; 1 − Ω m ], H 0 ∈ [50; 100]Km/sec/Mpc. For the EiBI parameter we set the range |κ 2 | ∈ [0.; 10 −33 ]. The posterior distribution of H 0 vs. the EiBI parameter is presented in Fig. 14 and 15. The Hubble parameter is H 0 = 64.99 ± 0.556km/sec/Mpc which is in between the Planck estimation of the Hubble function [2] and the latest SH0ES measurement [89]. The dark matter component is Ω m = 0.291 ± 0.0074 and the dark energy component is being Ω Λ = 0.689 ± 0.0084. The fit for the EiBI parameter gives κ 2 = (4.764 ± 3.074) · 10 −33 . We point out that the value of = 0 (which corresponds to GR) is very close to the mean value. Therefore, for | | (5.76 ± 3.53) · 10 −8 m 2 the theory is able to fit to the data while being able to successfully remove the Big Bang singularity. We point out again that this improves by several orders of magnitude the previous strongest constraints for EiBI gravity's parameter as obtained from particle physics experiments [62][63][64].
We finally point out that from the Bayes factor difference between the models that reads ∆B ij = 0.51 we get an indistinguishable difference for the ΛCDM model. Therefore, EiBI gravity with the constraint above for its parameter when coupled to a scalar yields very similar properties for the ΛCDM model, but prevents the initial singularity.

Conclusion and discussion
In this paper we have analyzed homogeneous and isotropic cosmological solutions in the context of Eddington-inspired Born-Infeld gravity coupled to a single scalar field. In the massless case we discussed the existence and properties of such solutions by direct resolution of the modified gravitational field equations, which depend on the sign of the EiBI parameter. For the bouncing solutions (s = −1) one finds that at the maximum achievable density of the scalar field the expansion factor gets to a minimum non-vanishing value (vanishing Hubble factor) undergoing a transition from a contracting phase to an expanding one (or vice versa). In the loitering solutions (s = +1) the cosmological evolution starts from an asymptotically Minkowski past with a minimum constant value of the expansion factor, such that at a given time a canonical cosmological expansion phase is triggered. These results are consistent with the expectancy raised in the seminal paper [35] and the numerical evidence found in more recent works [93] when considering perfect fluids as the matter source.
Though the massless case was known to have only two types of solutions, namely, bouncing and loitering, some kind of instability was expected for the loitering case beyond the massless scenario. Until now it was not known if they would lead to singular solutions of if they would remain nonsingular. Our analysis puts forward that the model is absolutely robust against singularities and that the loitering case gives rise to a new type of solutions in which an unstable Minkowskian phase (almost constant expansion factor) is possible in between a contraction phase and an expansion phase. The duration of this Minkowskian period (see Fig. 9) depends on model parameters and the initial conditions. No singular solutions are found in our analysis, neither in the massless nor in the massive case. All the nonsingular solutions recover the late-time GR cosmological evolution.
Getting into the details, for the s = −1 (bouncing) branch we found that for low masses the solutions closely resemble the behaviour of those of the massless case. However, for masses above a certain threshold the parametric representation of the Hubble function versus the energy density develops a kind of fish-like asymmetric fins at certain time in the cosmological evolution, a feature not present in the massless case (compare Fig. 1  that during the loitering (almost Minkowskian) phase the energy density can grow significantly while the Hubble function remains close to zero and positive (see Fig. 8). This represents a kind of intermediate state in a continuous transition from a contracting phase to a final expansion one. This behaviour is triggered by the energy density stored in the massive scalar field, whose drop marks the end of the loitering phase and the beginning of the expansion phase (see Figs. 12 and 13). To conclude, the results obtained in this paper further support the effectiveness of EiBI gravity to remove singularities in cosmological and astrophysical scenarios when coupled to different matter sources. In the present scenario with scalar fields this can be done while being consistent with the latest cosmological observations provided that the EiBI parameter is kept roughly within | | 5 · 10 −8 m 2 . One important aspect of the obtained solutions not tackled here is whether they are stable under tensorial perturbations or not. We point out that in the present case the equation of state is of the form ω = ω(t) since the energy density and the pressure of the scalar field (in the massive case) are not trivially related. In turn this should have an impact on the behaviour of the quantities prone to the development of instabilities in the tensorial perturbations equation, such as the sound speed [93,94]. Therefore, one would need to carry out such an analysis in order to further support the feasibility of this theory and their associated solutions to replace the Big Bang singularity by observationally viable nonsingular cosmologies.