Compatibility of big bang nucleosynthesis in some modified gravities

This paper is devoted to investigate the implications of Einstein-Aether and modified Hořava–Lifshitz theories of gravity to the formation of light elements in the early universe named as big bang nucleosynthesis. We choose different models from these theories for a detailed investigation of big bang nucleosynthesis epoch and compare it with the observational bounds. That is, we compare the deviation of freeze-out temperature Tf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{_f}$$\end{document} with the Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document}CDM paradigm and use observational bounds on ΔTfTf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left| \frac{\Delta T_{_f}}{T_{_f}}\right| $$\end{document} to inspect constraints on the involved free parameters of these models. We apply Chi-square test on the Hubble parameter H in each model to analyze the compatibility of model parameters with the observations and find consistent results. We find that chosen models of Einstein-Aether gravity and modified Hořava–Lifshitz gravity can satisfy big bang nucleosynthesis constraints and thus constitute a viable cosmology since they can be source for dark energy sector and late-time accelerated expansion.


Introduction
The CDM cosmological model also known as standard cosmological model has been successfully explained the dynamic and evolution of the universe from its beginning (big bang) to the present accelerated phase assuming that Einstein's general relativity (GR) defines gravity. But standard cosmological model has no proofs about dark energy (DE) and dark matter (DM). It is also unable to explain matterantimatter asymmetry [1]. Researchers are of the opinion on the basis of observational data coming from Supernova Ia [2,3], cosmic microwave background [4,5] and large scale a e-mails: ams@uo.edu.pk; maliksultan23@gmail.com b e-mails: abduljawad@cuilahore.edu.pk; jawadab181@yahoo.com (corresponding author) structure formation [6][7][8] that presented accelerated expansion of the universe is due the presence of DE. The DE sector commands over the cold dark matter and navigate the universe towards accelerated expansion due to its enough negative pressure [2,9,10]. The evolution of the universe starting with big bang expanded rapidly from its initial singularity and its density exceeded from critical density in t < 10 −43 seconds (s) called Planck epoch. After that separation between fundamental forces such as gravitational and electro-nuclear forces took place in t < 10 −36 s followed by a very rapid expansion called inflationary phase completed in t < 10 −32 s. Hereafter the electroweak force decomposed into the weak nuclear and electromagnetic forces when t = 10 −12 s and at the same time there was a decrease in the energy density of the universe such that matter can exist in the form of quarks [11].
For physicists, there are many puzzles related to the early universe among which one is the formation of light elements after the explosion. Bethe in 1939 first time presented the idea of big bang nucleosynthesis (BBN) which is also known as primordial nucleosynthesis [12]. It is a prime prediction of big bang cosmology which took place after baryogenesis. Most cosmologists believe that it took place about 10-20 s after the big bang and concluded in just 20 min approximately [13]. It was time when formation of stable neutrons and protons came into exist as universe became cool enough. This cooling allowed protons and neutrons to fuse to become isotopes of Hydrogen 1 H and Helium 4 He. This phenomenon elaborates the relative abundance of Helium in the universe.
The primordial 4 He came into exist when temperature of the cosmos was T ∼ 100 MeV, while energy and number density were under the effects of relativistic leptons (electrons, protons, neutrinos) and light emitting particles photons (where protons and neutron do not contribute in the total energy density) [14]. These particles remain in thermal equilibrium retaining their collisions, while the following relationships to protons and neutrons with lepton ν e + n ↔ p + e − , (1) e + + n ↔ p +ν e , (2) n ↔ p + e − +ν e , hold the particles in thermal equilibrium (where e represents number of electrons, p describes number of protons, n is number of neutrons and ν e gives number of neutrinos in the system). Moreover, in accelerated universe, the conversion rate of protons into neutrons is used to estimate the neutron abundance.
Moreover, there are some other evolutionary era of the universe which are of much importance such as radiation dominated era, matter dominated era and the formation of large scalae structures. In the beginning, after big bang explosion the universe was so hot that energy density of electromagnetic radiation was greater than the density of matter and was driving the expansion called radiation dominated era. It take 47,000 years when energy density of matter overtaken the density of radiation and hence matter became driver of universe's evolution called matter dominated era. Furthermore, the source for the formation of larger structure was smaller structure which come into exist first and built up into the larger structure. Scientists are of the opinion that at t = 700 Myr stars (population III stars) come into exist and then dwarf galaxies and quasars yet to be detected observationally [11].
To analyze the historical background and recent picture of universe, the modified theories of gravity have much prominent role as compared to Einstein's GR. The above motivation can be fulfilled by numerous gravitational modifications. Cosmological models which deal with higher-order corrections to the Einstein-Hilbert Lagrangian have additional motivation of improving the renormalizability of GR [15,16]. These higher-order cosmological models come into exist when an additional term in Einstein-Hilbert action is added such as f (R) gravity [17] with R as Ricci curvature, Weyl gravity [18,19], Lovelock cosmology [20], Galileon theory of gravity [21,22], Einstein cubic gravity [23,24] etc. For current analysis on DE issue including modified theories of gravity, a lot of work have been done in the literature .
Many physicists worked on constraining BBN under the realm of various cosmological models to analyze the formation of light in early universe. Capozziello et al. [13] used BBN observational bounds on the primordial abundance of photons to constrain f (T ) cosmology (where T is torsion scalar) and studied three different models in order to examine BBN constraints on their free parameters. Barrow et al. [54] used BBN observational data to impose the BBN constrains on the exponent of Barrow entropy. Bhattacharjee [55] gravity (where Q is non-metricity) to investigate the viability in cosmology and found that this cosmological model can explain the observed abundances of Helium and Deuterium while Lithium problem persists. Ghoshal and Lambiase [14] used Tsallis cosmology to examine the bound on Tsallis parameter to be β < 2 by using the constraints coming from the formation of light element from the BBN observational data which permit a very small deviation from GR. Asimakis et al. [56] used BBN data to impose observational constraints on higher-order modified theories of gravity particularly Gauss-Bonnet f (G) gravity, f (P) cubic gravity and running vacuum gravity. They have given a detailed discussion of BBN epoch and investigated the deviation of freeze-out time with a comparison to CDM regime. They found BBN constraints by using the observational data on involved free parameters of various models.
The aim and motivation of this work is to address the implications of Einstein-Aether and modified Hořava-Lifshitz theories of gravity on BBN which is a source to form the light element in the early universe as the cosmological models are assumed viable if and only if they satisfy some suitable conditions imposed by BBN. Since BBN occurred in first 10-20 s after big bang and completed in approximately twenty minutes when the universe was hot and dense (indeed BBN, together with cosmic microwave background radiation, provides the strong evidence about the high temperatures characterizing the primordial universe). It reports that nuclear reactions happened in a sequence that yield the synthesis of photons [57][58][59] and therefore drive the observed universe. Thus from BBN physics, one can deduce stringent constraints on a given theory of gravity related models. Hence, in this paper, we shall confront Einstein-Aether and modified Hořava-Lifshitz theories of gravity models with BBN depending upon recent observational data. We shall consider different models related to Einstein-Aether gravity and modified Hořava-Lifshitz theories of gravity that mimic CDM model.
The arrangement of this paper is as follows: In upcoming section, we review the Einstein-Aether and modified Hořava-Lifshitz theories of gravity. In Sect. 3, we calculate the BBN constraints. In Sect. 4, we review the chi-square test and discussed observational values of the Hubble parameter H . In Sect. 5, we analyze BBN in Einstein-Aether gravity. In Sect. 6, we analyze BBN constraints in modified Hořava-Lifshitz gravity. In the last section, we conclude our findings.

