**2.** *f* (*T***,** *TG*) **Gravity**

In this section, we briefly review *f*(*T*, *TG*) gravity [26–28]. As usual in torsional formulation of gravity, we use the tetrad field as the dynamical variable, which forms an orthonormal basis at the tangent space. In a coordinate basis, one can relate it with the metric through *gµν*(*x*) = *ηABe A µ* (*x*)*e B ν* (*x*), where *ηAB* = diag(−1, 1, 1, 1), and with Greek and Latin letters, denoting coordinate and tangent indices, respectively. Applying the Weitzenböck connection *W<sup>λ</sup> νµ* ≡ *e λ A ∂µe A ν* [25], the corresponding torsion tensor is

$$T^{\lambda}\_{\mu\nu} \equiv \mathcal{W}^{\lambda}\_{\nu\mu} - \mathcal{W}^{\lambda}\_{\mu\nu} = e^{\lambda}\_{A} \left( \partial\_{\mu} e^{A}\_{\nu} - \partial\_{\nu} e^{A}\_{\mu} \right) \tag{1}$$

and then the torsion scalar is obtained through the contractions

$$T \equiv \frac{1}{4} T^{\rho\mu\nu} T\_{\rho\mu\nu} + \frac{1}{2} T^{\rho\mu\nu} T\_{\nu\mu\rho} - T\_{\rho\mu}{}^{\rho} T\_{\nu\nu}^{\nu\mu} \tag{2}$$

and incorporates all information of the gravitational field. Used as a Lagrangian, the torsion scalar gives rise to exactly the same equations with general relativity, which is why the theory was named the teleparallel equivalent of general relativity (TEGR).

Similarly to curvature gravity, where one can construct higher-order invariants such as the Gauss–Bonnet one, in torsional gravity one may construct higher-order torsional invariants, too. In particular, since the curvature (Ricci) scalar and the torsion scalar differ by a total derivative, in [26] the authors followed the same recipe and extracted a higher-order torsional invariant that differs from the Gauss–Bonnet one by a boundary term, namely

$$\begin{split} T\_{\rm G} &= \left( \mathcal{K}^{\kappa}{}\_{\rho\tau} \mathcal{K}^{\theta\lambda}{}\_{\rho} \mathcal{K}^{\mu}{}\_{\chi\sigma} \mathcal{K}^{\chi\nu}{}\_{\tau} - 2 \mathcal{K}^{\kappa\lambda}{}\_{\tau} \mathcal{K}^{\mu}{}\_{\rho\rho} \mathcal{K}^{\rho}{}\_{\chi\sigma} \mathcal{K}^{\chi\nu}{}\_{\tau} \right. \\ &+ 2 \mathcal{K}^{\kappa\lambda}{}\_{\pi} \mathcal{K}^{\mu}{}\_{\rho\rho} \mathcal{K}^{\rho\nu}{}\_{\chi\sigma} \mathcal{K}^{\chi}{}\_{\sigma\tau} + 2 \mathcal{K}^{\kappa\lambda}{}\_{\pi} \mathcal{K}^{\mu}{}\_{\rho\rho} \mathcal{K}^{\rho\nu}{}\_{\sigma\tau} \right) \delta^{\pi\rho\tau\tau}\_{\kappa\lambda\mu\nu} \tag{3} \end{split}$$

where *K µν <sup>ρ</sup>* ≡ −<sup>1</sup> 2 *T µν <sup>ρ</sup>* − *T νµ <sup>ρ</sup>* − *T µν ρ* is the contortion tensor, and the generalized *δ πρστ κλµν* denotes the determinant of the Kronecker deltas. Note that similarly to the Gauss–Bonnet term, the teleparallel equivalent of the Gauss–Bonnet term *T<sup>G</sup>* is also a topological invariant in four dimensions.

Using the above torsional invariants, one can construct the new class of *f*(*T*, *TG*) gravitational modifications, characterized by the action [26]

$$S = \frac{M\_P^2}{2} \int d^4 \mathbf{x} \, e \, f(T\_\prime \, T\_G) \, \, \, \, \tag{4}$$

