**4. BBN Constraints on** *f* (*T***,** *TG*) **Gravity**

In this section, we will apply the BBN analysis in the case of *f*(*T*, *TG*) gravity. Let us mention here that in general, in modified gravity, inflation is not straightaway driven by an inflaton field, but the inflaton is hidden inside the gravitational modification, i.e., it is one of the extra scalar degrees of freedom of the modified graviton. Hence, in such frameworks reheating is usually performed gravitationally, and the reheating and BBN temperatures may differ from standard ones. Nevertheless, in the present work we make the assumption that we do not deviate significantly from the successful concordance scenario, in order to examine whether *f*(*T*, *TG*) gravity can at first pass BBN constraints or not. Clearly a more general analysis should be performed in a separate project, to cover more radical cases too. In the following, we will examine five specific models that are considered to be viable in the literature.

#### *4.1. Model I: f* = −*T* + *β*<sup>1</sup> p *T*<sup>2</sup> + *β*2*T<sup>G</sup>*

Firstly, we investigate the model *f* = −*T* + *β*<sup>1</sup> p *T*<sup>2</sup> + *β*2*T<sup>G</sup>* [28]. Since in our analysis we focus on the radiation era where the Hubble parameter *H*(*t*) = <sup>1</sup> 2*t* , we can express the derivatives of the Hubble parameter as powers of the Hubble parameter itself, e.g., *H* ˙ <sup>=</sup> <sup>−</sup>2*H*<sup>2</sup> and *<sup>H</sup>* ¨ = 8*H*<sup>3</sup> . Additionally, in order to eliminate one model parameter we will apply the Friedmann equation at present time, requiring

$$
\Omega\_{\rm DE0} \equiv \rho\_{\rm DE0} / (\Im M\_{\rm P}^2 H\_0^2) \,\text{.}\tag{27}
$$

where Ω*DE* is the dark energy density parameter and with the subscript "0" denoting the value of a quantity at present time. Doing so, and inserting *f* = −*T* + *β*<sup>1</sup> p *T*<sup>2</sup> + *β*2*T<sup>G</sup>* into (13) and then into (25), we finally find

$$\begin{split} \frac{\Delta T\_f}{T\_f} = (10c\_q T\_f^3)^{-1} \zeta H\_0 \Omega\_{DE0} (3 - 2\mathfrak{H}\_2)^{-3/2} \\ &\quad \cdot \left( 9 - 15\mathfrak{H}\_2 + 6\mathfrak{H}\_2^2 \right) \left[ (3 + 2\mathfrak{H}\_2) H\_0^2 + 2\mathfrak{H}\_2 H\_0 \right]^{3/2} \\ &\quad \cdot \left[ \left( 9 + 3\mathfrak{H}\_2 - 2\mathfrak{H}\_2^2 \right) H\_0^4 + 9\mathfrak{H}\_2 H\_0^2 H\_0 + \mathfrak{H}\_2^2 H\_0 \tilde{H}\_0 \right]^{-1}, \end{split} \tag{28}$$

where

$$\zeta \equiv \left(\frac{4\pi^3 g\_\*}{45}\right)^{\frac{1}{2}} M\_{Pl.}^{-1} \,. \tag{29}$$

In this expression, we insert [92]

$$
\Omega\_{\rm DE0} \approx 0.7, \quad H\_0 = 1.4 \times 10^{-42} \,\text{GeV}, \tag{30}
$$

and the derivatives of the Hubble function at present are calculated through *H* ˙ <sup>0</sup> <sup>=</sup> <sup>−</sup>*H*<sup>2</sup> 0 (1 + *q*0) and *H* ¨ <sup>0</sup> = *H*<sup>3</sup> 0 (*j*<sup>0</sup> + 3*q*<sup>0</sup> + 2) with *q*<sup>0</sup> = −0.503 the current decceleration parameter of the Universe [92], and *j*<sup>0</sup> = 1.011 the current jerk parameter [93,94]. Hence, *H* ˙ <sup>0</sup> ≈ −9.7 × 10 <sup>−</sup><sup>85</sup> GeV 2 and *H* ¨ <sup>0</sup> ≈ 4.1 × 10 <sup>−</sup><sup>126</sup> GeV 3 .

