Inflation driven by non-linear electrodynamics

We investigate the inflation driven by a nonlinear electromagnetic field based on an NLED lagrangian density Lnled=-FfF\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {L}}_{\text {nled}} = - {F} f \left( {F} \right) $$\end{document}, where fF\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f \left( {F}\right) $$\end{document} is a general function depending on F. We first formulate an f-NLED cosmological model with a more general function fF\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f \left( {F}\right) $$\end{document} and show that all NLED models can be expressed in this framework; then, we investigate in detail two interesting examples of the function fF\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f \left( {F}\right) $$\end{document}. We present our phenomenological model based on a new Lagrangian for NLED. Solutions to the field equations with the physical properties of the cosmological parameters are obtained. We show that the early Universe had no Big-Bang singularity, which accelerated in the past. We also investigate the qualitative implications of NLED by studying the inflationary parameters, like the slow-roll parameters, spectral index ns\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_s$$\end{document}, and tensor-to-scalar ratio r, and compare our results with observational data. Detailed phase-space analysis of our NLED cosmological model is performed with and without matter source. As a first approach, we consider the motion of a particle of unit mass in an effective potential. Our systems correspond to fast-slow systems for physical values of the electromagnetic field and the energy densities at the end of inflation. We analyze a complementary system using Hubble-normalized variables to investigate the cosmological evolution before the matter-dominated Universe.


Introduction
The inflationary paradigm [1][2][3][4][5][6] of the early Universe has become a crucial part of the standard cosmological model since it has received tremendous support from the latest observational data [7][8][9][10][11].According to this scenario, the early Universe underwent an accelerated expansion that could solve several puzzles of the hot Big Bang cosmology, such as the flatness problem and the horizon problem, and provide a mechanism to generate primordial cosmological perturbations [12].The most straightforward approach to describe the inflationary era is to use a canonical scalar field with self-interacting potential.A variety of inflationary models have been proposed, such as non-minimal Higgs inflation [13], Starobinsky inflation [1] and colour red others (see [14]).Despite the impressive success of inflation, the standard cosmological model has a cosmological singularity at a finite time in the past where the curvature and energy density are not finite [15,16].Some proposals of cosmological models are free of any singularities based on various distinct mechanisms.For an incomplete list of non-singular cosmological models, see [17][18][19][20][21][22].
Another approach developed by Born and Infeld (BI) [23][24][25] as a way to cure the divergences of selfenergy of charged particles is to replace the original Maxwell Lagrangian with a nonlinear electrodynamics (NLED) Lagrangian.Similarly, Plebanski studied different models of NLED Lagrangians and proved that the BI model satisfies physically acceptable requirements [26].There are various applications of NLED in the literature, including cosmology and astrophysics [27][28][29][30][31][32][33][34], high power laser technology, plasma physics, nonlinear optics [35][36][37][38] and the field nonlinear exponential growth due to chiral plasma instability [39].In this framework, the standard cosmological model based on Friedmann-Lemaître-Robertson-Walker (FLRW) geometry with the nonlinear electromagnetic field as its source leads to a cosmological model without primordial singularity.Other interesting NLED models have been introduced in the literature .
In [66], Benaoum and Övgün have proposed a phenomenologically viable cosmological model based on NLED that could address some open cosmological problems such as the absence of primordial singularity, an early acceleration of the Universe, and the generation of matter-antimatter.One of the exciting features of nonlinearity is the removal of the initial singularity.We have assumed that a stochastic magnetic field background fills the Universe.Magnetic fields are believed to have played a crucial role in the evolution of the Universe, and it is not surprising that our Universe is teeming with magnetic fields.Magnetic fields are everywhere in our Universe [67,68].Constraints on the magnetic fields depend on the generation mechanism of the primordial magnetic field [69,70].In particular, a lower bound on the strength of the magnetic field of the order of B ≥ 3 × 10 −16 G has been obtained for intergalactic magnetic fields [71] whereas the Planck satellite in 2015 gives an upper limit to be of the order of B < 10 −9 G [72].However, very little is known about the existence and origin of magnetic fields in the early Universe [73][74][75].Finding primordial magnetic fields would transform our understanding of how our Universe evolved.
Our main aim is to study a new generalized case of NLED Lagrangian density which can be important in the very early Universe, leading to the avoidance of the singularity.To do so, in the present work, we investigate the inflation driven by a nonlinear electromagnetic field based on an NLED lagrangian density L nled = −F f (F), where f (F) is a general function of F. The nonlinearity is encoded in the function f (F).We first formulate an f -NLED cosmological model with a more general function f (F) and show that all NLED models can be expressed in this framework; then, we investigate in detail two interesting examples of the function f (F).The outline of the paper is as follows.In section 2, we present our phenomenological model based on a new Lagrangian for NLED.Solutions to the field equations with the physical properties of the cosmological parameters are obtained.Here, we show that the early Universe had no Big-Bang singularity and tended to accelerate in the past.We also investigate the qualitative implications of NLED by studying the inflationary parameters, like the slow-roll parameters, spectral index n s , and tensor-to-scalar ratio r, and compare our results with observational data.Detailed phase-space analysis of our NLED cosmological model is performed in section 3 with and without matter source.Finally, we devote section 4 to our conclusions.

