Eddington-Born-Infeld cosmology: a cosmographic approach, a tale of doomsdays and the fate of bound structures

The Eddington-inspired-Born-Infeld scenario (EiBI) can prevent the Big Bang singularity for a matter content whose equation of state is constant and positive. In a recent paper we showed that, on the contrary, it is impossible to smooth a big rip in the EiBI setup. In fact the situations are still different for other singularities. In this paper we show that a big freeze singularity in general relativity (GR) can in some cases be smoothed to a sudden or a type IV singularity under the EiBI scenario. Similarly, a sudden or a type IV singularity in GR can be replaced in some regions of the parameter space by a type IV singularity or a loitering behaviour, respectively, in the EiBI framework. Furthermore, we find that the auxiliary metric related to the physical connection usually has a smoother behaviour than that based on the physical metric. In addition, we show that bound structures close to a big rip or a little rip will be destroyed before the advent of the singularity and will remain bound close to a sudden, big freeze or type IV singularity. We then constrain the model following a cosmographic approach, which is well-known to be model-independent, for a given Friedmann-Lema\^itre-Robertson-Walker geometry. It turns out that among the various past or present singularities, the cosmographic analysis can pick up the physical region that determines the occurrence of a type IV singularity or a loitering effect in the past. Moreover, to determine which of the future singularities or doomsdays is more probable, observational constraints on higher order cosmographic parameters is required.