with *M*<sup>2</sup> *P* the reduced Planck mass. The general field equations of the above action can be found in [26], where one can clearly see that the theory is different from *f*(*R*), *f*(*R*, *G*), and *f*(*T*) gravitational modifications, and thus it corresponds to a novel class of modified gravity.

In this work, we are interested in the cosmological applications of *f*(*T*, *TG*) gravity. Hence, we consider a spatially flat Friedmann–Robertson–Walker (FRW) metric of the form

$$ds^2 = -dt^2 + a^2(t)\delta\_{\dot{l}\dot{\jmath}}d\mathbf{x}^{\dot{l}}d\mathbf{x}^{\dot{\jmath}}.\tag{5}$$

with *a*(*t*) the scale factor, which corresponds to the diagonal tetrad

$$e^A\_{\;\;\mu} = \text{diag}(1, a(t), a(t), a(t)). \tag{6}$$

In this case, the torsion scalar (2) and the teleparallel equivalent of the Gauss–Bonnet term (3) become

$$T = 6H^2\tag{7}$$

$$T\_{\mathbb{G}} = 24H^2 \left(\dot{H} + H^2\right),\tag{8}$$

with *H* = *<sup>a</sup>*˙ *a* the Hubble parameter and where dots denote derivatives with respect to *t*. The general field equations for the FRW geometry are [27]

$$\begin{aligned} f - 12H^2 f\_T - T\_{\mathcal{G}} f\_{T\_{\mathcal{G}}} + 24H^3 f\_{T\_{\mathcal{G}}} &= 2M\_P^{-2} (\rho\_r + \rho\_m) \\ f - 4(3H^2 + \dot{H}) f\_T - 4H \dot{f}\_T - T\_{\mathcal{G}} f\_{T\_{\mathcal{G}}} \end{aligned} \tag{9}$$

$$+\frac{2}{3H}T\_{\!G}f\_{\!\!T\_{\!G}}^{\cdots} + 8H^2 f\_{\!\!T\_{\!G}}^{\cdots} = -2M\_{\!\!P}^{-2}(p\_r + p\_m) \tag{10}$$

with ˙ *<sup>f</sup><sup>T</sup>* <sup>=</sup> *<sup>f</sup>TTT*˙ <sup>+</sup> *<sup>f</sup>TT<sup>G</sup> T*˙ *G*, ˙ *fT<sup>G</sup>* = *fTT<sup>G</sup> <sup>T</sup>*˙ <sup>+</sup> *<sup>f</sup>TGT<sup>G</sup> T*˙ *<sup>G</sup>*, and ¨ *fT<sup>G</sup>* = *fTTT<sup>G</sup> <sup>T</sup>*˙ <sup>2</sup> <sup>+</sup> <sup>2</sup> *<sup>f</sup>TTGT<sup>G</sup> T*˙ *T*˙ *<sup>G</sup>* + *fTGTGT<sup>G</sup> T*˙ <sup>2</sup> *<sup>G</sup>* + *fTT<sup>G</sup> <sup>T</sup>*¨ <sup>+</sup> *<sup>f</sup>TGT<sup>G</sup> T*¨ *<sup>G</sup>*, and where *fTT*, *fTT<sup>G</sup>* ,... denote multiple partial differentiations with respect to *T* and *TG*. Note that in the above equations, we have also introduced the radiation and matter sectors, corresponding to perfect fluids with energy densities *ρr* , *ρ<sup>m</sup>* and pressures *p<sup>r</sup>* , *pm*, respectively. Lastly, we mention that the above equations for *<sup>f</sup>*(*T*, *<sup>T</sup>G*) = −*<sup>T</sup>* + <sup>Λ</sup> recover the TEGR and general relativity equations, where <sup>Λ</sup> is the cosmological constant.

As we can see, we can re-write the Friedmann Equations (9) and (10) in the usual form

$$\Im \mathcal{M}\_{\rm P}^{2} H^{2} = (\rho\_{\rm r} + \rho\_{\rm m} + \rho\_{\rm DE}) \tag{11}$$

$$-2M\_P^2 \dot{H} = (\rho\_r + p\_r + \rho\_m + p\_m + \rho\_{DE} + p\_{DE})\_\prime \tag{12}$$

where we have defined the effective dark energy density and pressure as