General Relativity Coupled to Nonlinear Electrodynamics
The action of Einstein's gravity coupled with NLED is given as follows: where R is the Ricci scalar, L nled is the NLED Lagrangian, and we use geometrized units, where 8πG = 1, c = 1.
In general, the NLED Lagrangian can be expressed as a function of F = 1 4 F µν F µν and G = 1 4 F µν F µν , where F µν is the field strength tensor and F µν is dual.Since the classical Maxwell theory is valid in the low-energy/weakcoupling limit, the NLED Lagrangian reduces to the Maxwell one, i.e.L = −F in the corresponding limit.Here, we restrict ourselves to the case of an NLED Lagrangian depending on the electromagnetic field strength invariant F where the classical Maxwell's Lagrangian density is replaced by where f ≡ f (F) is a general function depending on F. Variation of the action for the metric and the NLED fields leads to the following field equations, and where T µν is the energy-momentum tensor of the NLED fields, From the above NLED Lagrangian density, the energy-momentum tensor can be written as: where f F = d f dF .The energy density ρ and pressure p can be obtained as follows: Assuming that the stochastic magnetic fields are the cosmic background with a wavelength smaller than the curvature, we can use the averaging of EM fields which are sources in GR, to obtain an FLRW isotropic spacetime [76].The averaged EM fields are as follows: where the averaging brackets is used for simplicity.In what follows, we consider the case where the electric field vanishes, i.e.E 2 = 0, and a non-zero averaged magnetic field leads to a magnetic Universe.Such a purely magnetic case is relevant in cosmology, where the charged primordial plasma screens the electric field, and the Universe's magnetic field is frozen for the magnetic properties to occur.
The energy density and pressure for E 2 = 0 becomes where F = 1 2 B 2 .In [66], Benaoum and Övgün have proposed a function depending on two real parameters α and β given by: β F α is dimensionless, β is the nonlinearity parameter, and the usual Maxwell's electrodynamics Lagrangian is recovered when β = 0.The energy density and pressure are: The equation of state (EoS) satisfied by this NLED Lagrangian is: which clearly shows that when the non-linearly is turned off (i.e.β = 0), it reduces to a radiation EoS.
In the context of the inflationary paradigm, we choose the background spacetime to be described by a homogeneous, isotropic and spatially flat metric, which takes the following form: where a(t) is the scale factor that governs spacetime evolution.For such a metric, the Friedmann equations can be easily computed, which results in, where H = ȧ/a is the Hubble parameter.
Using the conservation of the energy-momentum tensor ∇ µ T µν = 0, the continuity equation of the NLED is derived as, The above equation can be readily integrated, yielding the following relation between the electromagnetic field strength F and the scale factor a, It follows that: where a end (a 0 ) is the value scale factor and ) is the value of the electromagnetic field at the end of inflation (at the current time) respectively.
Notice that in geometrized units, all the quantities have a dimension of the power of length [L].In this system of units, a quantity which has L n T m M p in ordinary units converse to L n+m+p .To recover nongeometrized units, we have to use the conversion factor c m (8πG/c 2 ) p .Thus, the dimension of B 0 and H 0 is [L −1 ] in geometrized units, and the conversion factors are 1Gauss = 1.44 × 10 −24 cm −1 and H 0 = h 1.08 × 10 −30 cm −1 , where h = (67.4± 0.5) × 10 −2 and N eff = 2.99 ± 0.17 according to the Planck 2018 results [8,77].Then, we are dealing with magnetic fields of the order 10 −40 cm −1 B 0 10 −33 cm −1 in the present epoch.Then, we can obtain where we have introduced the redshift z, such that where a 0 is the present value of the scale factor (we assume a 0 = 1), a end = a(t end ) is the scale factor evaluated at the end of inflation, and z end is the redshift at the end of inflation.
Assuming that a Grand Unified Theory (GUT) describes Nature, the grand unified epoch was the period in the evolution of the early Universe that followed the Planck epoch, beginning about 10 −43 seconds after the Big Bang, in which the temperature of the Universe was comparable to the characteristic temperatures of the GUT.If the grand unification energy is taken as 10 15 GeV, this corresponds to temperatures above 10 27 K.During this period, three out of four fundamental interactions, electromagnetism, the strong, and the weak, were unified into the electronuclear force.Gravity had separated from the electronuclear force at the end of the Planck era.Physical characteristics such as mass, charge, flavour, and colour were meaningless during the grand unification epoch.The GUT epoch ended at approximately 10 −36 seconds after the Big Bang.At this point, several key events took place.The strong force separated from the other fundamental forces.Some parts of this decay process violated the conservation of baryon number and gave rise to a slight excess of matter over antimatter.It is also believed that this phase transition triggered the cosmic inflation process that dominated the Universe's evolution during the following inflationary epoch.The inflationary epoch lasted from 10 −36 seconds after the conjectured Big Bang singularity to some time values between 10 −33 and 10 −32 seconds after the singularity.During the inflationary period, the Universe continued to expand, but at a slower rate.Then, we can take as characteristic values for z end at Grand Unified Theory (GUT) scale is z end z GUT 10 28 , or two orders of magnitude less, say z end 10 −2 z GUT 10 26 , or z end 10 10 , i.e. two orders of magnitude above nucleosynthesis.Hence, for the theoretical prior z end z GUT , we have 5 × 10 31 cm −2 F end 5 × 10 45 cm −2 .For the theoretical prior z end 10 10 , we have 5 × 10 −41 cm −2 F end 5 × 10 −27 cm −2 .Taking a prior z end 10 −2 z GUT 10 26 as an educated guess we have 5 × 10 23 cm −2 F end 5 × 10 37 cm −2 .
Moreover, the energy density and the pressure in terms of the scale factor a can be expressed as, where ρ 0 = ρ B (a = 0) is the energy density at the early phase of the Universe.Plugging back ρ B in equation ( 14), it can be readily integrated, yielding the following solution for the Hubble parameter, where H end = H(a end ) is the value of the Hubble parameter at the end of the inflation.
From the above equations, it follows, lim lim It is easy to see that the Ricci scalar curvature, the Ricci tensor and the Kretschmann scalar has no singularity at early/late stages, lim lim The absence of singularities at early/late is an attractive feature peculiar to NLED.Now we concentrate on the evolution of the EoS parameter ω B = p B /ρ B , and the deceleration parameter, q, which are: It follows that, at small scale a a end , and at large scale a a end , Figure 1 shows the behaviour of the equation of state parameter ω B and the deceleration q as a function of the scale factor a for different α = 0.1, 0.5, 1.5.The square speed of sound is Assuming β F α > 0, a requirement of classical stability and causality, i.e. 0 < c 2 s ≤ 1, leads for Now, we relax the condition β ≥ 0, the region where classical stability and causality are required in Figure 2. 2 Regions of parameter space β F α vs α where a requirement of classical stability and causality is satisfied, i.e. 0 < c 2 s ≤ 1.