Modified theories of gravity
In this section, we present two different theories of gravity. The first one is Einstein-Aether gravity in which Lorentz invariance is breakup by unit time like vector field (Aether) [60]. The other one is modified Hořava-Lifshitz gravity also known as f (R) gravity. We consider homogeneous, flat and isotropic FRW universe for which line element is where N is lapse variable depends on time t (for Einstein-Aether gravity N = 1) and a = a(t) is scale factor of the universe [61,62]. The energy momentum tensor for this model is where U a is the four velocity, ρ is the energy density and p describes pressure of the universe respectively. For both theories, continuity equation is given byρ +3H ( p +ρ) = 0, where H is the Hubble parameter. In upcoming subsections, we analyze both theories of gravity separately.

Einstein-Aether gravity
As it is believed by the physicists that Aether is a physical medium that exists everywhere in this universe homogeneously. It is the medium which is reason to travel the light from one place to the other even in vacuum. Physicists are of the opinion that Aether provides a specific static frame of reference in which every traveling object has absolute relative velocity and is well suitable for Newtonian dynamics. But Einstein in his theory of relativity rejected it by performing various experiments on optics. When the concept of cosmic microwave background come into exists, many scientists considered it as modern form of Aether. Gasperini [63] boosted up the Einstein-Aether theory which is said to be covariant moderation of general relativity [60]. The action of Einstein-Aether theory is given as [64,65] where g represents the determinant of metric tensor g μν and G is the gravitational constant. Moreover, L m is the Lagrangian density of matter while L E A is the Lagrangian density for the vector field which can be given by where M is a coupling constant, λ referred as Lagrangian multiplier and A a is a vector (tensor of rank one) to have time-like direction which satisfies the relation A a A a = −1. Furthermore, F(K ) is an arbitrary function involved due to theory with argument K which is given by where c 1 , c 2 , c 3 are constants without dimensions. Einstein field equations for this theory can be obtained from Eq. (6) as where prime denotes the derivative w.r.t the argument K , T m ab is the energy momentum tensor for matter field and is given in Eq. (5), J a b = −2K ad bc ∇ d A c and T E A ab represents energy momentum-tensor for vector field which is given by where indices (ab) represents symmetry, A a is time-like unitary vector and is defined as A a = (1, 0, 0, 0) and Y (ab) is given as The Friedmann equations for Einstein-Aether theory of gravity become where overhead dot means derivative w.r.t cosmic time 't', the constant parameter γ = c 1 + 3c 2 + c 3 and the argu- If ρ E A is the effective energy density and p E A is the effective pressure in Einstein-Aether gravity then Eqs. (14) and (15) becomes as 2K where