I. INTRODUCTION
With no doubt general relativity (GR) is an extremely successful theory about to become centenary [2]. Nevertheless, it is expected to break down at some point at very high energies where quantum effects can become important, for example in the past evolution of the Universe where GR predicts a big bang singularity [3]. This is one of the motivations for looking for possible extension of GR. Moreover, it is hoped that modified theories of GR, while preserving the great achievements of GR, would shed some light over the unknown fundamental nature of dark energy or whatsoever stuff that drives the present accelerating expansion of the Universe (see Ref. [4] and references therein), said in other words: What is the hand that started recently to rock the cradle?
Indeed, several observations, ranging from type Ia Supernovae (SNeIa) [5] (which brought the first evidence) to the cosmic microwave background (CMB) [6], the baryon acoustic oscillations (BAO) [7], gamma ray bursts (GRB) [8] and measures of the Hubble parameter at different * mariam.bouhmadi@ehu.es † b97202056@ntu.edu.tw ‡ pisinchen@phys.ntu.edu.tw redshifts [9] among others, showed that the Universe has entered in the recent past a state of acceleration if homogeneity and isotropy is assumed on its largest scale. Actually, observations show that such an accelerating state is fuelled by an effective matter whose equation of state is pretty much similar to that of a cosmological constant but which could as well deviate from it by leaving room for quintessence and phantom behaviours, the latter being known to induce future singularities (see Ref. [10] and references therein). Therefore, it is of interest to formulate consistent modified theories of gravity that could appease the cosmological singularities and could shed some light over the late-time acceleration of the Universe. Of course, an alternative way to deal with dark energy singularities is to invoke a quantum approach as done on Refs. [11][12][13][14][15] A very interesting theory at this regard has been reformulated recently: the Eddington-inspired-Born-Infeld theory (EiBI) [16][17][18], as its name indicates, is based on the gravitational theory proposed by Eddington [19] with an action similar to that of the non-linear electrodynamics of Born and Infeld [20]. Such an EiBI theory is formulated on the Palatini approach, i.e., the connection that appears in the action is not the Levi-Civita connection of the metric in the theory. For a metric approach to the EiBI theory see Ref. [21]. Like Eddington theory [19], EiBI theory is equivalent to GR in vacuum, however, it differs from it in the presence of matter. Indeed while GR cannot avoid the Big Bang singularity for a universe filled with matter with a constant and a positive equation of state (with flat and hyperbolic spatial section), the EiBI setup does as shown in [18,22]. The EiBI scenario was as well proposed as an alternative to the inflationary paradigm [23] through a bounce induced by an evolving equation of state fed by a massive scalar field. This model comes with the bonus of overcoming the tensor instability previously found in the EiBI model in Ref. [24] (see also [25] for an analysis of the scalar and vectorial perturbations for a radiation dominated universe and the studies on the large scale structure formation in Ref. [26]). Black hole solutions with charged particles and the strong gravitational lensing within the EiBI theory are studied in Ref. [27]. Besides, the fulfilment of the energy conditions in the EiBI theory was studied in Ref. [28] and a sufficient condition for singularity avoidance under the fulfilment of the null energy condition was obtained. Additionally, it was shown that the gravitational collapse of non-interacting particles does not lead to singular states in the Newtonian limit [29,30]. Furthermore, the parameter characterizing the theory has been constrained using solar models [31], neutron stars [32], and nuclear physics [33]. Very recently, neutron stars and wormhole solutions on the EiBI theory were analysed in [34][35][36]. Especially in [36], the authors showed that the universal relations of the f-mode oscillation [37], which is the fundamental mode of the pulsation modes in the neutron stars, and the I-Love-Q relations [38] ,which refers to the relation among the moment of inertia, tidal Love numbers (which are parameters measuring the rigidity of a planetary body and the susceptibility of its shape to change in response to a tidal potential) and the quadrupole moment of the neutron stars, found in GR are also valid in the EiBI theory. A theory which combines the EiBI action and the f(R) action is also analysed in Refs. [39,40] (see also Ref. [41]) . A drawback of this theory is that it shares some pathologies with Palatini f(R) gravity such as curvature singularities at the surface of polytropic stars [42] (see also [43]).
We showed recently that despite the Big Bang avoidance in the EiBI setup, the Big Rip [44][45][46][47][48][49][50][51] is unavoidable in the EiBI phantom model [1]. In this paper, we will assume an EiBI model and we will carry a thorough analysis of the possible avoidance of the other dark energy related singularities, known as: (i) Sudden, Type II, Big Brake or Big Démarrage singularity [52][53][54][55], (ii) Type III or Big Breeze singularity [54][55][56][57][58], and (iii) Type IV singularity [54,55,59,60]. Those singularities can show up in GR when a Friedmann-Lemaître-Robertson-Walker (FLRW) universe is filled with a Generalized Chaplygin gas (GCG) [55] (more precisely, a phantom Generalised Chaplygin gas, or pGCG for short) which has a rather chameleonic behaviour despite its simple equation of state [55,56]. Indeed, the GCG can unify the role of dark matter and dark energy [61,62] (for a recent update on the subject see Ref. [63]), avoid the Big Rip singularity [64], describe some primitive epoch of the Universe [65,66] and alleviate the observed low quadruple of the CMB [67]. We will complete our analyses by considering as well the possible avoidance of the Little Rip event [68] on the above mentioned setup.
In the EiBI theory, there are two metrics, the first one g µν appears in the action and couples to matter, the second one is the auxiliary one which is compatible with the connection Γ [18]. The two metrics reduce to the original one in GR when the curvature term is small. Therefore, we will analyse the singularity avoidance with respect to both metrics. Furthermore, we will use the geodesic equations compatible with both metrics to study the behaviour of the physical radius of a Newtonian bounded system near the singularities. For an exhaustive analysis of the geodesics close to the dark energy related singularities in GR see Refs. [69][70][71][72]. As a result, we find that the asymptotic behaviour of g µν , more precisely the Hubble parameter and its cosmic time derivatives as defined from the metric g µν , near the singularities is consistent with that of the geodesic behaviour dictated by the same metric g µν . However, the events corresponding to the singularities with respect to g µν are usually well behaved as observed by the connection, and therefore the auxiliary metric, and so do the geodesic equations defined from the physical connection. In addition, we show that bound structures close to a big rip or little rip will be destroyed before the advent of the singularity and will remain bound close to a sudden, big freeze or type IV singularity. This result is independent of the choice of the physical or auxiliary metric.
We will further complete our analyses by getting some observational constraints on the model through the use of a cosmographic approach [73][74][75][76]. This analysis will show that the EiBI model when filled with the matter content mentioned on the previous paragraph on top of the dark and baryonic matter is compatible with the current acceleration of the Universe. The cosmographic approach relies on putting constraints on some parameters which quantify the time derivatives of the scale factor and which are called the cosmographic parameters [73][74][75][76]. These parameters depend exclusively on the space-time geometry, in this case on the geometry of a homogeneous and isotropic space-time, and not on the gravitational action or the equations of motion that describe the model (see Ref. [74] for a nice review on the subject). Hence, this approach is quite useful because given a set of constraints on the cosmographic parameters [73,76], it can be applied to a large amount of models in particular to those with relatively messy Friedmann equations like the one we need to deal with [23]. The drawback of this approach is that with the current observational data the errors can be quite large [73,[76][77][78][79][80][81]. Nevertheless, we think it is a fear enough approach for the analysis we want to carry out. Essentially, We will show that among the various birth events or past singularities predicted by the theory, the cosmographic analyses pick up the physical region which determines the occurrence of a type IV singularity (or a loitering effect) in the past, which is the most unharmful of all the types of dark energy singularities. While among the various possible future singularities or doomsdays predicted, the use of observational constraints on higher order cosmographic parameters is necessary to predict which future singularity is more probable.
The paper is outlined as follows. In section II, we shortly review the idea of the EiBI theory and present a thorough analysis on the avoidance of various singularities in this theory, through deriving the asymptotic behaviours of the Hubble parameter and its cosmic time derivatives near the singularities for both metrics (physical and auxiliary). In section III, we analyse the effects of the cosmological expansion on local bound systems in the EiBI scenario by analysing the geodesics of test particles for both metrics close to a massive body. In section IV, we use a cosmographic approach to constrain the model, and calculate the cosmic time elapsed since now to the possible, past or future, singularities. The conclusions and discussions are presented in section V. In the appendix A and B, we discuss the influence of the addition of radiation and matter on top of a pGCG fluid [64], on the destiny of the Universe in GR; as well as the singularity avoidance of the EiBI theory if the Universe is filled with a plain GCG fluid [55] which fulfils the null, strong and weak energy conditions but does not fulfil the dominant energy condition, respectively, with the purpose being shown in the following section.

II. THE EIBI MODEL AND DARK ENERGY RELATED SINGULARITIES
We start reviewing the EiBI model whose gravitational action in terms of the metric g µν and the connection Γ α µν reads [18] S EiBI (g, Γ, Ψ) = 2 κ d 4 x |g µν + κR µν (Γ)| − λ |g| The theory is formulated within the Palatini approach and therefore the Ricci tensor is purely constructed from the connection Γ. In addition, R µν (Γ) in the action (2.1) is chosen to be the symmetric part of the Ricci tensor and the connection is also assumed to be torsionless. Within the Palatini formalism we are assuming here, the connection Γ α µν and the metric g µν are treated as independent variables. The parameter κ is a constant with inverse dimensions to that of a cosmological constant (in this paper, we will work with Planck units 8πG = 1 and set the speed of light to c = 1), λ is a dimensionless constant and S m (g, Ψ) stands for the matter Lagrangian in which matter is assumed to be coupled covariantly to the metric g only. Therefore, the energy momentum tensor derived from Eq. (2.1) is conserved like in GR [18]. One can also note that the action (2.1) will recover the Einstein-Hilbert action as |κR| gets very small with an effective cosmological constant Λ = (λ − 1)/κ [18]. From now on we will assume a vanishing effective cosmological constant, i.e., λ = 1. In addition, we will restrict our analysis to a positive κ, in order to avoid the imaginary effective sound speed instabilities usually present in the EiBI theory with negative κ [32].
For a FLRW universe filled with a perfect fluid with energy density ρ and pressure p, the Friedmann equation reads [23] whereH ≡ √ κH, H is the Hubble parameter as defined from the physical metric,ρ = κρ,p = κp, and dp/dρ ≡ c 2 s denotes the derivative of the pressure with respect to the energy density. For simplicity, we will also use the following dimensionless cosmic time:t ≡ t/ √ κ where t corresponds to the cosmic time as defined from the physical metric g µν . When the curvature gets very small, i.e., |κR| ≪ |g|, the Friedmann equation (2.2) becomes where a constant equation of statep = wρ is considered 1 .
Recall that the EiBI theory recovers GR when |κR| is very small as shown in [18]. On the other hand, the conservation equation, as mentioned previously, takes the standard form It can be easily verified that the Big Bang singularity can be avoided in this theory for a radiation dominated universe [18]; i.e.p =ρ/3, and in general a universe filled with a perfect fluid with a constant and positive equation of state w; i.e., fulfilling the null energy conditions [3], bounces in the past for κ < 0 or has a loitering behaviour in the infinite past for κ > 0 [22].
Aside, we can define an auxiliary metric q µν which is compatible with the connection Γ [18]: 1 The leading order in the expansion of the scalar curvature with respect toρ satisfies κR ∝ρ at the low energy density limit, thus we can expand with respect to the energy density when the low curvature assumption is considered. and a is the scale factor of the physical metric g µν . From the auxiliary metric q µν we can define as well an auxiliary Hubble parameter H q whose rescaled dimensionless value can be expressed asH q ≡ √ κH q and reads (2.10) Notice thatH 2 q does not depend on c 2 s , unlikeH 2 in Eq. (2.2). One can see that this auxiliary Hubble parameter also recovers the Hubble parameter in standard GR as the curvature gets small: (ρ) 2 + higher order ofρ, (2.11) where w is also a constant equation of state parameter. This auxiliary metric which is compatible with the physical connection cannot avoid the Big Bang singularity in the past because both H q and dH q /dt diverge at a vanishingã and at a finite pastt, and so does the Ricci scalar defined in Eq. (2.9). We will next analyse the possible avoidance of dark energy singularities in the EiBI setup. Those singularities, as we will next review, are characterised by possible divergence of the Hubble parameter and its cosmic time derivatives at some finite cosmic time. This translates into possible divergences of the scalar curvature and its cosmic time derivatives. The EiBI model we are considering is formulated within the Palatini formalism and therefore there are two ways of defining the Ricci curvature: (i) R µν (Γ) as presented in the action (2.1) and (ii) R µν (g) constructed from the metric g µν . There are in addition four ways of defining the scalar curvature: g µν R µν (Γ), g µν R µν (g), q µν R µν (Γ) and q µν R µν (g). Therefore whenever one refers to singularity avoidance, one must specify the specific curvature one is referring to. For the dark energy singularities the important issue is the behaviour of the Hubble parameter and its cosmic time derivatives and in this case we have two possible quantities for the Hubble parameter: H related to the physical metric and H q related to the physical connection as defined in Eq. (2.8) (see also Refs. [71,72] for an alternative way of classifying dark energy related singularities).
In general, the Universe is filled with radiation, dark and baryonic matter, and dark energy: whereρ r = κρ r ,ρ m = κρ m ,ρ de = κρ de , andp de = κp de are the energy density of radiation, matter, dark energy and the pressure of dark energy, respectively. Note that p de (ρ de ) means that the equation of state of dark energy is purely a function of the dark energy density. For the sake of completeness, we will assume a universe filled with a matter contents as shown in Eq. (2.12) to go through the analysis in this paper. Note that even though dark matter and radiation are unimportant for the analysis of future singularities, they are not for the analysis of past singularities. Before starting our analysis, we will review the definition of these dark energy related singularities: • The Big Rip singularity happens at a finite cosmic time with an infinite scale factor where the Hubble parameter and its cosmic time derivative diverge [44][45][46][47][48][49][50][51] • The Sudden singularity takes place at a finite cosmic time with a finite scale factor, where the Hubble parameter remains finite but its cosmic time derivative diverges [52][53][54].
• The Big Freeze singularity happens at a finite cosmic time with a finite scale factor where the Hubble parameter and its cosmic time derivative diverge [54][55][56][57][58].
• Finally Type IV singularity occurs at a finite cosmic time with a finite scale factor where the Hubble parameter and its cosmic time derivative remain finite, but higher cosmic time derivatives of the Hubble parameter still diverge [54,[56][57][58][59][60].
To analyse the Big Freeze, Sudden, and Type IV singularities, we regard the phantom Generalized Chaplygin Gas (pGCG) as the dark energy component in this model [55,64]. Its equation of state takes the form: where α and A > 0 are two dimensionless constants. In GR, this kind of phantom energy will drive a past sudden singularity for α > 0, a future big freeze singularity for α < −1, and a past type IV singularity for −1 < α < 0 except for some quantized values of α in which the Hubble rate and its higher order derivatives are all regular in the finite past [55]. Note that the last case is different from the results shown in Ref. [55] because in that reference the authors assumed a universe filled only with a pGCG instead of the matter content given in Eq. (2.12) to which we will stick in this paper. Actually, the addition of radiation and matter contributions does not make any comparable difference to the cases in which past Sudden and future Big Freeze occur in GR, i.e., α > 0 and α < −1, respectively. However, the conclusion is different when −1 < α < 0. For the sake of convenience, we will discuss this case within the GR setup in the appendix A and in the latter subsection about the topic on Type IV singularity. After integrating the conservation equation (2.4) and assuming α > −1, one can derive the energy density of this kind of pGCG which drives the finite past Sudden or Type IV singularity in GR [55]: where a min is the scale factor corresponding to the singularity.
For later convenience, we also rewrite the energy density in terms of the scale factor as .
(2.15) Note here that we have set the scale factor at present, a 0 , as a 0 = 1 and we will use this convention in the rest of this paper. A subscript 0 stands for quantities evaluated today. On the other hand, if α < −1 and A > 0, the energy density of this pGCG which drives the finite future Big Freeze singularity in GR reads [55]: where a max is the scale factor corresponding to the future singularity.
We also rewrite the energy density in terms of the scale factor as follows for the sake of later convenience: . (2.17) Additionally, there are some special case in which the phantom character shares the same equation of state (2.13) while does not imply A > 0, as shown in Refs. [55,64]. This special pGCG will drive a finite future big freeze singularity in GR and its energy density and pressure arē where A < 0 and 1+α = 1/(2m) with m being a negative integer [55]. We will also discuss this special case within the EiBI scenario in the upcoming subsection.
A. The EiBI scenario and the Big Rip

The physical metric gµν
We showed recently that despite the Big Bang avoidance in the EiBI setup, the Big Rip singularity [44,45] is unavoidable in the EiBI phantom model [1]. Indeed, we have shown analytically and numerically that in the EiBI theory, a universe filled with matter and phantom energy with a constant equation of state w < −1 will still hit a big rip singularity; i.e. the Hubble parameter H and dH/dt blow up in a finite future cosmic time and at an infinite scale factor. Essentially, the square of the dimensionless Hubble parameterH and its cosmic time derivative near the singularity are almost linear functions of the energy density: Therefore, at very large scale factor and energy density (which grows asρ ∝ a −3(1+w) for w < −1 and constant), H and dH/dt get equally large. This happens at a finite future cosmic time [1].

The auxiliary metric qµν
As for the quantities defined by the auxiliary metric, it can be shown that and second and higher order derivatives ofH q with respect tot vanish whenρ → ∞ because their leading order in the expansion onρ is inversely proportional toρ. Furthermore, we also find that the energy density blows up andã ∝ e Hqt (2.21) whent → ∞, which corresponds to a finite t. Therefore, there is no singularity when the auxiliary metric is considered to be on the form of a FLRW metric in the EiBI theory. Indeed, the Universe approaches a de Sitter state as described by the auxiliary metric in this case. Note that according to Eq. (2.9), q µν R µν (Γ) ≈ 4/κ asρ → ∞, which is in concordance with our previous results [1].
B. The EiBI scenario and the Sudden singularity

The physical metric gµν
We seek now the possibility of smoothing the Sudden singularity that can appear in GR. We consider a pGCG fulfilling the equation of state (2.13) with α > 0 [55]. Note that in GR a universe filled with this fluid hits a past sudden singularity. The presence of matter or radiation cannot remove the occurrence of this cosmic birth on the past of the Universe. After integrating the conservation equation (2.4), one can derive the energy density of this kind of pGCG [55] which is shown in Eqs. (2.14) and (2.15).
As the Universe is filled with radiation, matter and pGCG with α > 0, the asymptotic behaviour ofH 2 and its cosmic time derivatives near the singularity (a → a min ,ρ de → 0, andρ →ρ r +ρ m →ρ ini whereρ ini is the initial dimensionless energy density at a = a min ) are the following: for α = 2; more precisely, we find that (2.25) One can see that the first cosmic time derivative of the Hubble rate blows up if α > 2 and the second order derivative of the Hubble rate blows up if α = 2 and b 4 = 0 becauseρ de vanishes at a = a min . If α = 2 and b 4 = 0, the second order cosmic time derivative is finite, but the third order derivative diverges asρ de → 0. The scale factor dependence on the cosmic time since the Universe expands from a min to a given size (at a given t) can be approximated as follows: for α = 2. One find that the Universe starts expanding from a finite past in these cases. Therefore the Universe hits a sudden singularity for α > 2 and a type IV singularity for α = 2. Furthermore, if 0 < α < 2 and α = 4/(3n + 2) in which n is a natural number, the (n + 2)-th derivative ofH will diverge even though the 1,...,(n + 1)-th derivatives are all regular. The reason is the following: the (n + 1)-th derivative ofH behaves as with C n+1 being a finite non-vanishing constant. Then the next order becomes which diverges because α > 0 and implies a type IV singularity in the finite past. If, however, 0 < α < 2 and α = 4/(3n + 2), we find that as long as α satisfies 4/(3p + 2) < α < 4/(3p − 1) with p being a positive integer, the (p+1)-th derivative of H blows up while the 1,...,p-th derivatives are all finite. This indicates a type IV singularity again.
Hence, the past Sudden singularity originally driven by a pGCG in GR will induce the following behaviours in the EiBI scenario: • If α > 2, the Universe expands from a finite past sudden singularity.
• If 0 < α ≤ 2, the Universe expands from a finite past type IV singularity.
Actually, the Sudden singularities can also occur in the future if the Universe is filled with a GCG which fulfils the strong, null, and weak energy conditions within the GR setup [55]. However, the Universe will not get into a late-time accelerating expansion phase that is observationally corroborated, so that the theory will only be worth analysing from a mathematical point of view. For completeness we present our analysis in the appendix B.

The auxiliary metric qµν
On the other hand, it can be shown for α > 0 and and second and higher order derivatives ofH q with respect tot vanish when a → a min , as well asρ de → 0, because their leading order in the expansion onρ de is proportional to (ρ de ) α/2 . Notice that in this case the auxiliary Hubble rateH q is negative whenρ de → 0 in the past becauseã = √ V a → +∞ because V → ∞ (see Eqs. (2.7) and (2.13) for a → a min ) and this happens at an infinite pastt. Furthermore, we also find that a ∝ e Hqt whent → −∞.
(2.31) Indeed, the Universe approaches a contracting de Sitter state as described by the auxiliary metric in this case. Therefore, there is no singularity of the auxiliary metric when the auxiliary metric is considered to be into a FLRW form within the EiBI theory. Note that q µν R µν (Γ) ≈ 4/κ when a → a min andã → ∞ in this case.
C. The EiBI scenario and the Big Freeze

The physical metric gµν
We seek now the possibility of smoothing the Big Freeze singularity that can appear in GR. We consider a pGCG fulfilling the equation of state (2.13) with α < −1 [55]. Note that in GR a universe filled with this fluid hits a future big freeze singularity. The presence of matter or radiation cannot remove the occurrence of this cosmic doomsday on the future of the Universe. After integrating the conservation equation (2.4), one can derive the energy density of this kind of pGCG [55] which is shown in Eqs. (2.16) and (2.17).
The asymptotic behaviour ofH 2 and the cosmic time derivatives of the Hubble parameter as a → a max ,ρ ≈ ρ de → ∞ within the EiBI setup reads for α = −3; more precisely, we find that in whichρ m andρ r denote the dimensionless energy density of matter and radiation at a = a max . It can be shown that c 0 and c 3 are always positive for any physical value of A (note thatρ r ≪ 1 at a = a max ). Nevertheless, c 2 can vanish but the first derivative ofH with respect to the cosmic time is still finite, more precisely, it vanishes whenρ blows up. It can also be shown that the scale factor dependence on the cosmic time since the Universe has a given size (at a givent) till it reaches a max is the following: One can find from Eqs. (2.36) and (2.37) that the cosmic time till the scale factor approaches a max is finite for α < −1. Therefore, a pGCG with α < −1, that leads to a big freeze in GR, fuels the following behaviour in the EiBI setup: • If α < −3, the Universe will end up into a finite future sudden singularity.
• If −3 < α < −1, the Universe will end up into a finite future big freeze singularity.
• If α = −3, the Universe will end up into a finite future type IV singularity.
In summary, as with respect to GR (α < −1) the Big Freeze singularity is smoothed in general except for −3 < α < −1 which maintains its GR character.
Additionally, there is also a finite future big freeze singularity in GR, which is driven by a very special pGCG whose energy density and pressure are shown in Eqs. (2.18), where A < 0 and 1 + α = 1/(2m) with m being a negative integer [55]. The asymptotic behaviour of H 2 and dH/dt on this case are also given by Eqs. (2.32). One can easily see that −3 < α < −1, thus the Big Freeze singularity cannot be avoided in this case.
Actually, the Big Freeze singularity can also occur in the finite past if the Universe is filled with a GCG which fulfils the strong, null, and weak energy conditions [55]. However, the Universe will not get into an accelerating expansion phase at the present time as implied by astrophysical and cosmological observations, so that the theory will only be worth to be analysed from a mathematical point of view. For completeness, we present this analysis in the appendix B.

The auxiliary metric qµν
On the other hand, it can also be shown that for α < −1 and A > 0H and second and higher order derivatives ofH q with respect tot also vanish when a → a max , i.e.,ρ → ∞, because their leading order in the expansion onρ is proportional to (ρ) (α−1)/2 and α < −1. Note that in this caseã = √ V a → +∞ and this happens at an infinite futuret. Furthermore, we also find that a ∝ e Hqt whent → ∞. (2.39) Indeed, the Universe approaches a de Sitter state as described by the auxiliary metric in this case. Therefore, there is no singularity of the auxiliary metric when the auxiliary metric is on the form of a FLRW metric in the EiBI theory. Note that q µν R µν (Γ) ≈ 4/κ when a → a max andã → ∞ in this case.
D. The EiBI scenario and the Type IV singularity

The physical metric gµν
To analyse the possibility of smoothing a type IV singularity within the EiBI theory, we consider the same kind of dark energy pGCG as shown in Eqs. (2.13) and (2.14), with −1 < α < 0. Indeed, this fluid drives a past type IV singularity in GR except for some quantized cases, i.e., if α = −n/(n + 1) with n being natural numbers, the Hubble rate and all of its cosmic time derivatives are all regular in the finite past. Note that this result is different from the one proposed in Ref. [55] because on that case the authors assumed a purely pGCG dominated universe for the analysis, which is not the case in this paper (see Eq. (2.12)). For the sake of convenience, we will briefly discuss in the appendix A how the results shown in Ref. [55] are altered under the addition of matter and radiation for −1 < α < 0 in GR.
First, if −1/2 < α < 0, the asymptotic behaviour of H 2 and the derivatives ofH asρ de → 0, a → a min , ρ →ρ r +ρ m →ρ ini , andp →ρ r /3 arē where n is a natural number and whereρ ini is the initial dimensionless energy density and ρ r is the dimensionless radiation energy density evaluated at a = a min . Furthermore, we can derive the asymptotic cosmic time behaviour near the singularity through the conservation equation Eq. (2.4) to confirm that the Universe starts to expand from a min at a finite past cosmic time for this case. Actually, a universe will start from a finite past cosmic time as long as −1 < α < 0 in the EiBI theory becausē ρ de ∝ (t −t min ) −α + higher order of (t −t min ), (2.41) for −1/2 < α < 0, and ρ de ∝ (t −t min ) 1+α + higher order of (t −t min ), (2.42) for −1 < α ≤ −1/2 whenρ de → 0 as well as a → a min . According to Eqs. (2.40), one can show that if −1/2 < α < −1/3, the first order cosmic time derivative ofH goes to infinity andH is finite, implying a finite past sudden singularity.
If α = −1/(n + 2) where n is a positive integer (note that in this case −1/3 ≤ α < 0), the Hubble rate and its higher order derivatives are all regular. The reason is the following: All the derivatives of the Hubble rate can be written as with D n being a finite non-vanishing constant, which is finite at a = a min whereρ de = 0. The next order derivative becomes which still remains finite at a = a min . We can then conclude that all the derivatives ofH are regular at a min if α = −1/(n + 2) with n being a positive integer.
On the other hand, if −1 < α ≤ −1/2 and α cannot be written as −n/(n + 1) with n being a natural number, we find that as long as α satisfies −(p + 1)/(p + 2) < α < −p/(p + 1), the p-th derivative ofH blows up while the 1,...,(p − 1)-th derivatives are all finite. This implies that a finite past type IV singularity except for a finite past sudden singularity in which the first order cosmic time derivative of the Hubble rate diverges if −2/3 < α < −1/2.
However, if α = −n/(n + 1), the Hubble rate and its higher order derivatives are all regular. The reason is the following: All the derivatives of the Hubble rate can be written as with E n being a finite non-vanishing constant, are finite at a = a min . The next order derivative hence becomes which still remains finite at a = a min . We can then conclude that all the derivatives of H are regular at a min if α = −n/(n + 1) with n being a positive integer.
On the above discussion, we have assumed that the total pressurep < 1 during the evolution of the Universe so that the left hand side of the modified field equation (2.2) is always positive. However, in some cases there may exist a particular scale factor a b satisfying a b > a min where the total pressurep = 1 at a b . Then, for this case the non-vanishing leading orders of the Hubble parameter and its cosmic time derivative in the expansion near a b are:H

The auxiliary metric qµν
On the other hand, for A > 0 it can also be shown that near a min for −1 < α < 0 and both the auxiliary Hubble rate and its firstt derivative are regular. Interestingly, we start from the conservation equation Eq. (2.4) and find that the Universe starts from a finitet as long as −1 < α < 0. More precisely, we havē ρ de ∝ (t −t min ) 1+α + higher order of (t −t min ), (2.49) for −1 < α < 0 whenρ de → 0. This fact implies it is necessary to analyse the asymptotic behaviours of the higher ordert derivatives of the auxiliary Hubble rate to see whether there is a finite past type IV singularity of the auxiliary metric or not. If −1/2 < α < 0 and α cannot be written as −1/(n+2) where n is a positive integer, we find that as long as α satisfies −1/(p + 1) < α < −1/(p + 2), the (p + 1)-th derivative ofH q blows up while the 1,...,p-th derivatives are all finite. This indicates a type IV singularity of the auxiliary metric.
However, if α = −1/(n + 2) where n is a positive integer, the auxiliary Hubble rate and its higher ordert derivatives are all regular. The reason is the following: all thet derivatives of the auxiliary Hubble rate can be written as d n dt nH q ∝ F n + (ρ de ) −α + higher order ofρ de , (2.50) with F n being a finite non-vanishing constant, which are finite whenρ de = 0. The next order derivative hence becomes d n+1 dt n+1H q ∝ F n+1 + higher order ofρ de , (2.51) which still remains finite. We can then conclude that all thet derivatives of H q are regular asρ de → 0 if α = −1/(n + 2) with n being a positive integer.
On the other hand, if −1 < α ≤ −1/2 and α cannot be written as −n/(n + 1) with n being a natural number, we find that as long as α satisfies −(p + 1)/(p + 2) < α < −p/(p + 1), the (p + 1)-th derivative ofH q blows up while the 1,...,p-th derivatives are all finite. This also implies that a past type IV singularity of the auxiliary metric.
However, if α = −n/(n + 1), the auxiliary Hubble rate and its higher ordert derivatives are all regular. The reason is the following: all thet derivatives of the auxiliary Hubble rate can be written on this case as d n dt nH q ∝ G n + (ρ de ) 1+α + higher order ofρ de , (2.52) with G n being a finite non-vanishing constant, which are finite whenρ de = 0. The next order derivative hence becomes d n+1 dt n+1H q ∝ G n+1 + higher order ofρ de , (2.53) which still remains finite. We can then conclude that all thet derivatives of H q are well defined asρ de → 0 if α = −n/(n + 1) with n being a positive integer. Thus, considering the auxiliary metric for −1 < α < 0, the results can be summarized as follows: • If α cannot be written as −1/(n + 2) or −n/(n + 1) with n being a positive integer, the Universe expands from a type IV singularity of the auxiliary metric in whichH q and dH q /dt are regular, while higher ordert derivatives of H q blow up at a finitẽ t.
As for the case in which the singularities are replaced with a loitering effect of the physical metric (p → 1) discussed in the end of previous subsubsection, the loitering effect of the physical metric corresponds to a big bang singularity of the auxiliary metric compatible with the physical connection because both H q and dH q /dt diverge at a vanishingã at a finite pastt, and so does the Ricci scalar defined in Eq. (2.9).
E. The EiBI scenario and the Little Rip

The physical metric gµν
We conclude the analysis of this section by considering as well the possibility of smoothing a little rip event within the EiBI formalism. The Little Rip event is quite similar to the Big Rip singularity except that the former happens at an infinite future while the latter at a finite cosmic time. Such an event, despite avoiding a future singularity at a finite cosmic time, will still lead to the destruction of all structures in the Universe like the Big Rip. The Little Rip has been previously analysed under four-dimensional (4D) standard cosmology [68], later on rediscovered on [54,57,82]. It can be found in dilatonic brane-world models [83] or other kind of braneworld models [84,85]. Forty years later after its discovery, the event has been baptised and named the "Little Rip" [86][87][88][89].
The simplest and the most common-used dark energy equation of state driving the Little Rip in GR is [54,82,86]p fluid (2.54), one can easily check that its energy densitȳ ρ de → ∞ as the scale factor a → ∞. The asymptotic future behaviour ofH 2 and the cosmic time derivative of the Hubble rate as a → ∞ arē Besides, for a finite a c (at a given cosmic timet c ) very close to the Little Rip event (note that a c is large enough so that the asymptotic equations (2.55) are valid), the scale factor dependence on the cosmic timet, can be approximated by whereρ dec is the dimensionless dark energy density when a = a c . As could be expected the radiation and dark matter components have no effect on the asymptotic behaviour and therefore where the Little Rip could take place. Therefore, like in GR the scale factor, Hubble parameter and its cosmic time derivatives blow up in an infinite cosmic time where the Universe would hit a little rip.

The auxiliary metric qµν
Similarly it can be shown that for a matter content given by Eq. (2.54), the asymptotic behaviours ofH 2 q and dH q /dt read and other higher ordert derivatives ofH q approach zero whenρ → ∞. Note again that q µν R µν (Γ) ≈ 4/κ when ρ → ∞. Therefore, there is no Little Rip in the EiBI theory if one regards the auxiliary metric as the FLRW metric. Actually, the Universe approaches a de Sitter state as described by the scale factor (2.57).
In summary, the asymptotic behaviours of the Universe filled with various kinds of dark energy, i.e., a phantom energy with a constant equation of state, a phantom energy driving a little rip event, and pGCG with different values of α on the basis of the EiBI theory are shown in TABLE. I.

III. THE GEODESIC ANALYSES OF A NEWTONIAN OBJECT WITHIN THE EIBI SETUP
In this section, we will consider a spherical Newtonian object with mass M and a test particle rotating around the mass M with a physical radius r. We assume that both of them are embedded in a spherically symmetric FLRW background. In Ref. [70], the authors have shown that the bound systems with a strong enough coupling in a de-Sitter background will not comove with the accelerating expansion of the Universe. However it is not the case when general accelerating phases are considered, such as the various singularities we have analysed in this paper. Therefore, we will analyse the evolution equations of its physical radius, or the geodesic equations, when the Universe approaches those singularities. In the Palatini formalism, there are two metrics, the first one g µν couples to matter, the second one q µν is the auxiliary one which is compatible with the connection and fixes the curvature of the space-time. If we regard the first metric as the one used to define the distances, then the geodesic equation is then defined by the Levi-Civita connection of g µν . On the other hand, if we consider the curvature, therefore q µν , responsible for the geodesic equations then we can define another geodesic equation expressed by the coordinatest andã defined in Eq. (2.8).
First, we regard the first metric g µν as the physical metric and the evolution equation of the physical radius reads [69,70] where the overdot denotes the cosmic time derivative and L is the constant angular momentum per unit mass of the test particle. Essentially, L satisfies the angular conservation equation [69,70]: in spherical coordinate. According to Ref. [70], theär/a term can be treated as a perturbation when the object is embedded in the de-Sitter background. However, this is not the case as the Universe approaches the Big Rip, Little Rip, Big Freeze, and the Sudden singularities because of the divergence ofä/a. In these cases, the evolution equation (3.1) takes the formr because theär/a dominates over the other terms in Eq. (3.1) [70]. There are two solutions to Eq. (3.3) because it is a second order differential equation. One solution r 1 = a(t) is the trivial solution, and the other solution can be derived directly through the following integration: past Type IV past Sudden (−2/3 < α < −1/3) past Type IV (−1 < α < 0) (1)past Type IV (α = −n/(n + 1)) (2)finite past without singularity finite past without singularity past loitering effect (a b > amin) Big Bang finite past without singularity finite past without singularity finite past without singularity (α = −n/(n + 1)) (−1 < α < 0) past loitering effect (a b > amin) Big Bang Little Rip Little Rip expanding de-Sitter TABLE I. This table summarizes how the asymptotic behaviour of a universe near the singularities in GR is altered in the EiBI theory when the Universe is filled with matter, radiation as well as phantom energy. The row labelled by (1) corresponds to −1/3 < α < 0 or −1 < α < −2/3, and where α cannot be written as −1/(n + 2) or −n/(n + 1), with n being a natural number. If α = −1/(n + 2) (−1/3 ≤ α < 0 naturally) which is labelled by (2), there is no singularity while the Universe starts to expand from a finite size at a finite cosmic time. Note that it is possible for the Universe to start from a loitering phase of the physical metric instead of a past singularities, as long as the total pressure reaches the valuep = 1 at some particular scale factor a b such that a b > amin, and it corresponds to a past big bang singularity of the auxiliary metric.
Therefore, the solution to Eq. (3.3) is the linear combination of r 1 and r 2 : Similarly, the angular motion of the particle can also be obtained by integrating Eq. (3.2), as shown in Ref. [70].

A. Dark energy with a constant equation of state: Big Rip case
As the Universe approaches the Big Rip singularity in which the asymptotic behaviours of the Hubble rate and its cosmic time derivative take the form in Eqs (2.19), the cosmic time dependence of the scale factor is similar to that in GR [64,70] Note that the exact analytical form of the previous equation was provided in Ref. [64]. Thus, we have the first trivial solution r 1 = a(t) and after integrating Eq. (3.4) we can also derive the total solution One can see that the evolution of the physical radius of the bound system is governed by the first solution in Eq. (3.7) because the second one becomes negligible as the Big Rip is approached. Therefore, the physical radius of the object will comove and diverge with the scale factor. Next, the angular motion can be obtained by integrating Eq. (3.2): where φ 0 is a constant angle from now on. Hence, one can see that φ(t) → φ 0 as the Big Rip is approached, which means that the angular motion slows down and freezes near the singularity. These qualitative descriptions of the asymptotic behaviour of the bound system near the singularity confirm the existence of the Big Rip singularity in the EiBI theory.

B. Phantom-GCG with α > 2: Sudden singularity case
If the Universe approaches a finite past sudden singularity in which the asymptotic behaviours of the Hubble rate and its cosmic time derivative take the form in Eqs. (2.22) for α > 2 (for the sake of convenience, we will only consider α > 2 even if there are other regions of the parameter space in which the Sudden singularity occurs), the cosmic time dependence of the scale factor which can be derived from Eq. (2.26) is where C S is a positive constant. Following a similar procedure as in the previous subsection, we have the first trivial solution r 1 = a(t) and after integrating Eq. (3.4) we can also derive the total solution One can see that the evolution of the physical radius of the bound system is also governed by the first solution in Eq. (3.10) because the second one becomes negligible as t → t min : which can be shown to be similar to the behaviour near the Sudden singularity in GR [55].
On the other hand, the angular motion in this case iṡ Thus, the particle starts its motion from r(t min ) = A 1 , φ(t min ) = φ 0 , with an infinite radial accelerationr at the past singularities.
C. Phantom-GCG with −3 < α < −1: Big Freeze case If the Universe approaches a finite future big freeze singularity in which the asymptotic behaviours of the Hubble parameter and its cosmic time derivative take the form in Eqs. (2.32) for −3 < α < −1 (for the sake of convenience, we will only consider −3 < α < −1 even if there are other regions of the parameter space in which the Big Freeze singularity occurs), the cosmic time dependence of the scale factor which can be derived from Eq. (2.36) is where C BF is a positive constant. Thus, we have the first trivial solution r 1 = a(t) and after integrating Eq. (3.4) we can also derive the total solution (3.14) One can see that the evolution of the physical radius of the bound system is also governed by the first solution in Eq. (3.14) because the second one becomes negligible as t → t max : which can be shown to be similar to the behaviour near the Big Freeze in GR [55]. Similarly, the angular motion of this particle near the singularities isφ consequently, Thus, the particle will remain its bound structure at r(t max ) = A 1 and φ(t max ) = φ 0 , with an infinite radial accelerationr at the future singularity.

D. Dark energy driving the Little Rip event
If the Universe approaches a little rip singularity in which the asymptotic behaviours of the Hubble rate and its cosmic time derivative take the form in Eqs. (2.55), the cosmic time dependence of the scale factor which can be derived from Eq. (2.56) is where a c is the scale factor when the Universe is close enough to the Little Rip. Thus, we have the first trivial solution r 1 = a(t) and after integrating Eq. (3.4) we can also derive the total solution where Ei[z] is the exponential integral functions [90]. We have numerically found that the second term in Eq. (3.18) proportional to A 2 vanishes near the Little Rip event, that is, whent → ∞. Thus, the evolution of the physical radius of the bound system is also governed by the first solution in Eq. (3.18) Similarly, the angular motion of this particle near the Little Rip isφ consequently, Therefore, one can see that the angular motion slows down and freezes near the Little Rip event while the radius of the bound system blows up near the Little Rip event. These qualitative descriptions of the asymptotic behaviour of the bound system near the singularity confirm the existence of the Little Rip event in the EiBI theory.
E. The Type IV singularity and the geodesic defined by the auxiliary metric As the Universe approaches a type IV singularity where the Hubble rate and its first cosmic time derivative remain finite, all the terms in Eq. (3.1) are also finite. Therefore, a bound system will remain bounded near the singularity.
As for the geodesic equations defined by the auxiliary metric q µν , we can use the coordinatet andã defined in Eq. (2.8) and similarly express the radial evolution or the geodesic equations as wherer ≡ √ V r is the physical radius corresponding to the auxiliary metric. As mentioned in the previous sections, the Hubble parameter and its derivative with respect to the cosmic timet are all finite near the singularities with respect to the metric g µν , such as the Big Rip, Big Freeze, and Sudden singularities. Therefore, the second ordert derivative ofã remains finite, so do the geodesic equations. This result, again, confirms that the auxiliary metric and the physical connection have a much smoother behaviour close to the singularities, as shown in previous section.
Before concluding this section, we would like to mention briefly an alternative method to analyse the fate of bound structure. More precisely, the motion of a test particle moving around a massive object of mass M is described by the equation of motion (3.1), or alternatively one can invoke the effective potential [69]

IV. A COSMOGRAPHIC APPROACH OF THE EIBI SCENARIO
In this section, we will use the cosmographic approach to constrain the parameters of our model, especially in the cases in which the singularities are driven by the pGCG introduced in the previous sections. The cosmo-graphic approach does not assume any particular form of the Friedmann equations and only depends on the assumption that the space-time is described by a FLRW metric. This makes this approach completely model independent so that we can use it to constrain the parameters of the pGCG model in the EiBI framework [73][74][75][76]. In Section II, one can see that the parameters in this theory have a profound influence on the doomsday or the birth of the Universe. Therefore, if the parameters in this theory are somehow constrained, one can further forecast the future evolution of the Universe and the possibility of past singularities different from the Big Bang.
The starting point of the cosmographic approach is the Taylor expansion of the scale factor a(t) with respect to the cosmic time t around the present time t 0 [73][74][75][76]: It is convenient to define the following cosmographic parameters: which are commonly called the Hubble, deceleration, jerk, snap and lerk parameters [73][74][75][76]. Furthermore, one can use the definitions in Eqs. (4.2) to derive the relations between these parameters and the redshift z derivatives of the square of the Hubble rate [76]: × (4qj + 3qs + 3q 2 j − j 2 + 4s + l). (4.3) Next, one can evaluate these quantities at the present time: where the subscript 0 denotes the quantities at the present time.
With the above equations and definitions, we can basically use the matter content given in Eq. (2.12), regarding the pGCG as the dark energy component, and rewrite the modified Friedmann equation (2.2) as a function of the redshift z then taking its z derivatives. There are six parameters in our model: κ, α, a max (or a min ), Ω m , Ω de , and Ω r where the last three are the density parameters of dark and baryonic matter, dark energy, and radiation, respectively. For the remainder of this paper, we will assume Ω r = 8.48 × 10 −5 according to Ref. [91]. Therefore, we are left with five parameters and we can in principle use Eqs. (4.4) and (H/H 0 ) 2 | z=0 = 1 to close our system and constrain our model as long as all the cosmographic parameters are given. However, one has to keep in mind that the past evolution of the Universe has imposed some physical constraints on the parameters of the model. For example, when one considers the past singularities, i.e., α > −1, the minimum scale factor a min should be very small to make this model in accordance with the wellknown evolution of the Universe. With this assumption, one may expect that these cases should be very close to the ΛCDM version of the EiBI theory as the dark energy density approaches a constant at the present time (see Eq. (2.14)). On the other hand, when one considers the future Big Freeze singularities, these physical restrictions are loosened.
In the following analysis, we will use two different methods to constrain our model, depending on which kind of singularity we analyse: (1) we will define a new dimensionless parameter in which x ≡ a s 3(1+α) where a s corresponds to the location of the singularity; i.e., it corresponds to a min or a max depending on the value of α, then we leave Y as a free variable for the sake of convenience of the computations.
(2) we can also assume that Ω m is model independent, that is, Ω m = 0.315 according to the Planck mission [6] when the models in which future singularities occur are dealt with, because there is no physical constraint on the maximum scale factor a max . Furthermore, we can assume in both approaches Ω κ ≡ 3κH 0 2 is very small according to the results in Refs. [32,33], where the authors showed that Ω κ is much smaller than the other density parameters Ω m and Ω de . In this way, we only need two cosmographic parameters q 0 and j 0 to constrain our model and we can avoid suffering from the large error bars when other cosmographic parameters s 0 and l 0 are taken into account [76,78]. In summary, our strategy will be the following: (i) method A: Even though we have five free parameters (note that Ω r has been fixed), Ω κ can be chosen as a small number [32,33]. Therefore, we are left with only four free parameters which are constrained by the Friedmann equation H/H 0 , q 0 and j 0 , leaving Y as the free parameter. (ii) method B: Again, and even though, we have five parameters, Ω κ can be chosen as a small number [32,33] and Ω m can be fixed by Planck data. Therefore, we are left again with only three free parameters which are constrained by the Friedmann equation H/H 0 , q 0 and j 0 .
Note that the cosmographic approach is a kinematic approach very useful when combined with the observational data of the current universe. In addition, the EiBI theory is very close to GR because Ω κ is very small [32,33]. In GR, one can derive the following relations from Eqs. (2.15), (2.17), and Eqs. (4.4): where X ≡ x/(1 − x) and 0 < x = a s 3(1+α) < 1. If we insert the dimensionless parameter Y = x/(1 − j 0 + 2Ω r ) defined previously and keep it as a free parameter whose value changes between 0 and 1/(1 − j 0 + 2Ω r ), in GR we can further express α, a s , and Ω m as functions of Y , q 0 , and j 0 analytically: (4.9) Note as well that Y and 1 − j 0 + 2Ω r have the same sign because 0 < Y (1 − j 0 + 2Ω r ) = x < 1. Before concluding, we would like to stress that we have 4 parameters on the GR setup: α, a s , Ω m and Ω de and three constraints: the Friedmann equation evaluated at present, the observational values of q 0 and j 0 . Therefore we are left with a unique degree of freedom or free parameter that we have chosen as Y . Before solving numerically the cosmographic constraints in the EiBI theory, we will provide some qualitative behaviours of α, a s and Ω m as functions of Y in GR (see Eqs. (4.7), (4.8), and (4.9)). First of all, one can see from Eq. (4.7) that α → +∞ (−∞) for Y → 0 + (0 − ) if 1 − j 0 + 2Ω r is positive (negative), and α → 0 for Y → 1/(1 − j 0 + 2Ω r ). Note that q 0 is always negative in an accelerating universe as it is in our case. Second, from Eq. (4.8) one can see that a s → 1 for the limits Y → 0 and Y → 1/(1 − j 0 + 2Ω r ). Note that the right hand side of Eq. (4.8) is always smaller than 1 if Y and α are positive, corresponding therefore a s to a min (See the bottom figure in FIG. 3). For negative Y and α, the values of a s defined in Eq. (4.8) can be divided into a min and a max by a particular Y whose absolute value reads |Y p |, which corresponds to 1 + α = 0. Besides, we find that a min has a local minimum for 1 − j 0 + 2Ω r > 0. Furthermore, there is a positive, divergent a max at |Y p | − corresponding to 1 + α → 0 − and a vanishing a min at |Y p | + corresponding to 1 + α → 0 + for 1 − j 0 + 2Ω r < 0. Finally, one can see from Eq. (4.9) that Ω m is a straight line ranging from (2 + 2q 0 − 4Ω r )/3 (Y → 0), which corresponds exactly to Ω m in the radiation+ΛCDM model, to 1 − Ω r (Y → 1/(1 − j 0 + 2Ω r )), which is exactly a pure radiation+CDM model.
In the following subsections, we will apply the two approaches enumerated previously just after Eqs. (4.4).
A. The first method: introducing Y The most recent cosmographic analysis based on SNeIa observational data has been carried out in Ref. [78] (as far as we know) 2 . The authors amended the conventional methodology of cosmography employing Taylor expansions of observables by an alternative method using Padé approximations, and claimed that the numerical fitting analysis for the cosmographic parameters is improved substantially by this mean. Their analysis is based on Type Ia supernovae data from the Union 2.1 compilation of the supernova cosmology project. They performed several fits distinguished by numbers (1)  Besides, the authors also showed that fits (2), (3) and (7) seem to have the most reasonable results. Note that the results in fit (7) are nearly identical to the ΛCDM model (see TABLE I of Ref. [78]) because in the ΛCDM model, q 0 = −1 + 3Ω m /2 and j 0 = 1, so we will be mainly using fit (7) in this subsection, especially for the cases in which past singularities could happen.

The analyses for positive Y in the EiBI theory
First, we use fit (7) of Ref. [78] in which q 0 = −0.561, and j 0 = 0.999 to evaluate Ω m for different Y and Ω κ in the EiBI theory. According to the definition of Y given in Eq. (4.5), only positive Y need to be considered here because 1 − j 0 + 2Ω r is positive. The results are shown in FIG. 2. The circle, block, star symbols, and the black solid line correspond to the numerical results imposing Ω κ = 10 −6 , 10 −7 , 10 −8 , and the analytical GR result described by Eq. (4.9) (Ω κ = 0), respectively. Additionally, we also include the 1σ errors of Ω m from the constraint of q 0 in this fit on the basis of the ΛCDM model, which is shown in the region between the blue lines.
From these two figures, one can obtain two simple conclusions: (i) We do not see much difference between using EiBI and GR, the reason of course is that Ω κ is very small as predicted in Refs. [32,33]. (ii) In order to obtain values of Ω m compatible with the ΛCDM model, which we will consider as a guiding line of our analysis, we will stick to small values of Y which we will consider to be smaller than 5. In FIG. 3, one can see again that the EiBI theory does not make any distinguishable difference on the results of α and a min as functions of Y from GR. Besides, though a min approaches 1 both at Y → 0 and Y → 1/(1 − j 0 + 2Ω r ) (α → ∞ and 0, respectively), there is a local minimum in the middle. We can determine the location of this minimum by taking the derivative of Eq. (4.8) with respect to Y and equating it to zero. After some calculations, we obtain ln(a min )| min = − 1 2 where Y min fulfils in which the subscript min denotes the local minimum. According to Eqs. (4.10), (4.11) and (4.7), one can see that Y min increases as j 0 gets closer to 1 + 2Ω r , with α getting closer to zero and a min approaching 0 for a given q 0 in fit (7). For example, we have found that the local minimum of a min ≈ 0.04, corresponding to α ≈ 0.1, 1 − j 0 + 2Ω r = 10 −5 and q 0 = −0.561; while the local minimum of a min ≈ 0.002 corresponding to α ≈ 0.056, 1 − j 0 + 2Ω r = 10 −9 and q 0 = −0.561. Note that the values of j 0 are compatible with fit (7) in Ref. [78] and to get our above estimations we have simply set Ω κ to zero. According to these constraints and the asymptotic behaviours analyses in previous sections (the Universe would start from a finite past sudden singularity if α > 2 and a finite past type IV singularity if 0 < α ≤ 2), a universe based on the EiBI theory may start its expansion from a past type IV singularity with both a min and α being very small, as long as j 0 is very close to 1 + 2Ω r , i.e., the radiation+ΛCDM model. It is, however, unlikely that the Universe starts from a sudden singularity on this case because this would require a relatively large value of α which would imply a too large value of a min which is incompatible with the history of the Universe.

The analyses for negative Y in the EiBI theory
In this subsubsection, we will carry out the analysis for negative Y , thus we have to assume that 1 − j 0 + 2Ω r is negative according to Eq. (4.5). We make a prior assumption that q 0 and j 0 are independent for the fit and j 0 deviates in absolute value by the same amount from the ΛCDM model as in the model discussed on the previous subsubsection, that is, q 0 = −0.561 and j 0 = 1.001. Note that with this assumption j 0 is within the 1σ errors of fit (7) of Ref. [78]. those corresponding to the Y positive case, analysed on the previous subsubsection, i.e.: (i) we do not see much difference between using EiBI and GR, the reason is that Ω κ is very small as predicted in Refs. [32,33]. (ii) In order to obtain values of Ω m compatible with the ΛCDM model, which we will consider as a guiding line of our analysis, we will stick to small values of |Y | which we will consider to be small, roughly smaller than 5 (see FIG. 4). For the sake of presenting our results in a clear and suitable way, we will highlight the region where |Y | is smaller than 3 (see FIG. 5). Similar to what we have done in the previous subsubsection, we can also numerically evaluate α, a s for different |Y | and Ω κ , then compare these results with the analytical results in Eqs. (4.7) and (4.8) in GR. The results are shown in FIG. 5. All these figures indicate again that the results in the EiBI theory are almost indistinguishable from those in GR for the Ω κ values of interest [32,33]. The first (top) figure shows that α → 0 for large |Y |, while it approaches −∞ for small |Y |. Note that because of the presence of radiation, there may be some particular scale factor a b > a min in which the total pressure satisfiesp = 1 thus implying a loitering effect at an infinite past. This particular scale factor value a b depends almost only on Ω κ and its location is shown with the horizontal blue lines in FIG. 5: lines from top to bottom correspond to decreasing Ω κ . One has to keep in mind that the value of α can determine what kind of singularity would happen in the Universe, no matter in GR or in the EiBI theory (see TABLE. I The conclusions are the following: (i) If α < −1, which implies the existence of future singularities, the values α < −3 (sudden singularities), α = −3 (type IV singularities) or −3 < α < −1 (big freeze singularities) are all compatible with the fit (7) in [78]. However, we cannot tell which of these singularities are preferred from an observational point of view and observational constraints on higher cosmographic parameters are necessary. It is worth mentioning that for a fixed α, the closer to 1 + 2Ω r the jerk parameter at present j 0 is, the larger the maximum of the scale factor at the doomsday a max would be.
(ii) If −1 < α < 0, the Universe will start either from a big loitering effect or a type IV singularity. A sudden singularity would not be allowed observationally because it will take place at a too large value of a min incompatible with the history of our universe (see FIG. 6). Moreover, it is also worth mentioning that if j 0 is getting much closer to 1 + 2Ω r , that is, the radiation+ΛCDM model, the allowable region of α enlarges because a min decreases as j 0 gets closer to 1 + 2Ω r for a fixed α. Finally, we evaluate Ω de , Ω m a s and the dimensionless cosmic time between the singularities and the current time for various α and Ω κ in TABLE. II. Note that the cases in which the past singularities are replaced with a big loitering effect are also shown.

B. The second approach: assuming Ωm
In the previous subsection, we only consider fit (7), which is the closest fit in Ref. [78] to the radiation+ΛCDM model. Additionally, we can also use other fits to constrain our model. One can see that other fits from (1) to (6) have a significant difference from fit (7): The jerk parameter j 0 is different from 1 + 2Ω r by a comparable amount. This fact makes these data sets deviate a lot from the radiation+ΛCDM model and when applied to the model we are analysing we get a too large a min which is incompatible with the history of the Universe, as we mentioned previously. Therefore, we will only analyse the models in which future singularities happen with the data sets from fits (1) to (6).
To analyse the future singularities with these data sets, we use another approach different from the one we followed on the previous subsection: we fix the value of Ω m according to the Planck mission [6] and assume it to be model independent, then we assume that Ω κ is roughly within the range 10 −7 to 10 −4 . We find that only fits (2) and (3) are compatible with the analyses of the cases in which α < −1. The reason is that for the cases in which α < −1, the derivative ofp de /ρ de with respect to the scale factor should be negative. Furthermore, the value ofp de /ρ de should be smaller than −1. On the basis of GR, these criteria are only valid in these two data sets fits (2) and (3). Interestingly, the authors of [78] also claimed that these two fits, in addition to fit (7), are the most reasonable results of their analyses. Hence, we use the data in fits (2) and (3), set the values of Ω m = 0.315, Ω r = 8.48 × 10 −5 , and Ω κ from 10 −7 to 10 −4 , and numerically solve the resulting α, Ω de , a max as well as the dimensionless cosmic time elapsed from the current time to the doomsday H 0 (t max − t 0 ). The results are shown in TABLE III and TABLE IV.  According to TABLE III and TABLE IV, one can see that fit (2) prefers the parameter space −3 < α < −1, implying a big freeze singularity for the death of the Universe, while fit (3) prefers the parameter space α < −3 where the occurrence of a sudden singularity is preferred.

V. CONCLUSIONS
The Eddington-inspired-Born-Infeld theory (EiBI) proposed recently is characterised by being equivalent to Einstein theory in vacuum but differing from it in the presence of matter. Most importantly, it also features the ability to avoid some singularities such as the Big Bang singularity in the finite past of the Universe, and the singularity formed after the collapse of a star. It is hence interesting to see whether this ability to avoid/smooth other kinds of singularities, especially, those driven by the phantom dark energy which could be responsible for the current accelerating expansion of the Universe, is efficient enough or not.
In this paper, we give a thorough analysis on the avoidance of all dark energy related singularities by deriving the asymptotic behaviours of the Hubble rate and the cosmic time derivatives of the Hubble rate defined by the physical metric g µν coupled to matter, and by the auxiliary metric q µν compatible with the physical connection. For the physical metric g µν we find that though the Big Rip singularity and the Little Rip event driven by phantom dark energy are not cured in the EiBI theory, this theory to some extent smooth the other phantom dark energy related singularities present in GR by leaving some region of the parameter space in which the future Big Freeze singularity is altered into a future sudden or future type IV singularity. Additionally, the past singularity present in GR is also smoothed in this theory in some parameter space as a past type IV singularity. Note that a past type IV singularity present in GR is, in some parameter space, worsened into a past sudden singularity while smoothed as a regular birth of the Universe at some quantized parameter space or even as a loitering effect in an infinite past. As for the auxiliary metric q µν compatible with the physical connection (we remind that the EiBI setup we are dealing with is formulatedà la Palatini formalism), all the dark energy related singularities of interest are avoided except for some very specific parameter space in which the past Type IV singularities of the auxiliary metric still exist (See TABLE. I  for a summary).
Furthermore, we analysed the fate of a bound structure near the singularities of the EiBI theory. We find that the bound structure would be destroyed before the Universe approaches a big rip singularity and a little rip event, while remains bounded at a sudden, big freeze and type IV singularities.
Besides, we also use the cosmographic approach, which is characterised by its theoretical model-independence, to constrain the parameters present in our model, so that we in principle can forecast the doomsdays and describe the birth of the Universe based on our model. As a result, it turns out that the cosmographic analyses pick up the physical region which determines the occurrence of a type IV singularity in the finite past or the loitering effect in an infinite past. While it is necessary to impose more conditions, such as the use of higher order cosmo-   TABLE II. Using the first approach in Section IV, here we show the constraints on the parameters derived for various α and Ωκ based on fit (1) of Ref. [78]. Here we assume q0 = −0.561, Ωr = 8.48 × 10 −5 , and j0 = 0.999 or 1.001, corresponding to positive or negative α. The parameter constraints under the GR framework (Ωκ = 0) are shown. Note that the cases in which the past singularities are replaced with a loitering effect in an infinite past are also shown. The additional analyses for Ωκ = 10 −40 and 10 −43 indicate that as long as Ωκ is small enough, and α is close enough to −1, it is possible to derive a small enough as to stand for the existence of a past type IV singularity in the EiBI theory. We have to stress that the allowable region of α in which a small enough as can be obtained also increases as j0 gets close to 1 + 2Ωr, as mentioned in this section.   (2) in Ref. [78] where H0 = 70.25, q0 = −0.683, j0 = 2.044. Here we use the second approach presented in Section IV in which we assume Ωm = 0.315 according to the Planck data and Ωκ = 10 −4 , 10 −5 , 10 −6 , 10 −7 . The parameter constraints under the GR framework (Ωκ = 0) are also shown.
graphic parameters with more accurate observations or other physical constraints, to forecast the future doomsdays of the Universe in this model. According to these results, the EiBI theory is indeed a reliable theory which is able to cure or smooth the singularities predicted orig-   (3) in Ref. [78] where H0 = 70.09, q0 = −0.658, j0 = 2.412. Here we use the second approach presented in Section IV in which we assume Ωm = 0.315 according to the Planck data and Ωκ = 10 −4 , 10 −5 , 10 −6 , 10 −7 . The parameter constraints under the GR framework (Ωκ = 0) are also shown.
inally in GR, thus it makes the theory a convincing alternative to GR as a way to smooth singularities.
If α < −1, the pressure has an upper bound because of the square root of the modified Friedmann equation (2.2) in the EiBI theory, so does the energy density, i.e., p ≤ 1 andρ ≤ |A| 1/α . Hence, the square of the Hubble parameter and its cosmic time derivative as the energy density approaches the maximum value arē whereρ max = |A| 1/α and δρ ≡ρ −ρ max . Note that δp = −α|A|(ρ) −α−1 δρ has been used to derive the above asymptotic equations. Because the square of the Hubble rate scales as (δρ) 2 ; i.e. in a similar way to that of a radiation dominated universe in the EiBI scenario and for which it was shown that the Big Bang singularity is replaced by a big loitering effect at an infinite past [18,22], we can conclude that the past Big Freeze singularity driven by a plain GCG in GR can also be cured through a big loitering effect in the EiBI theory.
If α > 0, the energy density is bounded from below because the pressure is bounded from above, i.e.,p ≤ 1 andρ ≥ |A| 1/α . Hence, the square of the Hubble parameter and its cosmic time derivative as the energy density approaches the minimum value can be approximated as whereρ min = |A| 1/α and δρ ≡ρ −ρ min . Note that δp = −α|A|(ρ) −α−1 δρ has again been used to derive the above asymptotic equations. Because the square of the Hubble rate also scales as (δρ) 2 ; i.e. in a similar way to that of a radiation dominated universe in the EiBI scenario and for which it was shown that the Big Bang singularity is replaced by a big loitering effect at an infinite past [18,22], we can also conclude that the future Sudden singularity driven by a plain GCG in GR can also be cured through a future big loitering effect in the EiBI theory. If −1/2 < α < 0, the situation is quite different in the EiBI setup because in this case the pressure vanishes as the singularity approaches. The asymptotic behaviours ofH 2 and the cosmic time derivative of the Hubble rate in the EiBI scenario arē Furthermore, we also find that all the higher order derivatives vanish asρ → 0 when t → ∞ because their leading  order in the expansion onρ is proportional to (ρ) α+1 (Note that 1 + α > 0 in this case). This implies the avoidance of the Type IV singularity. The results are shown in TABLE. V.
According to these results, the singularities driven by a plain GCG [55]; i.e., fulfilling the null, strong, and weak energy conditions, except for the dominant energy condition in GR, are all cured in the EiBI setup.