Cosmological Parameters
In this section, we will demonstrate that it is indeed possible to have a proper inflationary phase in the early Universe described by NLED coupled with Einstein's gravity.To describe inflation, we use the e-folds number left to the end of inflation, Then, In terms of e-folds number N, the energy density and the pressure read, where the function f 1 (N) is given by: The Hubble parameter will be: In this formalism, the slow-roll parameters are defined as: Where the first slow-roll parameter ε relates to the acceleration measure during inflation, the second slow-roll parameter η tells us how long the acceleration expression will be sustained.
The slow-roll parameters become, Therefore the tensor-to-scalar r and the scalar spectral index n s can be expressed as [78]: In figure 3, we plot the behaviour of r and n s as a function of the number of e-fold for different values of α.The Planck 2018 bounds on the spectral index and the tensor-to-scalar ratio are the following, It can be seen from the figures that the Planck 2018 result rules out the model and that the e-folding number must be large and α extremely small to achieve successful inflation, away from the theoretical prior 50 < N < 60.
To overcome the major drawback of the model given by (10), we propose the below general function f (F) depending on three real parameters A, α and β as an alternative to the function given by (10), It is clear that A = 1 3 reproduces the NLED model given by (10).The energy density of this model has the following e-folding number dependence, where From the above equation, we obtain the Hubble parameter as, One can now get the scalar spectral index and the tensor-to-scalar as, Figure 4 displays the behavior of the tensor-to-scalar ratio r versus the parameter A by fixing α = 0.5 and the parameter α by fixing A = −0.985for the e-folding number N = 50, 60 and 70.We note from the figure that the value of r increases till a maximum value and then decreases.We see that the bound r < 0.064 is achieved for −1 < A < −0.995 and −0.967 < A < 0. The figure also indicates that when α increases, the value of r decreases and the bound r < 0.064 is satisfied for α > 0.8.
In figure 5, we plot the tensor-to-scalar ratio r versus the e-folding number N for fixed α = 0.5 with A = −0.985,−0.965, −0.945 and for fixed A = −0.985with α = 0.5, 1., 1.5.The main observation from the figure is that r decreases with the e-folding number.
Similarly, figure 6 shows the variation of the spectral index n s as a function of the parameter A by fixing α = 0.5 and the parameter α by fixing A = −0.985for the e-folding number N equal to 50, 60 and 70.From the left figure in 6, the spectral index n s decreases as A increases.In the figure to the right, the spectral index n s increases to a maximum and decreases as α increases.
In figure 7, we draw the spectral index n s as a function of the e-folding number N for α = 0.5 with A = −0.985,−0.965, −0.945 and for A = −0.985with α = 0.5, 1., 1.5.In the figure to the left, for fixed value α = 0.5, increasing A leads to lower spectral index values n S .A similar tendency in the figure to the right is observed by fixing A = −0.985and increasing α.Moreover, in order to present our results more transparently, we plot in figure 8, r versus n s by fixing α = 0.5 with −1 < A ≤ −0.9 and A = −0.985with 0 < α ≤ 2 for e-folding number N equal to 50, 60 and 70.It shows that the current bounds on n s and r are satisfied.8 The spectral index n s as function of the tensor-to-scalar ratio r for α = 0.5 with A ∈ (−1, −0.9] and for fixed A = −0.985with α ∈ (0, 2].