Using the BBN constraint (26), we conclude that *β*<sup>2</sup> ∈ (−2.98, −2.93) ∪ (0.99, 1.01), where we have used (27) to find

$$\begin{split} \boldsymbol{\beta}\_{1} = \sqrt{3} H\_{0} \boldsymbol{\Omega}\_{\rm DE0} \Big[ (\mathbf{3} + 2\boldsymbol{\beta}\_{2}) \boldsymbol{H}\_{0}^{2} + 2\boldsymbol{\beta}\_{2} \dot{\boldsymbol{H}}\_{0} \Big]^{3/2} \\ &\quad \cdot \Big[ \Big( 9 + 3\boldsymbol{\beta}\_{2} - 2\boldsymbol{\beta}\_{2}^{2} \Big) \boldsymbol{H}\_{0}^{4} + 9\boldsymbol{\beta}\_{2} \boldsymbol{H}\_{0}^{2} \dot{\boldsymbol{H}}\_{0} + \boldsymbol{\beta}\_{2}^{2} \boldsymbol{H}\_{0} \dot{\boldsymbol{H}}\_{0} \Big]^{-1} . \end{split} \tag{31}$$

Using the above range of *β*2, we find that *β*<sup>1</sup> ∈ 2.09 × 10 −26 , 0.001 ∪ (1.380, 1.384). In Figure 1, we depict |∆*T<sup>f</sup>* /*T<sup>f</sup>* | appearing in (28) versus the model parameter *β*2. As we can see, the allowed range is within the vertical dashed lines.

β2

**Figure 1.** |∆*T<sup>f</sup>* /*T<sup>f</sup>* | vs. the model parameter *β*<sup>2</sup> (blue solid curve), for Model I: *f* = −*T* + *β*1 p *T*<sup>2</sup> + *β*2*TG*. The allowed range of *β*2, where (26) is satisfied (horizontal red dashed line), is within the vertical dashed lines.

*4.2. Model II: f* = −*T* + *a*1*T* <sup>2</sup> + *a*2*T* p |*TG*|

Let us now study the case *f* = −*T* + *a*1*T* <sup>2</sup> + *a*2*T* p |*TG*|, where *a*1, *a*<sup>2</sup> are the free parameters of the theory [28]. In this case, we find

$$\frac{\Delta T\_f}{T\_f} \quad = \frac{3}{10} c\_q^{-1} \zeta^3 T\_f \left\{ \frac{\Omega\_{DE0}}{3H\_0^2} - \sqrt{6}a\_2 \left[ \frac{\sqrt{H\_0^2 + \dot{H}\_0}}{6H\_0} \left( 6 - \frac{2\dot{H}\_0^2 - H\_0 \dot{H}\_0}{H\_0^2 + \dot{H}\_0} \right) - 1 \right] \right\}. \tag{32}$$

Using the constraint (26), and according to (32), ∆*T<sup>f</sup>* /*T<sup>f</sup>* is linear in *a*2; we deduce that (32) is valid for a small region around 2.7 <sup>×</sup> <sup>10</sup><sup>83</sup> GeV−<sup>2</sup> , where we have used the constraint from current cosmological era (27)

$$a\_1 \quad = \quad \frac{\Omega\_{\rm DE0}}{18H\_0^2} - \sqrt{6}a\_2 \frac{\sqrt{H\_0^2 + \dot{H}\_0}}{36H\_0} \left[ 6 - \frac{2\dot{H}\_0^2 - H\_0 \dot{H}\_0}{\left(H\_0^2 + \dot{H}\_0\right)^2} \right]. \tag{33}$$

Using the above value of *<sup>a</sup>*2, we find that *<sup>a</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup>1.1 <sup>×</sup> <sup>10</sup><sup>83</sup> GeV−<sup>2</sup> .

*4.3. Model III: f* =−*T* + *β*<sup>1</sup> p *T*2+*β*2*T<sup>G</sup>* + *a*1*T* <sup>2</sup> + *a*2*T* p |*TG*|

Now, we analyze the model *f* = −*T* + *β*<sup>1</sup> p *T*<sup>2</sup> + *β*2*T<sup>G</sup>* + *a*1*T* <sup>2</sup> + *a*2*T* p |*TG*|, where we have four free parameters, namely, *β*1, *β*2, *a*1, *a*<sup>2</sup> [28]. In order to simplify the analysis, we will impose the constraint <sup>−</sup>2.99 <sup>&</sup>lt; *<sup>β</sup>*<sup>2</sup> <sup>&</sup>lt; <sup>3</sup> 2 , obtained above.

In this case, we find

$$\begin{split} \frac{\Delta T\_f}{T\_f} &= -\left(60c\_qT\_f^3\right)^{-1} \Big[ 3\sqrt{12}\beta\_1 (3-2\beta\_2)^{-1/2} (1+\beta\_2-2\beta\_1) \\ &- 18\left[\frac{\Omega\_{DE0}}{3H\_0^2} + \frac{\sqrt{12}\beta\_1}{18H\_0^3} \Big[ (3+2\beta\_2)H\_0^2 + 2\beta\_2 H\_0 \Big]^{-1/2} \\ &\cdot \Big[ (3-6\beta\_1+2\beta\_2)H\_0^2 + 2\beta\_2 H\_0 \Big] \\ &- \frac{a\_2}{\sqrt{6}H\_0} \sqrt{H\_0^2 + H\_0} \Big[ 6 - \frac{2H\_0^2 - H\_0H\_0}{(H\_0^2 + H\_0)^2} \right] + \sqrt{6}a\_2 \\ &- \frac{\sqrt{12}\beta\_1\beta\_2}{18H\_0^3} \Big[ (3+2\beta\_2)H\_0^2 + 2\beta\_2 H\_0 \Big]^{-3/2} \\ &- \left[ (3+2\beta\_2)H\_0^4 + (9+8\beta\_2)H\_0^2 H\_0 \right] \\ &+ \beta\_2 \Big( 4\dot{H}\_0^2 + H\_0\dot{H}\_0 \Big) \Big] \Big\{ \xi^2 T\_f^4 \Big] \zeta. \end{split}$$

Observing that expression (34) is linear in *a*2, and using the constraint (26) and two values for *β*<sup>1</sup> from the aforementioned range we extracted in model I, i.e., *β*<sup>1</sup> = 1.4 and *<sup>β</sup>*<sup>2</sup> <sup>=</sup> 1, we find that (32) is valid for a small region around the point <sup>−</sup>3.5 <sup>×</sup> <sup>10</sup><sup>83</sup> GeV−<sup>2</sup> . Using another set of values (*β*<sup>1</sup> = 0.001, *β*<sup>2</sup> ≈ −2.96), we find that (32) is valid for a small region around the point <sup>−</sup>5.3 <sup>×</sup> <sup>10</sup><sup>83</sup> GeV−<sup>2</sup> , where we have used

$$\begin{split} a\_{1} = \frac{\Omega\_{\rm DE}}{18H\_{0}^{2}} &+ \frac{\sqrt{12}\beta\_{1}}{108H\_{0}^{3}} \Big[ (3+2\beta\_{2})H\_{0}^{2} + 2\beta\_{2}H\_{0} \Big]^{-1/2} \\ & \qquad \qquad \qquad \qquad \qquad \qquad \Big[ (3-6\beta\_{1}+2\beta\_{2})H\_{0}^{2} + 2\beta\_{2}H\_{0} \Big] \\ &- \frac{\sqrt{6}}{36} \frac{a\_{2}}{H\_{0}} \sqrt{H\_{0}^{2} + H\_{0}} \Big( 6 - \frac{2H\_{0}^{2} - H\_{0}H\_{0}}{\left(H\_{0}^{2} + H\_{0}\right)^{2}} \Big) \\ & \qquad \qquad \qquad \qquad - \frac{\sqrt{12}\beta\_{1}\beta\_{2}}{108H\_{0}^{3}} \Big[ (3+2\beta\_{2})H\_{0}^{2} + 2\beta\_{2}H\_{0} \Big]^{-3/2} \\ &\qquad \qquad \qquad \qquad \qquad \times \Big[ (3+2\beta\_{2})H\_{0}^{4} + (9+8\beta\_{2})H\_{0}^{2}H\_{0} + \beta\_{2} \Big(4H\_{0}^{2} + H\_{0}H\_{0}\Big) \Big], \end{split} \tag{34}$$

from (27). Imposing the above range of *<sup>a</sup>*2, we find that *<sup>a</sup>*<sup>1</sup> <sup>=</sup> 1.4 <sup>×</sup> <sup>10</sup><sup>83</sup> GeV−<sup>2</sup> for the first case and *<sup>a</sup>*<sup>1</sup> <sup>=</sup> 2.2 <sup>×</sup> <sup>10</sup><sup>83</sup> GeV−<sup>2</sup> for the second.

*4.4. Model IV: f* = −*T* + *β*<sup>1</sup> *T* <sup>2</sup> + *β*2*T<sup>G</sup> n*

As a next model, we consider the power-law model *f* = −*T* + *β*<sup>1</sup> *T* <sup>2</sup> + *β*2*T<sup>G</sup> n* , where the free parameters are *β*1, *β*2, *n*. In this model, we use values of *β*1, *β*<sup>2</sup> in order to constrain the power *n*. In this case, repeating the above steps, we find

$$\frac{\Delta T\_f}{T\_f} = \left(10c\_q\right)^{-1} \Omega\_{DE0} H\_0^{2(1-n)} \xi^{4n-1} T\_f^{8n-7} (3 - 2\beta\_2)^{n-2}$$

$$\cdot \left[ (3 + 2\beta\_2) H\_0^2 + 2\beta\_2 H\_0 \right]^{2-n} \left[ \left( 9 - 12\beta\_2 + 4\beta\_2^2 \right) \right.$$

$$\left. - 2n \left( 18 - 39\beta\_2 + 18\beta\_2^2 \right) + 16n^2 \beta\_2 (2\beta\_2 - 3) \right]$$

$$\cdot \left\{ \left( 9 + 12\beta\_2 + 4\beta\_2^2 \right) H\_0^4 + 4\beta\_2 (3 + 2\beta\_2) H\_0^2 H\_0 + 4\beta\_2^2 H\_0^2$$

$$\right. \left. - 2n \left[ \left( 18 + 15\beta\_2 + 2\beta\_2^2 \right) H\_0^4 + \beta\_2 (27 + 12\beta\_2) H\_0^2 H\_0 \right. \right.$$

$$\left. + 6\beta\_2^2 H\_0^2 + 2\beta\_2^2 H\_0 H\_0 \right] + 2n^2 \beta\_2 \left[ 4(3 + 2\beta\_2) H\_0^2 H\_0 \right. \tag{35}$$

$$\left. + 4\beta\_2 H\_0^2 + 2\beta\_2 H\_0 H\_0 \right] \right]^{-1}. \tag{36}$$

We use the constraint (26) and four values for *β*<sup>2</sup> from the range we extracted in model I above. For *β*<sup>2</sup> ≈ −2.9, we find that the constraint (26) is valid for *n* . 0.5. Similarly, using the value *β*<sup>2</sup> = −2, we find *n* . 0.47, while for *β*<sup>2</sup> = −1 we find *n* . 0.46. Finally, for *β*<sup>2</sup> = 1, we find *n* . 0.47. We mention that we have used the relation

$$\begin{split} \beta\_{1} = -6(12)^{-n} H\_{0}^{2(1-n)} \Omega\_{\rm DED} [(3+2\beta\_{2}) H\_{0}^{2} + 2 \beta\_{2} H\_{0}]^{2-n} \\ \quad \cdot \left\{ \left(9 + 12\beta\_{2} + 4\beta\_{2}^{2} \right) H\_{0}^{4} + 4\beta\_{2} (3 + 2\beta\_{2}) H\_{0}^{2} t H\_{0} + 4\beta\_{2}^{2} H\_{0}^{2} \\ \quad - 2n \left[ \left(18 + 15\beta\_{2} + 2\beta\_{2}^{2} \right) H\_{0}^{4} + \beta\_{2} (27 + 12\beta\_{2}) H\_{0}^{2} H\_{0} \right. \\ \left. + 6\beta\_{2}^{2} H\_{0}^{2} + 2\beta\_{2}^{2} H\_{0} \dot{H}\_{0} \right] + 2n^{2} \beta\_{2} \left[ 4(3 + 2\beta\_{2}) H\_{0}^{2} \dot{H}\_{0} \right. \\ \left. \left. + 4\beta\_{2} H\_{0}^{2} + 2\beta\_{2} H\_{0} \dot{H}\_{0} \right] \right\}^{-1}, \end{split} \tag{36}$$

which arises from (27).

Now, taking *β*<sup>2</sup> ≈ −2.9, *n* . 0.5 we find *β*<sup>1</sup> ∈ - <sup>−</sup>6.1 <sup>×</sup> <sup>10</sup>−82, 0.0007 GeV2(1−2*n*) . Similarly, for *β*<sup>2</sup> = −2, *n* . 0.47 we find *β*<sup>1</sup> ∈ - <sup>−</sup>3.5 <sup>×</sup> <sup>10</sup>−74, 5.9 <sup>×</sup> <sup>10</sup>−<sup>6</sup> GeV2(1−2*n*) , while using *β*<sup>2</sup> = −1, *n* . 0.46 we find *β*<sup>1</sup> ∈ - <sup>−</sup>4.4 <sup>×</sup> <sup>10</sup>−58, 1.2 <sup>×</sup> <sup>10</sup>−<sup>6</sup> GeV2(1−2*n*) . Finally, for *β*<sup>2</sup> = 1, *n* . 0.47 we find *β*<sup>1</sup> ∈ - <sup>−</sup>6.4 <sup>×</sup> <sup>10</sup>−<sup>8</sup> , 9.0 <sup>×</sup> <sup>10</sup>−<sup>6</sup> GeV2(1−2*n*) .

In order to provide the above results in a more transparent way, in Figure 2, we present |∆*T<sup>f</sup>* /*T<sup>f</sup>* | from (35) in terms of the model parameter *n*. As we observe, *n* needs to be *n* . 0.5 to pass the BBN constraint (26).

**Figure 2.** |∆*T<sup>f</sup>* /*T<sup>f</sup>* | vs. the model parameter *n* (blue solid curve), for Model IV: *f* = −*T* + *β*1 *T* <sup>2</sup> + *β*2*T<sup>G</sup> <sup>n</sup>* with *<sup>β</sup>*<sup>2</sup> ≈ −2.90, and the upper bound for <sup>|</sup>∆*T<sup>f</sup>* /*T<sup>f</sup>* | from (26) (red dashed line). As we observe, constraints from BBN require *n* . 0.5.

*4.5. Model V: f* = −*T* + *α* ln *β*<sup>1</sup> *T* <sup>2</sup> + *β*2*T<sup>G</sup> n*

In the last model we examine is the logarithmic one, characterized by *f* = −*T* + *α* ln *β*<sup>1</sup> *T* <sup>2</sup> + *β*2*T<sup>G</sup> n* , where *β*1, *β*2, *n* are the free parameters. Repeating the above analysis, we find

$$\begin{split} \frac{\Delta T\_f}{T\_f} &= \left( 10c\_6 \xi T\_f^7 \right)^{-1} H\_0^2 \Omega\_{D \to 0} \Big\{ \ln \beta\_1 + n \left[ \ln 12 \right. \\ &\left. + 4 \ln \left( \zeta T\_f^2 \right) + \ln (3 - 2 \beta\_2) \right. \\ &\left. - 2(3 - 2 \beta\_2)^{-2} \left( 18 - 39 \beta\_2 + 18 \beta\_2^2 \right) \right] \Big\} \\ &\cdot \left\{ \ln \beta\_1 + n \left\{ \ln 12 + 2 \ln (H\_0) + \ln \left[ (3 + 2 \beta\_2) H\_0^2 + 2 \beta\_2 H\_0 \right] \right. \\ &\left. - 2 \left[ (3 + 2 \beta\_2) H\_0^2 + 2 \beta\_2 H\_0 \right] \right]^{-2} \Big[ \left( 18 + 15 \beta\_2 + 2 \beta\_2^2 \right) H\_0^4 \\ &\left. + \beta\_2 (27 + 12 \beta\_2) H\_0^2 H\_0 + 6 \beta\_2^2 H\_0^2 + 2 \beta\_2^2 H\_0 H\_0 \right) \Big\} \right\}^{-1} . \tag{37} \end{split}$$

where using relation (27) we find

$$\begin{split} \mathfrak{a} = -6H\_0^2 \Omega\_{DE0} \Big\{ \ln \beta\_1 + n \Big[ \ln 12 + 2\ln(H\_0) \\ \qquad \qquad \qquad \qquad + \ln \Big[ (3 + 2\beta\_2)H\_0^2 + 2\beta\_2 H\_0 \Big] \Big] \\ \qquad \qquad \qquad - 2 \Big[ (3 + 2\beta\_2)H\_0^2 + 2\beta\_2 H\_0 \Big]^{-2} \Big[ (18 + 15\beta\_2 + 2\beta\_2^2)H\_0^4 \\ \qquad \qquad \qquad \qquad \qquad \qquad + \beta\_2 (27 + 12\beta\_2)H\_0^2 H\_0 + 6\beta\_2^2 H\_0^2 + 2\beta\_2^2 H\_0 H\_0 \Big] \Big\} \Big]^{-1} . \end{split} \tag{38}$$

We consider the values *β*<sup>1</sup> = 0.001 GeV −4*n* , *β*<sup>2</sup> ≈ −2.9, and we find that *n* is allowed to take every value apart from −0.0003 and a very small region around it since (37) diverges. Moreover, *α* is allowed to take every value apart from 0, which is the value it obtains using the above narrow window for *n*. Using the same considerations as the above models, we find that for *β*<sup>1</sup> = 0.001 GeV −4*n* , *β*<sup>2</sup> = −2 the value of *n* is allowed to take every value apart from −0.012 and *a* every value but 0. Similarly, for *β*<sup>1</sup> = 0.001 GeV −4*n* , *β*<sup>2</sup> = −1 we find that *n* 6= −0.018 and *a* 6= 0, while for *β*<sup>1</sup> = 0.001 GeV −4*n* , *β*<sup>2</sup> = 1 we find *n* 6= −0.018 and *a* 6= 0.

As an example, in Figure 3 we present |∆*T<sup>f</sup>* /*T<sup>f</sup>* | from (37) as a function of the model parameter *n*. The model parameter *n* is allowed to take all possible values except those values around a very small region centered at −0.0003, in which (37) diverges. Hence, we conclude that the logarithmic *f*(*T*, *TG*) model can easily satisfy the BBN bounds.

n

**Figure 3.** |∆*T<sup>f</sup>* /*T<sup>f</sup>* | vs. the model parameter *n* (blue solid curve), for Model V: *f* = −*T* + *α* ln *β*<sup>1</sup> *T* <sup>2</sup> + *β*2*T<sup>G</sup> n* , choosing *β*<sup>1</sup> = 0.001 GeV −4*n* , *β*<sup>2</sup> ≈ −2.7. The vertical dashed line at *n* − 0.0003 denotes the point where (37) diverges.

### **5. Conclusions**

Modified gravity aims to provide explanations for various epochs of the universe evolution, and at the same time to improve the renormalizability issues of general relativity. Nevertheless, despite the specific advantages at a given era of cosmological evolution, one should be very careful not to spoil the other, well understood and significantly constrained, phases, such as the big bang nucleosynthesis (BBN) one.

In particular, there are many modified gravity models, which are constructed phenomenologically in order to be able to describe the late-time universe evolution at both the background and perturbation level. Typically, these models are confronted with observational data such as Supernovae Type Ia (SNIa), Baryonic Acoustic Oscillations (BAO), cosmic microwave background (CMB), cosmic chronometers (CC), gamma-ray bursts (GRB), growth data, etc. The problem is that although modified gravity scenarios, through the extra terms they induce, are very efficient in describing the late-time universe, quite often they induce significant terms at early times too, thus spoiling the early-time evolution, such as the BBN phase, in which the concordance cosmological paradigm is very successful. Hence, independently of the late-universe successes that a modified gravity model may have, one should always examine whether the model can pass the BBN constraints too.

In the present work, we confronted one interesting class of gravitational modification, namely, *f*(*T*, *TG*) gravity, with BBN requirements. The former is obtained using both the torsion scalar, as well as the teleparallel equivalent of the Gauss–Bonnet term, in the Lagrangian. Hence, one obtains modified Friedmann equations in which the extra torsional terms constitute an effective dark energy sector.

We started by calculating the deviations of the freeze-out temperature *T<sup>f</sup>* , caused by the extra torsion terms, in comparison to ΛCDM paradigm. We imposed five specific *f*(*T*, *TG*) models that have been proposed in the literature in phenomenological grounds, i.e., in order to be able to describe the late-time evolution and lead to acceleration without an explicit cosmological constant. Hence, we extracted the constraints on the model parameters in order for the ratio |∆*T<sup>f</sup>* /*T<sup>f</sup>* | to satisfy the BBN bound ∆*T<sup>f</sup> Tf* <sup>&</sup>lt; 4.7 <sup>×</sup> <sup>10</sup> −4 . As we found, in most of the models the involved parameters are bounded in a narrow window around their general relativity values, as expected. However, the logarithmic model can easily satisfy the BBN constraints for large regions of the model parameters, which acts as an advantage for this scenario.

We stress here that we did not fix the cosmological parameters to their general relativity values; on the contrary, we left them completely free and we examined which parameter regions are allowed if we want the models to pass the BBN constraints. The fact that in most models the parameter regions are constrained to a narrow window around their general relativity values was in some sense expected, but in general is not guaranteed or known a priori since many modified gravity models are completely excluded under the BBN analysis since for all parameter regions their early-universe effect is huge.

In conclusion, *f*(*T*, *TG*) gravity, apart from having interesting cosmological implications both in the inflationary and late-time phase, possesses particular sub-classes that can safely pass BBN bounds; nevertheless, the torsional modification is constrained in narrow windows around the general relativity values. This feature should be taken into account in future model building.

**Author Contributions:** Conceptualization, E.N.S.; methodology, E.N.S., P.A.; software, P.A.; validation, E.N.S., S.B.; formal analysis, P.A.; investigation, P.A.; data curation, P.A.; writing—original draft preparation, E.N.S.; writing—review and editing, P.A., S.B., K.Y.; visualization, P.A.; supervision, E.N.S., S.B.; project administration, P.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This research is co-financed by Greece and the European Union (European Social Fund-ESF) through the Operational Programme "Human Resources Development, Education and Lifelong Learning" in the context of the project "Strengthening Human Resources Research Potential via Doctorate Research" (MIS-5000432), implemented by the State Scholarships Foundation (IKY). The work of N.E.M is supported in part by the UK Science and Technology Facilities research Council (STFC) under the research grant ST/T000759/1. S.B., N.E.M., and E.N.S. also acknowledge participation in the COST Association Action CA18108 "Quantum Gravity Phenomenology in the Multimessenger Approach (QG-MM)".

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Universe* Editorial Office E-mail: universe@mdpi.com www.mdpi.com/journal/universe

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-7982-5