$$\rho\_{DE} \equiv \frac{M\_P^2}{2} \Big( 6H^2 - f + 12H^2 f\_T + T\_G f\_{T\_G} - 24H^3 f\_{T\_G}^\cdot \Big),\tag{13}$$

$$p\_{DE} \equiv \frac{M\_P^2}{2} \left[ -2(2H + 3H^2) + f - 4(H + 3H^2)f\_T \right]$$

$$-4H\dot{f}\_T - T\_G f\_{T\_G} + \frac{2}{3H} T\_G \dot{f}\_{T\_G} + 8H^2 \ddot{f}\_{T\_G} \right] \tag{14}$$

of gravitational origin.

### **3. Big Bang Nucleosynthesis Constraints**

Big bang nucleosynthesis (BBN) was a process that took place during radiation era. Let us first present the framework, which provides the BBN constraints through standard

cosmology [76–80]. The first Friedmann equation from Einstein–Hilbert action can be written as

$$\text{CH}^2 = M\_P^{-2} \rho\_\prime \tag{15}$$

where *ρ* = *ρ<sup>r</sup>* + *ρm*. In the radiation era, the radiation sector dominates; hence, we can write

$$H^2 \approx \frac{M\_P^{-2}}{3} \rho\_r \equiv H\_{\rm GR}^2. \tag{16}$$

In addition, it is known that the energy density of relativistic particles is

$$
\rho\_r = \frac{\pi^2}{30} \text{g}\_\* T^4 \text{.}\tag{17}
$$

where *g*<sup>∗</sup> ∼ 10 is the effective number of degrees of freedom and *T* is the temperature. Thus, if we combine (16) with (17) we obtain

*H*(*T*) ≈ 4*π* 3*g*∗ <sup>45</sup> 1/2 *<sup>T</sup>* 2 *MPl* , (18)

where *MPl* = (8*π*) 1 <sup>2</sup> *<sup>M</sup><sup>P</sup>* <sup>=</sup> 1.22 <sup>×</sup> <sup>10</sup><sup>19</sup> GeV is the Planck mass.