A route to a generalization
Interestingly, one may use the Lagrangian (2) in such a way to describe all nonlinear electrodynamics models in the literature.While doing so, we want a unified prescription of all models.In the following, we provide an incomplete list of models that can be recovered from (2).

Model
Lagrangian L f (F) Övgün's exponential correction model Table 1 Incomplete List of some NLED models.
Initially, the idea is to consider that f = f (F) is an arbitrary function.Then we proceed to specify some suitable conditions the model has to satisfy.From equation ( 9), the energy density, the pressure and the equation of state parameter for purely magnetic field (i.e.E = 0 and F = 1 2 B 2 ), we have where we use the notations f F = f (F), f FF = f (F), . ... The Ricci scalar, which represents the curvature of the spacetime, is calculated by using Einstein field equations and the energy-momentum tensor, The Ricci tensor squared R µν R µν and the Kretschmann R µναβ R µναβ are also obtained, The squared sound of speed is In the following, we define the conditions that should satisfy a viable NL f ED: 1. Removal of singularities at early/late phase of the Universe.In flat spacetime, a sufficient condition for the Ricci scalar, the Ricci tensor squared, and the Kretschmann scalar to not be singular is to choose f = f (F) such that the energy density and the pressure are finite in the limit of large F and small F. Interestingly, in our approach, the pressure can be expressed in terms of the energy density and its derivative as: That implies that in order to remove singularities to the early/late phase of the Universe, both the energy ρ B density and its derivative dρ B d ln F have to be finite at a → 0 (large F) and a → ∞ (small F).