Modified Hořava-Lifshitz gravity
In this subsection, we will discuss the modified Hořava-Lifshitz gravity. The generalized action term for this gravity is [61,66] where √ −g = g (3) N and S m is the matter part of action. Also, where μ, λ are real numbers, L R is a function depending on three dimensional metric g (3) ab and the covariant derivatives ∇ (3) a are defined by the metric given in Eq. (4). In this scenario the argumentR takes the form The argumentR reduces to R and hence usual f (R)-gravity is obtained if we choose parameters λ = μ = 1 with flat FRW metric. If we select μ = 0,R reduces to R H L (Ricci scalar for Hořava-Lifshitz gravity) [66] and thus action (20) becomes similar to the action term of Hořava-Lifshitz-like f (R)-gravity [67]. Hence, this assumption (μ = 0) corresponds to some degenerate limit of general f (R) Hořava-Lifshitz gravity. We call this limit degenerate as it is very difficult to obtain (might be impossible). Considering FRW cosmology for action (20) for which the spatial curvature R (3) ab = R (3) vanishes and thus, L (3) R does not contribute any thing (same FRW cosmology obtained for any choice of L . It is obvious that this situation varies when black hole or solutions with non-trivial dependence are considered. Suppose that universe is composed with perfect fluid, by varying (20) w.r.t g (3) ab and setting N = 1, the Friedmann equations for modified Hořava-Lifshitz gravity are given by where prime denotes the derivative of respective function with respect to argument. Here c is the constant of integration and the term ca −3 represents the dark matter part when c > 0 [68]. As during BBN, dark matter part vanishes (i.e. c = 0) [66], therefor Friedmann equations (23) and (24) become The value of argumentR from Eq. (22) reduces tõ

BBN constraints: basic scenario
In this section, we will develop the relation for deviation of freeze-out time T f and BBN constraints. Since BBN noticed in radiation era [58,69,70], thus during BBN for standard model radiation in context of GR, the first Friedmann equation can be approximated as where is reduced Planck mass, ρ r is energy density of a relativistic particles filling the universe and is given as Here T represents the temperature and g * = g(T ) is the effective number of the degree of freedom usually approximated as g * ∼ 10. Substituting the values of ρ r and M p in Eq. (28), we get where M p = √ 8π M P is referred to as Plank mass. Due to radiation conservation, the scale factor of universe evolves as where t is the cosmic time. The Hubble parameter H in terms of cosmic time t becomes as H = 1 2t which leads to relation between temperature T and cosmic time t Since number of neutron arises during BBN due to conversion of some protons-neutron rate [70,71] which is given as and its inverse np (T ). Thus the total rate can be given as Considering all these particles at same temperature which is low enough causes to use Boltzmann distribution instead of the Fermi-Dirac distribution. Assuming electron mass negligible as compared to electron and neutrino energies, some simple mathematics can be used to calculate neutron abundance via conversion rate of proton into neutron which gives the relation for total rate as [13,54,72,73] where Q is the mass difference between neutron and proton given as Q = m n − m p = 1.29 × 10 −3 GeV and gives the neutron to proton equilibrium ratio and τ = 8803 ± 1.1 s is the neutron average life time [74]. Moreover, the function λ(T f ) represents the portion of neutrons that decay into proton in the interval [T f , T n ]. Now we compare the universe expansion rate H −1 and the function tot (T ) given in (32) to calculate the freeze-out time T f . We can consider system in thermal equilibrium if interaction time is much larger then expansion time (i.e. H −1 tot (T )) [57,69]. On the other hand, the particles decouple if expansion time is much larger then interaction time (i.e. H −1 tot (T )). The temperature during which decoupling takes place is freeze-out temperature [13,54,72,73]. Using Eqs. (30) and (32), above requirement leads to the relation The Hubble parameter H will deviate from H G R in context of modified cosmology, due to which freeze-out time also presents a deviation T f from (28). Due to this deviation in fractional mass, Y p is given as where we take T (T n ) = 0 as T n is fixed by deuterium binding energy [72,73,75]. During BBN era, the observational estimation of mass fraction is [76][77][78][79][80][81] One can obtain an extra term in Friedmann equations by using any modified theories of gravity which needs to be small compared to the radiation sector of standard cosmology during BBN era so that observational facts may not be spoiled. Thus, from the general modified Friedmann equation, 3M 2 where H G R = M p ρ r 2 is the rate at which universe expands in standard cosmology. Thus, we obtain This deviation from standard cosmology will lead to a deviation in the freeze-out temperature T f . Since H G R = tot ≈ c q T 5 f . This relation along with (33) leads to the relation In the regime ρ DE ρ r , one can find finally The above theoretically evaluated relation should be compared with the following observational bound which is found by using observation estimation of the baryon mass fraction converted into 4 He [76][77][78][79][80][81] which is given in Eq. (35).

Chi-square test
A statistical test Chi-square developed by Pearson in 1900 [82] is used for a comparison between the observed values and expected values. Pearson also applied it to test the goodness of fit for frequency curves. This test exhibits, whether the difference between two data sets is due to chance or due to relationship of variables under consideration. The condition to perform Chi-square test is that data set must have points more than five. The outcome of Chi-square test is a single number which analyze the difference between the observed values and expected values of the data sets. There is no difference between two data sets (data sets are identical) if the value of Chi-square is zero. A larger value of χ 2 exhibits bigger difference between the two data sets as compared to a smaller value of χ 2 . To evaluate the limit and best fitting values for different models, we use the χ 2 static as where i counts the data points, O i represents the observed value, E i is the corresponding expected value, σ i is the error associated with ith observed value in the data and p denotes the set of model parameters. To observe the compatibility of models with the recent observational data, we will calculate the Hubble parameter H for each model in upcoming sections and compare the corresponding values of H with recent data. Bhardwaj et al. [83] presented the observed values of H for 46 different points against redshift parameter z along with their observational error evaluated by using various age approach by different cosmologists presented in the Table 1 as

BBN of various models in Einstein-Aether gravity
In this section, we apply the formalism obtained in the Sect. 3 to find the BBN constraints under the realm of Einstein-Aether gravity. We use DE relation (18) which holds in the contexts of Einstein-Aether gravity to find the BBN constraints as well as the Hubble parameter H . We will focus on three different models to examine the BBN constraints and Hubble parameter.

Model 1
The first model which we use to inspect BBN constraints in context of Einstein-Aether gravity is where f 0 is a positive constant as cosmic acceleration is observed only for f 0 > 0 [100] and n is the only free parameter. Moreover, f 0 can be expressed in terms of H 0 (i.e. present value of Hubble parameter) as n = 1. The present value of DE density can be given as where DE0 is the current matter density parameter and is approximately equivalent to 0.7 while H 0 = 73.02 ± 1.79 km/(sMpc), ∼ 2.1 × 10 −42 GeV.
is the Hubble parameter at present [13]. The best fit on other free parameters can be obtained by taking the CC + H 0 + SNeIa + BAO observational data [101]. During BBN era, assuming ρ DE = ρ DE0 and inserting its value along with F(K ) in Eq. (18), we get Inserting Eq. (42) along with (44) into (18) and then simplifying (39), we have where T f is given in (33) and We depict Fig. 1 against the free parameter n appearing in (45) with an upper bound coming in (40). It can be observed from the figure that the mathematical relation (45) satisfies the bound (40) The above equation is a complicated first order ordinary differential equation from which it is difficult to obtain the exact solution. To overcome this complexity, we solve the above equation numerically by using the software mathematica and extract the 43 different values of H against z given in serial number 1 to 43 in Table 1 and plots are in Fig. 2. The error bars representing observations given in [83] are The Chi-square analysis (to test the goodness of fit for the curves) depicts that the obtained results are closed to the observations when n = 1.085 and differ utmost from the observational data when n = 1.1. Moreover, the values of model parameters which are estimated to analyze the BBN compatibility are best fit according to χ 2 test.

Model 2
The second model which we inspect for BBN constraints under the realm of Einstein-Aether gravity is given as where α and β = 0 are real constants. This model is more generalized as the model given in (42) and reduces to the model (42) if we choose β = 0. We can find mathematical relation for the constant α for BBN era as Substituting Eq. (49) along with (50) into (18) and then simplifying for the constraints given in (39), we have We plot Substituting the corresponding values in the above equation and using the transformation from t to z, simplification leads to the equation We solved the above equation numerically using the software mathematica due to its complexity and extracted 43 different values of Hubble parameter H for various values of z given in serial number 1 to 43 in Table 1 and plotted in Fig. 4. The error bars added to the figure for comparison representing observations given in [83].

Model 3
The third model which we assume to investigate BBN constraints under the realm of Einstein-Aether gravity is much more generalized as compared to both of the previous models and is given as where α and β are real constants. As in BBN era, DM does not exist. Thus the mathematical relation for constant β in BBN era can be obtained as Using the transformation from t to z and differentiating w.r.t. redshift z, simplification leads to the relation

Fig. 5 Variation of
T f T f given in (56) against free parameter n for the model F(K ) = K + α K 2 + β K n in Einstein-Aether cosmology We proceed the above equation numerically by using mathematica software as it is a complex first order differential equation and extracted 43 different values of Hubble parameter H for various values of z given in serial number 1 to 43 in Table 1 and plot is in Fig. 6. The error bars representing the recent observational data given in [83] are also added to the figure for comparison. The The Chi-square analysis depicts that the values of Hubble parameter are closed to the observations when n = 1.0358 and differ utmost from the observational data when n = 1.0374. It can be seen from the χ 2 analysis that the model parameters estimate best to analyze BBN compatibility with the observational data and the CDM regime.

BBN in modified Hořava-Lifshitz gravity
In this section, we use the formalism to impose the BBN constraints obtained in Sect. 3 under the realm of modified Hořava-Lifshitz gravity. We use DE relation (25) which holds in context of modified Hořava-Lifshitz gravity. We will concentrate on distinct models to examine the BBN constraints.

Model 1
The first model with which we inspect BBN constraints in context of modified Hořava-Lifshitz gravity is given as where α is a real constant andR is given in (27). Since DM does not exist during BBN era, we can get expression for the constant α as Inserting Eqs. (59) and (60) into (25), the simplification of Eq. (39) yields  Substituting the transformation from t to z and differentiating w.r.t. redshift z, simplification leads to the relation The above equation is a complicated second order differential equation from which it is difficult to obtain the exact solution.
To overcome this complexity, we solve this equation numerically by using the software mathematica and extracted 43 different values of H against z given in serial number 1 to 43 in Table 1 and plot them in Fig. 8 for different choices of the parameter μ. The error bars representing observations given in [83] are  The second model which we use to investigate BBN constraints under the realm of modified Hořava-Lifshitz gravity is where α and β are real constants. During BBN era when DM does not exist, the value of constant β can be obtained as Inserting Eqs. (64) and (65) into (25), the simplification of Eq. (39) yields In Fig. 9, we present the graph of Substituting the transformation from t to z and differentiating w.r.t redshift z, simplification leads to the relation We solved the above equation numerically using the software mathematica due to its complexity and calculated first 43 values of H against z listed in serial number 1 to 43 in Table 1 and plotted in Fig. 10. The error bars representing observations given in [83] are also added to the figure for comparison.  The Chi-square analysis (to test the goodness of fit for the curves) gives that the obtained results are closed to the observations [83] when H 0 = 63 and differ utmost from the observational data when H 0 = 72. It is easy to deduct from the χ 2 analysis that the values of model parameters H (0) = 1, λ = 0.334, μ = 0.0001 and α = 0.001, ω = −1/3, DE0 = 0.7 estimate best to analyze BBN compatibility with the observational data and the CDM regime.

Concluding remarks
In this paper, BBN phenomenon has been investigated in Einstein-Aether and modified Hořava-Lifshitz theories of gravity in the presence of observational bounds on the primordial abundance of 4 He. These theories efficiently describe the late time accelerated expansion of the universe and do not spoil the behavior of early universe, particularly BBN era. Firstly, we have investigated BBN for three different models of Einstein-Aether gravity. We have attained the

BBN bound
T f T f < 4.7 × 10 −4 as per recent observations for F(K ) = f 0 K n (with n ≤ 0.824), F(K ) = α K n +β (with n ≤ 0.752) and F(K ) = K + α K 2 + β K n (with n ≤ 0.807, γ = 0.002 and α = 10 −9 ) as shown in Figs. 1, 3 and 5 respectively. It is of great interest that all three models described the CDM regime for free parameter n → 0. Moreover, Chi-square analysis has been performed for the same three models to estimate the values of model parameters and have been plotted in Figs. 2, 4 and 6 along with error bars of observational data [83]. The estimated values of the parameters for the model F(K ) = f 0 K n are H 0 = 0.69, M = 5, γ = 1, for the model F(K ) = α K n +β are n = 1.096, M = 5, β = 0.6, γ = 1 and for the model F(K ) = K + α K 2 + β K n are H 0 = 0.69, α = 10 −9 , M = 1.2, γ = 0.002. Chi-square test showed that obtained results of Hubble parameter H are very near to the observational data.
Secondly, we have investigated BBN phenomenon for two models of modified Hořava-Lifshitz gravity. We have obtained BBN bound for both models F(R) = αR n (with n ≤ 0.824) and F(R) = αR + βR n (with n ≤ 0.824) as shown in Figs. 7 and 9, respectively. Moreover, both models described CDM regime if we restrict the free parameter n near zero. Additionally, Chi-square has computed for the same two models for an estimation of parameters values existing in the model which have been plotted in Figs. 8  T f T f < 4.7 × 10 −4 confirms the primordial abundance of Helium ( 4 He) which is a confirmation to the idea of photon (light-element) production in the early history of universe. In Figs. 2, 4, 6, 8 and 10, the Chi-square test confirms the estimation of parameter values presented in the models of Einstein-Aether and modified Hořava-Lifshitz theories. In future, we can study the effects of Einstein-Aether and modified Hořava-Lifshitz theories of gravity on the primordial gravitational wave backgrounds. Since its observational evidence of number of existing polarizations are a successful tool for testing GR and extended theories of gravity. We expect that our contribution will be a source for cosmologists to explore more about BBN as it is an interesting addition in modified theories of gravity.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This work is theoretical, so no data has been generated or analyzed. Hence, data sharing is not applicable in this manuscript.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 . SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.