During the radiation era, the scale factor evolves as *a*(*t*) ∼ *t* 1/2. Therefore, using the relation of the Hubble parameter with the scale factor, we find that in the radiation era the Hubble parameter evolves as *H*(*t*) = <sup>1</sup> 2*t* . Combining the last one with (18), we find the relation between temperature and time. Thus, we have <sup>1</sup> *t* ≃ 32*π* 3*g*∗ <sup>90</sup> 1/2 *<sup>T</sup>* 2 *MPl* (or

$$T(t) \underset{-}{\simeq} (t/\text{sec})^{-1/2} \underset{-}{\text{MeV}}.$$

During the BBN, we have interactions between particles. For example, we have interactions between neutrons, protons, electrons, and neutrinos, namely, *n* + *ν<sup>e</sup>* → *p* + *e* −, *n* + *e* <sup>+</sup> <sup>→</sup> *<sup>p</sup>* <sup>+</sup> *<sup>ν</sup>*¯*<sup>e</sup>* , and *n* → *p* + *e* <sup>−</sup> + *ν*¯*<sup>e</sup>* . We name the conversion rate from a particle A to particle B as *λBA*. Hence, the conversion rate from neutrons to protons is *λpn*, and it is equal to the sum of the three interaction conversion rates written above. Therefore, the calculation of the neutron abundance arises from the protons-neutron conversion rate [78,79]

$$
\lambda\_{pn}(T) = \lambda\_{(n+\nu\_{\ell}\to p+\epsilon^{-})} + \lambda\_{(n+\epsilon^{+}\to p+\tilde{\nu}\_{\ell})} + \lambda\_{(n\to p+\epsilon^{-}+\tilde{\nu}\_{\ell})} \tag{19}
$$

and its inverse *λnp*(*T*), and therefore for the total rate we have *λtot*(*T*) = *λnp*(*T*) + *λpn*(*T*). Now, we assume that the various particle (neutrino, electron, and photon) temperatures are the same and low enough in order to use the Boltzmann distribution instead of the Fermi-Dirac one, and we neglect the electron mass compared to the electron and neutrino energies. The final expression for the conversion rate is [81–84]

$$
\lambda\_{tot}(T) = 4A \, T^3 (4!T^2 + 2 \times 3!QT + 2!Q^2) \,. \tag{20}
$$

where *<sup>Q</sup>* <sup>=</sup> *<sup>m</sup><sup>n</sup>* <sup>−</sup> *<sup>m</sup><sup>p</sup>* <sup>=</sup> 1.29 <sup>×</sup> <sup>10</sup>−<sup>3</sup> GeV is the mass difference between neutron and proton and *<sup>A</sup>* <sup>=</sup> 1.02 <sup>×</sup> <sup>10</sup>−<sup>11</sup> GeV−<sup>4</sup> .

We proceed in calculating the corresponding freeze-out temperature. This will arise comparing the universe expansion rate <sup>1</sup> *<sup>H</sup>* with *<sup>λ</sup>tot*(*T*). In particular, if <sup>1</sup> *<sup>H</sup>* ≪ *λtot*(*T*), namely, if the expansion time is much smaller than the interaction time, we can consider thermal equilibrium [76,77]. On the contrary, if <sup>1</sup> *<sup>H</sup>* ≫ *λtot*(*T*) then particles do not have enough time to interact so they decouple. The freeze-out temperature *T<sup>f</sup>* , in which the decoupling takes place, corresponds to *H*(*T<sup>f</sup>* ) = *<sup>λ</sup>tot Tf* ≃ *c<sup>q</sup> T* 5 *f* , with *c<sup>q</sup>* ≡ 4*A* 4! ≃ 9.8 <sup>×</sup> <sup>10</sup>−<sup>10</sup> GeV−<sup>4</sup> [81–84]. Now, if we use (18) and *H*(*T<sup>f</sup>* ) = *<sup>λ</sup>tot Tf* ≃ *c<sup>q</sup> T* 5 *f* , we acquire

$$T\_f = \left(\frac{4\pi^3 \text{g}\_\*}{45M\_{Pl}^2 c\_q^2}\right)^{1/6} \sim 0.0006\text{ GeV}.\tag{21}$$

Using modified theories, we obtain extra terms in energy density due to the modification of gravity. The first Friedmann Equation (11) during radiation era becomes

$$3\mathcal{M}\_\mathrm{P}^2 H^2 = \rho\_r + \rho\_{\mathrm{DE}} \tag{22}$$

where *ρDE* must be very small compared to *ρ<sup>r</sup>* in order to be in accordance with observations. Hence, we can write (22) using (16) as

$$H = H\_{GR}\sqrt{1 + \frac{\rho\_{DE}}{\rho\_r}} = H\_{GR} + \delta H\_\prime \tag{23}$$

where *HGR* is the Hubble parameter of standard cosmology. Thus, we have ∆*H* = q 1 + *ρDE ρr* − 1 *HGR*, which quantifies the deviation from standard cosmology, i.e., form *HGR*. This will lead to a deviation in the freeze-out temperature ∆*T<sup>f</sup>* . Since *HGR* = *λtot* ≈ *c<sup>q</sup> T* 5 *f* and <sup>q</sup> 1 + *ρDE ρr* <sup>≈</sup> <sup>1</sup> <sup>+</sup> <sup>1</sup> 2 *ρDE ρr* , we easily find

$$\left(\sqrt{1+\frac{\rho\_{DE}}{\rho\_r}}-1\right)H\_{GR} = 5c\_q \ T\_f^4 \Delta T\_{f\prime} \tag{24}$$

and finally

$$\frac{\Delta T\_f}{T\_f} \simeq \frac{\rho\_{DE}}{\rho\_r} \frac{H\_{GR}}{10 c\_q \ T\_f^5} \,\mathrm{}\,\mathrm{}\,\tag{25}$$

where we used that *<sup>ρ</sup>DE* << *<sup>ρ</sup><sup>r</sup>* during BBN era. This theoretically calculated <sup>∆</sup>*T<sup>f</sup> Tf* should be compared with the observational bound

$$
\left|\frac{\Delta T\_f}{T\_f}\right| < 4.7 \times 10^{-4} \,\text{.}\tag{26}
$$

which is obtained from the observational estimations of the baryon mass fraction converted to <sup>4</sup>*He* [85–91].