2.
Early/late time radiation/dark energy domination phase for large/small magnetic field: 3. Condition of the accelerated Universe ρ B + 3p B < 0 with the sources of NL f ED fields are, It gives, This inequality has to be satisfied for large B, that is, acceleration during the inflationary phase, and for small B, which corresponds to late-time acceleration.4. Classical stability: i) Causality of the Universe: the speed of the sound should be lower than the local light speed (c s < 1) [46].ii) To avoid the Laplacian instability, we require the conditions that must be positive (c 2 s > 0).Classical stability and causality give, which implies,

Integrability and connection with the observables
In this section, we comment on the integrability of the system at hand.Moreover, we calculate some observables in the e-folding number N.
From the first equation of Friedmann's equations (15), one can obtain an equation which shows the conservation of energy for a particle moving in an effective potential V eff (a): with For a positive scale factor, we get which can be solved by quadratures: For the general function given by (41), we obtain, As a function of N := ln(a end /a), the magnetic field B (N), the magnetic field strength F (N), the Hubble parameter H (N) and the deceleration q (N) are (recall N = 0 at the end of inflation): where Recall that at large scale a a end , N < 0. The relation between N and the redshift is Then, we have In figure 9, we illustrate the behaviour of the EoS parameter ω B and deceleration q as a function of the scale factor a for different values of α and A = −0.985.One aspect of the model is to obtain constraints on the model parameters allowed by current observations, deriving them from the constraints on the scalar spectral index n s and tensor to scalar ratio r.In slow-roll approximation, they can be expressed as [78], In our scenario, the slow-roll parameters defined in (37), can be expressed as where, according to the definition of N, (62), we have passed to derivatives with respect the redshift through For the scalar spectral index n s and tensor to scalar ratio r, we have Replacing H(z) from (63) in expressions (65), we obtain Using (41), we acquire , r(z) = 24(A + 1) which implies The Planck 2018 bounds on the spectral index and the tensor-to-scalar ratio (40), n s = 0.9649 ± 0.0042, r < 0.064 gives the parameter region {−0.992367< A ≤ −0.9897, 0.008(−375A − 251) < n s < 0.9691}, or {−0.9897 < A ≤ −0.989567, 0.008(−375A − 251) < n s < −3A − 2} or {−0.989567 < A < −0.9869, 0.9607 < n s < −3A − 2}.
Figure 10 shows constraints in the r vs n s plane for the Planck 2018 baseline analysis, and when also adding BICEP/Keck data through the end of the 2018 season plus BAO data to improve the constraint on n s (taken from [79]).The red dashed line represents a parametric plot of r vs n s given by our model (71) for the parameters A = −0.9895,α = 1.0, β = 1.0.The contour displays 1σ and 2σ confidence levels from darker to lighter.We have used data from the Planck collaboration + BICEP/Keck + BAO which is public at http://bicepkeck.org.The red dashed line represents our theoretical model that agrees within 1σ with the more stringent data.

Phase space analysis
In dynamical systems, phase space is a space in which all possible states of a system are represented.In phase space, every degree of freedom or parameter of the system is represented as an axis of a multidimensional space; a one-dimensional system is called a phase line, while a two-dimensional system is called a phase plane.Dynamical systems methods have been proven to be a powerful scheme for investigating the physical behaviour of cosmological models.As we know, there exist four standard ways of systematic investigation that can be used to examine cosmological models: (i) Obtaining and analyzing exact solutions; (ii) Heuristic approximation methods; (iii) Numerical simulation, and (iv) Qualitative analysis [80].The last case can be used with three different approaches: (a) Piecewise approximation methods, (b) Hamiltonian methods, (c) Dynamical systems methods.In approach iv (a), the evolution of the model universe is approximated through a sequence of epochs in which specific terms in the governing differential equations can be neglected, leading to a more straightforward system of equations.This heuristic approach is firmly based on the existence of heteroclinic sequences, which is a concept from iv (c).In approach iv (b), Einstein's equations are reduced to a Hamiltonian system dependent on time for a particle (point universe) in two dimensions.This approach has been used mainly for modelling and analyzing the dynamics of the Universe, nearly the Big Bang singularity (one of the approaches we will follow).In the approach iv (c) Einstein's equations for homogeneous cosmologies can be described as an autonomous system of first-order ordinary differential equations plus certain algebraic constraints.Specifically, Einstein's field equations of Bianchi's cosmologies and their isotropic subclass (FLRW models) can be written as an autonomous system of first-order differential equations whose solution curves partitioned to R n in orbits, defining a dynamical system in R n .In the general case, singular points, invariant sets, and other elements of the phase space partition can be listed and described.This study consists of several steps: determining singular points, the linearization in a neighbourhood of them, the search for the eigenvalues of the associated Jacobian matrix, checking the stability conditions in a neighbourhood of the singular points, the finding of the stability and instability sets and the determination of the basin of attraction, etcetera.On some occasions, to do that, it is necessary to simplify a dynamic system.Two approaches are applied to this objective: one, reduce the dimensionality of the system, and two, eliminate the nonlinearity.Two rigorous mathematical techniques that allow substantial progress along both lines is the centre manifold theory and normal forms.Using this approach, in [80], many results have been obtained concerning the possible asymptotic cosmological states in Bianchi and FLRW models, whose material content is a perfect fluid (usually modelling "dark matter", a component that plays an important role in the formation of structures in the Universe, such as galaxies and clusters of galaxies) with linear equation of state (with the possible inclusion of a cosmological constant).Also, several classes of inhomogeneous models are examined, comparing the results with those obtained using numerical and Hamiltonian methods.This analysis is extended in [81], to other contexts, having considered other material sources, such as the scalar fields.
Moreover, one can use tools of the theory of averaging in nonlinear differential equations and the qualitative analysis of dynamical systems to obtain relevant information about the solution's space of cosmological models.The averaging methods were applied extensively in [82][83][84][85][86][87] to single field scalar field cosmologies, and for scalar field cosmologies with two scalar fields which interact only gravitationally with the matter in [88].Within this context, one deal with perturbation problems of differential systems expressed in Fenichel's normal form [89][90][91][92][93][94][95].That is, given (x, y) ∈ R n+m and f , g smooth functions, equations can be written as: ẋ = f (x, y; ε), ẏ = εg(x, y; ε), x = x(t), y = y(t). (72 The system ( 72) is called "fast system" as opposed to which is obtained after the re scaling τ = εt, and is called "slow system".Notice that for ε > 0, the phase portraits of ( 72) and ( 73) coincide.It follows two problems that manifestly depend on two scales: (i) the problem in terms of the "slow time" variable, whose solution is analogous to the outer solution in a boundary layer problem; (ii) the fast system: a change of scale on the system which describes the rapid evolution that occurs in shorter times; analogous to the inner solution of a boundary layer problem.The solution of each subsystem will be sought in the form of a regular perturbation expansion.The subsystems will have simpler structures for singularly perturbed problems than the complete problems.Then, the slow and fast dynamics are characterized by reduced phase line or phase plane dynamics.Therefore, information on the dynamics for small values of ε is obtained.This technique is used to construct uniformly valid approximations of the solutions of perturbation problems using seed solutions that satisfy the original equations in the limit of ε → 0 [93].

Phase space analysis: pure NLED
The equation ( 55) represents the motion of a particle of the unit mass in the effective potential.This equation is satisfied on the zero-energy level, where ρ B plays the role of effective energy density parameterized through the scale factor a(t). Therefore the standard cosmological model can be represented in terms of a dynamical system of a Newtonian type: The scale factor a plays the role of a positional variable of a fictitious unit mass particle, miming the Universe's expansion.
We assume F end > 0, and introduce the variables and the time variable This system can be written in the form with effective particle-potential Thus, v 2 2 +W (u) = E, is the constant of energy.From the above system, we see that, generically, the equilibrium points of the system (77) are situated on the axis u (v = 0), and they satisfy ∂W (u)  ∂ u = 0. From the characteristic equation, it follows that just three types of equilibrium points are admitted: The eigenvalues of the Jacobian matrix evaluated at the equilibrium point with coordinate u c are −i W (u c ), i W (u c ) , such that the condition for having periodic solutions is W (u c ) > 0. Since W (u c ) = 0 at the equilibrium points, we end up with the condition 2F end f F end e −4u c + 3e 4u c f F end e −4u c < 0 as a sufficient condition for having a cyclic Universe.
Due to the relation F = F end e −4u , the above conditions can be summarized as follows: 1.The equilibrium points are given by u c = 1/4 ln(F end /F c ), where The above equation must be considered an algebraic (in most cases transcendent) equation of F c for a given f and not a differential equation for f .Furthermore, evaluated at the equilibrium point, we obtain ρ B + 3p B = 0.That is zero acceleration point.That is not unexpected since the condition for obtaining the equilibrium points is ä = 0. Evaluating the equilibrium point, we obtain The conditions for classical stability at the equilibrium points will be 2. The equilibrium point is a saddle for or, equivalently: 3. The equilibrium point is a focus for This condition leads to the existence of periodic solutions.4. It is degenerate for Using (41), the effective potential is written as System (77) becomes The equilibrium points of the system (77) are (u c , 0) such that u c = 1/4 ln(F end /F c ) where F c are the roots of (80) which is reduced to We assume F c = 0, then, we have either (for all values of α and (1 + 3A)/β > 0) or A+1) (for all values of α and . Hence, the equilibrium points are P 0 exists for (1 + 3A)/β > 0 and P 1 exists for −1 < α < 0, β < 0. From physical conditions, we assume β ≥ 0. That is, the equilibrium point P 1 is discarded.Then, for P 0 does exists we require A > −1/3, β > 0.
On the other hand, using the generating function f given by is (41), condition for classical stability and causality is (81) Then, assuming α > 0, A > −1/3, β > 0, we have that P 0 can be a nonlinear centre or nonlinear spiral (because it is nonhyperbolic).The classical stability and causality condition is also given by (91).Due to the non-hyperbolicity, we use numerical methods to investigate the stability.
For the numerics it is convenient to use the variables (F, v) through the redefinition u = 1/4 ln(F end /F).Hence, where According to our previous estimations, for the theoretical prior z end 10 −2 z GUT , we have 5 × 10 23 cm −2 F end 5 × 10 37 cm −2 .Then, we have 1.17851 × 10 −20 cm ε := 1 12 √ F end 1.17851 × 10 −13 cm, therefore, we are in the presence of a fast-slow system.The system (92) and (93) gives the dynamics in the fast manifold.That corresponds to the "horizontal motion" v = v 0 (constant), and This fact is confirmed numerically by considering the parameters  92) and ( 93) for some values in the parameter space (96).
Using the slow time T = ετ, we have the system We can easily see that the equilibrium governs the dynamics at the slow manifold points (F, v), which satisfies v = 0 and β F

Phase space analysis: NLED including matter
For the metric (13), the Friedmann equations for NLED with an extra matter source (a matter fluid with density ρ m , pressure p m and a barotropic equation of state p m = w m ρ m ) can be easily computed which results in, where H = ȧ/a is the Hubble parameter.
The conservation of the energy-momentum tensor ∇ µ T µν = 0 leads to the continuity equation of the NLED being given by (16).The continuity equation for barotropic matter is given by ρm + 3(1 + w m )Hρ m = 0.
(100) Therefore, where We have where We assume F end > 0, and using the variables (75), and the time variable τ given by (76), the system is then equivalent to where the effective potential is Now, the equilibrium points are found by solving numerically As in the previous section, a given equilibrium point u c is of one of the following types: We have the expressions As before, the eigenvalues of the Jacobian matrix evaluated at the equilibrium point with coordinate u c are −i W (u c ), i W (u c ) , such that the condition for having periodic solutions is W (u c ) > 0. Since W (u c ) = 0 at the equilibrium points, u c satisfies the equation 2F end e uc(3wm−5) (2Fend f (Fende −4uc )+e 4uc f (Fende −4uc )) + ρ m,end = 0, we end up with the condition 2F end e 2u c e 4u c (3w m − 7) f F end e −4u c − 4F end f F end e −4u c + e 10u c (3w m − 1) f F end e −4u c > 0 as a sufficient condition for having a cyclic universe.
Due to the relation F = F end e −4u , the above conditions can be summarized as follows: 1.The equilibrium points are given by u c = 1/4 ln(F end /F c ), where F c are the zeroes of the algebraic (most of the cases transcendent) equation for a given f .2. The equilibrium point is a saddle for 3. The equilibrium point is a focus for This condition leads to the existence of periodic solutions.4. It is degenerated for We apply the procedure to the present model (41).That is, using (41), the effective potential is written as System (77) becomes For the numerics it is convenient to use the variables (F, v) through the redefinition u = 1/4 ln(F end /F).Hence, where According to our previous estimations, for the theoretical prior z end 10 −2 z GUT , we have 5 × 10 23 cm −2 F end 5 × 10 37 cm −2 .Furthermore, from equations ( 101) and ( 33), we have We consider dust matter (w m = 0) and assume β = 1.Next, using H 0 = h 1.08 × 10 −30 cm −1 , where h = (67.4± 0.5) × 10 −2 , Ω m0 = 0.315 ± 0.007 and N eff = 2.99 ± 0.17 according to the Planck 2018 results [8], and considering the theoretical prior z end 10 0.840896 therefore, we are in the presence of a fast-slow system.
As before, (118) and (119) gives the dynamics in the fast manifold.That corresponds to the "horizontal motion" v = v 0 (constant), and F(τ) = F(τ 0 )e −4v 0 (τ−τ 0 ) .This fact is confirmed numerically by considering the parameters In figure 12 is presented the horizontal flow of the system (118) and ( 119) for some values in the parameter space (122).Using the slow time T = ετ, we have the system We can easily see that the equilibrium governs the dynamics at the slow manifold points (F, v), which satisfies v = 0 and

Evolution of normalized energy densities
Defining and taking F as a dynamical variable, we obtain the dynamical system  defined on the invariant surface (127).The above equation can be solved globally for Ω m , and we obtain a 2D dynamical system given by Considering the energy condition ρ m ≥ 0, ρ B ≥ 0, the phase-plane is defined by This heuristic analysis determines the dynamics of the slow manifold.Indeed, in the vacuum case, we have obtained the system (92) and (93).Moreover, according to our estimations, for the theoretical prior z end 10 −2 z GUT , we have 5 × 10 23 cm −2 F end 5 × 10 37 cm −2 .Then, we have 1.17851 × 10 −20 cm ε := 1 12 √ F end 1.17851 × 10 −13 cm, therefore, we are in the presence of a fast-slow system.The system (92) and (93) gives the dynamics in the fast manifold.That corresponds to the "horizontal motion" v = v 0 (constant), and This fact is confirmed numerically by considering the parameter region (96).For vacuum, the dynamics at the slow manifold are governed by the equilibrium points (F, v) which satisfies v = 0 and β F Depending on the parameter values, the attractors are P 0 or P 1 .They were analyzed in the coordinates (u, v).Notice that that the axis (F, v) = (0, v c ) is a line of equilibrium points, with eigenvalues {0, −4v c }. Analogously, we analyzed the matter case with dust matter (w m = 0), and we assumed β = 1.Using H 0 = h 1.08 × 10 −30 cm −1 , where h = (67.4± 0.5) × 10 −2 , Ω m0 = 0.315 ± 0.007 and N eff = 2.99 ± 0.17 according to the Planck 2018 results [8], and considering the theoretical prior z end 10 −2 z GUT , we have 4.8 × 10 17 cm −2 ρ m,end 5.2×10 17 cm −2 .We select the best-fit value ρ m,end 5.0×10 17 .Then, we have 1.17851×10 −20 cm ε := Due to the difficulties in analyzing the fast-slow dynamics and the dependence on the experimental parameters F end , ρ m,end etcetera, we have considered alternative Hubble-normalized variables (126).In the general case, the system admitted the equilibrium points (at the finite region of the phase space): (i) (Ω , F) = (0, 0), whose eigenvalues are {−4, −1+3w m }.It corresponds to the FLRW matter-dominated solution that it is a sink for w m < 1 3 or a saddle for w m > 1 3 ; and (ii) (Ω , F) = 1 f (0) , 0 .The eigenvalues are {−4, 1−3w m } that it is a saddle for w m < 1 3 or a sink for w m > 1 3 .In the particular case of f (F) in (41) we have obtained the system (134) and (135) where the late-time attractor is the FRW matter dominated solution (Ω , F) = (0, 0).The last analysis complemented the analysis of the fast-slow dynamics.

Fig. 1
Fig.1EoS of state parameter ω B and deceleration q as function of the scale factor a for different values of α.

Fig. 3
Fig.3The tensor-to-scalar r and the scalar spectral index n s as function of the e-folds number N for different values of α.

Fig. 9
Fig. 9 EoS of state parameter ω B and deceleration q as function of the scale factor a for different values of α and A = −0.985.

Fig. 10
Fig.10Constraints in the r vs n s plane for the Planck 2018 baseline analysis, and when also adding BICEP/Keck data through the end of the 2018 season plus BAO data to improve the constraint on n s (taken from[79]).The red dashed line represents a parametric plot of r vs n s given by our model(71) for the parameters A = −0.9895,α = 1.0, β = 1.0.