Model-independent Reconstruction of f(T) Gravity from Gaussian Processes

Modified gravity is an effective approach to describe the acceleration of the early- and late-time universe (Capozziello & Laurentis 2011a; Nojiri & Odintsov 2011), except for the introduction of inflation and/or dark energy components (Peebles & Ratra 2003; Cai et al. 2010; Elizalde et al. 2018; Li & Shafieloo 2019; Li et al. 2019). In particular, interest in studying modified gravity has been curbed by the recently reported measurements of the H₀ tension that has failed to be addressed within the standard ΛCDM cosmology (Freedman et al. 2019; Vagnozzi 2019; Wong et al. 2019). Among various constructions of modified gravity one can find the interesting class that is based on the torsional formulation (for a review see Cai et al. 2016). In particular, starting from the simplest torsional gravity, namely the teleparallel equivalent of general relativity (Einstein 1928; Hayashi & Shirafuji 1979, 1982; Aldrovandi & Pereira 2013; Maluf 2013), one can construct modifications such as the f(T) gravity (Bengochea & Ferraro 2009; Linder 2010; Bamba et al. 2011, 2013, 2016; Cai et al. 2011, 2016, 2018; Capozziello et al. 2011; Chen et al. 2011; Dent et al. 2011; Myrzakulov 2012; Wei et al. 2012; Wu & Geng 2012; Amoros et al. 2013; Li et al. 2013, 2018; Ong et al. 2013; Otalora 2013; Haro & Amoros 2014; Nashed & El Hanafy 2014; Darabi et al. 2015; El Hanafy & Nashed 2015; Farrugia & Said 2016; Guo et al. 2016; Bahamonde et al. 2017; Malekjani et al. 2017; Qi et al. 2017; Iosifidis & Koivisto 2018; Karpathopoulos et al. 2018; Krššák et al. 2019; Nunes 2018), the f(T, T_G) gravity (Kofinas & Saridakis 2014a, 2014b), scalar-torsion theories (Geng et al. 2011; Hohmann et al. 2018), etc. Finally, one can proceed beyond the teleparallel framework and construct more complicated torsional theories.

One key question for any theory of modified gravity is how to determine a viable choice from the arbitrary functions that have been involved. Some general features may be determined by theoretical arguments such as the theoretical requirements for a ghost-free theory that possesses stable perturbations etc., or the desire for the action to possess Noether symmetries; however, the main tool for constraining the possible forms of modification remains that of confrontation with observations. The general recipe is to consider by hand a variety of specific forms inside some general class, apply them in a cosmological framework, predict the dynamical behaviors at both the background and perturbation levels, and then use observational data to constrain the involved parameters or exclude the examined form (a similar procedure can also be followed to confront with local/solar system data). For the case of torsional gravity, cosmological confrontations have been performed in Wu & Yu (2010), Nesseris et al. (2013), Capozziello et al. (2015), Basilakos (2016), Nunes et al. (2016, 2017), Nunes (2018), Basilakos et al. (2018), and Anagnostopoulos et al. (2019), and the solar system tests can be found in Iorio & Saridakis (2012), Iorio et al. (2015), and Farrugia et al. (2016), while the latest limit from galaxy lensing has been presented in Chen et al. (2019). Hence, in the literature there exist at least three viable scenarios for f(T) gravity (Nesseris et al. 2013; Basilakos et al. 2018).

Although the above procedure of observational constraints is very useful in offering crucial information on the possible forms of modification, it is even more productive if the data can reconstruct the involved modification functions in a model-independent way without inserting an initial guess. Such a procedure has been successfully developed in the early-time inflationary cosmology (Copeland et al. 1993; Lidsey et al. 1997; Herrera 2018); however, concerning the late-time cosmology the involved complications allow only for a partial application, such as in the cosmography framework (Bamba et al. 2012; Capozziello et al. 2019) or backscattering procedure (Capozziello et al. 2017). Interestingly, a useful tool in the above reconstruction is the analysis of Gaussian processes (GPs) (Rasmussen & Williams 2005; Holsclaw et al. 2010; Seikel & Clarkson 2013), which allow one to investigate features of the form of the involved unknown functions in a model-independent way using only the given data sets. Such a procedure has been applied to dark energy models using various data sets in order to reconstruct the evolution of the Hubble function, of the dark energy equation-of-state (EoS) parameters, of the dimensionless comoving luminosity distance, of the dark interaction term, etc. (Seikel et al. 2012; Shafieloo et al. 2012; Kim et al. 2013; Seikel & Clarkson 2013; Nair et al. 2014; Santos-da-Costa et al. 2015; Yang et al. 2015; Zhang & Xia 2016; Wang & Meng 2017; Elizalde et al. 2018; Gomez-Valent & Amendola 2018; Melia & Yennapureddy 2018; Pinho et al. 2018; Zhang & Li 2018; Elizalde & Khurshudyan 2019; Gomez-Valent & Amendola 2019; Rau et al. 2019; Yin & Wei 2019).

In this article we develop the GP analysis for the case of f(T) cosmology, in order to reconstruct the form of the f(T) modification in a model-independent way, namely using as the only input the observational data sets of Hubble function measurements H(z). Such a procedure becomes easy in the case of f(T) cosmology, since the latter has the advantage that the torsion scalar is a simple function of H, namely T = −H², and thus eventually all cosmological equations can be expressed in terms of H(z) and its derivative. Hence, reconstructing H(z) and its derivative from the Hubble data through the GP analysis, leads to the reconstruction of the f(T) form itself, without any assumption.

The plan of the article is as follows. In Section 2 we provide a brief review of the cosmological equations of f(T) gravity. In Section 3 we describe the basic ingredients of GPs and then we apply them to observational data on the Hubble function to reconstruct H(z) and its derivative. Then in Section 4 we use this reconstruction to reconstruct the form of f(T) in a model-independent way, and to extract a constraint on various f(T) models in the literature. Finally, in Section 5 we summarize our results with a discussion.

In this section we briefly review f(T) gravity and cosmology. In the torsional formulation one uses as dynamical variables the vierbeins fields, which form an orthonormal basis at a manifold point. In a coordinate basis they are related to the metric through ${g}_{\mu \nu }(x)={\eta }_{{AB}}\,{e}_{\mu }^{A}(x)\,{e}_{\nu }^{B}(x)$, with Greek and Latin indices respectively used for the coordinate and tangent space. In the particular class of teleparallel gravity one introduces the Weitzenböck connection ${W}_{\nu \mu }^{\lambda }\equiv {e}_{A}^{\lambda }\ {\partial }_{\mu }{e}_{\nu }^{A}$ (Weitzenböck 1923), and thus the corresponding torsion tensor is

$$\begin{eqnarray}&&{T}_{\,\mu \nu }^{\lambda }\equiv {W}_{\nu \mu }^{\lambda }-{W}_{\mu \nu }^{\lambda }={e}_{A}^{\lambda }\ ({\partial }_{\mu }{e}_{\nu }^{A}-{\partial }_{\nu }{e}_{\mu }^{A}),\end{eqnarray} \tag{ 1 }$$

while the corresponding curvature tensor is zero. We mention here that the above formulation is constructed in a specific cosmological gauge where the spin connection components are vanishing, which simplifies the analysis and allows us to impose the specific vierbein choice below. Leaving a non-zero spin connection is also possible, but this would involve more complicated vierbeins (Krššák & Saridakis 2016; Hohmann et al. 2018).

The torsion tensor incorporates all the information on the gravitational field, and the torsion scalar arises from its contraction as

$$\begin{eqnarray}&&T\equiv \displaystyle \frac{1}{4}{T}^{\rho \mu \nu }{T}_{\rho \mu \nu }+\displaystyle \frac{1}{2}{T}^{\rho \mu \nu }{T}_{\nu \mu \rho }-{T}_{\rho \mu }{}^{\rho }{T}^{\nu \mu }{}_{\nu }.\end{eqnarray} \tag{ 2 }$$

This forms the Lagrangian of teleparallel gravity (similarly to general relativity where the Lagrangian is the Ricci scalar), and since variation in terms of the vierbeins gives the same equations as general relativity, the constructed theory was named teleparallel equivalent of general relativity (TEGR).

One can start from the TEGR to derive the torsional-based modifications. The simplest extension is to generalize T in the action to be T + f(T) (Cai et al. 2016):

$$\begin{eqnarray}&&S=\displaystyle \frac{1}{16\pi G}\int {d}^{4}x\,e\left[T+f(T)+{L}_{m}\right],\end{eqnarray} \tag{ 3 }$$

where $e=\det ({e}_{\mu }^{A})=\sqrt{-g}$, G is the gravitational constant, and where for completeness we have also included the matter Lagrangian L_m. Variation of the above action results in the field equations

$$\begin{eqnarray}&&\begin{array}{l}{e}^{-1}{\partial }_{\mu }({{ee}}_{A}^{\rho }{S}_{\rho }{}^{\mu \nu })[1+{f}_{T}]+{e}_{A}^{\rho }{S}_{\rho }{}^{\mu \nu }{\partial }_{\mu }(T){f}_{{TT}}\\ \qquad -\,[1+{f}_{T}]{e}_{A}^{\lambda }{T}^{\rho }{}_{\mu \lambda }{S}_{\rho }{}^{\nu \mu }+\displaystyle \frac{1}{4}{e}_{A}^{\nu }[T+f(T)]\\ \quad =\,4\pi {{Ge}}_{A}^{\rho }{\mathop{T}\limits^{{\boldsymbol{e}}{\boldsymbol{m}}}}_{\rho }{}^{\nu },\end{array}\end{eqnarray} \tag{ 4 }$$

with ${f}_{T}\equiv \partial f/\partial T$, ${f}_{{TT}}\equiv {\partial }^{2}f/\partial {T}^{2}$, and where ${\mathop{T}\limits^{{\boldsymbol{e}}{\boldsymbol{m}}}}_{\rho }{}^{\nu }$ denotes the total matter (namely dark and baryonic matter) energy–momentum tensor. Finally, ${S}_{\rho }{}^{\mu \nu }\equiv {\textstyle \tfrac{1}{2}}\left({K}^{\mu \nu }{}_{\rho }+{\delta }_{\rho }^{\mu }\,{T}^{\alpha \nu }{}_{\alpha }-{\delta }_{\rho }^{\nu }\,{T}^{\alpha \mu }{}_{\alpha }\right)$ is the superpotential, with ${K}^{\mu \nu }{}_{\rho }\equiv -\tfrac{1}{2}\left({T}^{\mu \nu }{}_{\rho }-{T}^{\nu \mu }{}_{\rho }-{T}_{\rho }{}^{\mu \nu }\right)$ the contorsion tensor.

In order to apply f(T) gravity in a cosmological framework we impose the spatially flat Friedmann–Lemaître–Robertson–Walker metric

$$\begin{eqnarray}&&{{ds}}^{2}={{dt}}^{2}-{a}^{2}(t){\delta }_{{ij}}{{dx}}^{i}{{dx}}^{j},\end{eqnarray} \tag{ 5 }$$

which corresponds to the vierbein form ${e}_{\mu }^{A}=\mathrm{diag}(1,a,a,a)$, where a(t) is the scale factor. Inserting this choice into the general field Equation (4) we obtain the Friedmann equations of f(T) cosmology, namely

$$\begin{eqnarray}&&{H}^{2}=\displaystyle \frac{8\pi G}{3}{\rho }_{m}-\displaystyle \frac{f}{6}+\displaystyle \frac{{{Tf}}_{T}}{3},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\dot{H}=-\displaystyle \frac{4\pi G({\rho }_{m}+{P}_{m})}{1+{f}_{T}+2{{Tf}}_{{TT}}},\end{eqnarray} \tag{ 7 }$$

with $H\equiv \dot{a}/a$ the Hubble function and where dots denoting derivatives with respect to t. Additionally, in the above equations ρ_m and P_m are respectively the energy density and pressure of the matter fluid. Note that the torsion scalar (2) in Friedmann–Robertson–Walker (FRW) geometry becomes

$$\begin{eqnarray}&&T=-6{H}^{2},\end{eqnarray} \tag{ 8 }$$

and this expression proves very useful for the purpose of the present work.

As a next step we define an effective dark energy sector with energy density and pressure respectively given by

$$\begin{eqnarray}&&{\rho }_{{\rm{DE}}}\equiv \displaystyle \frac{3}{8\pi G}\left[-\displaystyle \frac{f}{6}+\displaystyle \frac{{{Tf}}_{T}}{3}\right],\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{P}_{{\rm{DE}}}\equiv \displaystyle \frac{1}{16\pi G}\left[\displaystyle \frac{f-{f}_{T}T+2{T}^{2}{f}_{{TT}}}{1+{f}_{T}+2{{Tf}}_{{TT}}}\right],\end{eqnarray} \tag{ 10 }$$

and therefore its EoS parameter reads

$$\begin{eqnarray}&&w\equiv \displaystyle \frac{{P}_{{\rm{DE}}}}{{\rho }_{{\rm{DE}}}}=-\displaystyle \frac{f/T-{f}_{T}+2{{Tf}}_{{TT}}}{\left[1+{f}_{T}+2{{Tf}}_{{TT}}\right]\left[f/T-2{f}_{T}\right]}.\end{eqnarray} \tag{ 11 }$$

Hence, the first Friedmann Equation (6) effectively acquires the standard form ${H}^{2}=\tfrac{8\pi G}{3}({\rho }_{m}+{\rho }_{{\rm{DE}}})$. Finally, the system of cosmological equations closes by considering the conservation of matter,

$$\begin{eqnarray}&&{\dot{\rho }}_{m}+3H({\rho }_{m}+{P}_{m})=0,\end{eqnarray} \tag{ 12 }$$

which using (6) and (7) implies additionally the conservation of the effective dark energy,

$$\begin{eqnarray}&&{\dot{\rho }}_{{\rm{DE}}}+3H({\rho }_{{\rm{DE}}}+{P}_{{\rm{DE}}})=0.\end{eqnarray} \tag{ 13 }$$

In this section we first present the general steps of the GP approach, and then we apply it in the case where the inserted data come from observations of the Hubble function.

3.1. Gaussian Process

The GP is a powerful tool allowing one to reconstruct the behavior of a function (and its derivatives) directly from given data sets (Seikel et al. 2012). The basic ingredients of the GP techniques are the covariance function (kernel), and the feature that the parameters describing it can be estimated directly from observational data (Seikel et al. 2012). Hence, one does not need to consider any specific parameterization for the involved unknown function of the model, since it can be reconstructed from observational data directly by using the cosmological equations.

In the GP one assumes that the observations of the data set are sampled from a multi-variance Gaussian distribution. Moreover, the values of the function evaluated at different points are not independent, and the connection between neighboring points is due to the covariance functions chosen in advance. The Gaussian distribution corresponds to a random variable characterized by a mean value and a covariance. Similar to Gaussian distributions, GPs should be understood as distributions over functions, determined by a mean function and a covariance matrix. Since the covariance function, for a given set of observations, can determine the relation between independent and dependent variables, the GP correlates the function at different points by using the covariance function.

There exist a number of possible choices for the covariance function, i.e., for the kernel, e.g., squared exponential, polynomial, spline, etc., which are used in various applications (Rasmussen & Williams 2005). Although there is a discussion on possible effects of the choice of kernel on the results (Rasmussen & Williams 2005; Seikel & Clarkson 2013) (in a similar way that there is a discussion of the possible effect of the choice of covariance matrix on the usual observational fittings, Eifler et al. 2009; Morrison & Schneider 2013), one commonly used choice, with good theoretical justification, is the squared exponential function (Rasmussen & Williams 2005; Seikel & Clarkson 2013):

$$\begin{eqnarray}&&k(x,{x}^{{\prime} })={\sigma }_{f}^{2}{e}^{-\tfrac{{\left(x-{x}^{{\prime} }\right)}^{2}}{2{l}^{2}}},\end{eqnarray} \tag{ 14 }$$

where σ_f and l are parameters known as hyperparameters. These parameters represent the length scales in the GP. In particular, l corresponds to the correlation length along which successive f(x) values are correlated, while we need the parameter σ_f to control the variation in f(x) relative to the mean of the process. Therefore, the covariance between output variables will be written as a function of the input ones. We mention that the covariance is maximum for variables whose inputs are suitably close. Furthermore, as can be seen from (14), the squared exponential function is infinitely differentiable, which proves to be a useful property in the case of reconstructing higher-order derivatives. Additionally, the initial and effective approach adopted in the GP to estimate the values of the hyperparameters, is based on the training by maximizing the likelihood, showing that the reconstructed function has the measured values at the data points.

Finally, we mention that the consideration of the squared exponential kernel is the most natural choice among various possibilities, given the assumption that the error distribution for the observational data is the Gaussian one (Rasmussen & Williams 2005; Seikel & Clarkson 2013).⁸

3.2. H(z) Data

In this work we will use GP techniques with Hubble data H(z) (1 + z = a₀/a, with a₀ = 1 the present scale factor). In particular, we use 30-point samples of H(z) arising from the differential evolution of cosmic chronometers, alongside 10-point samples obtained from the method of radial baryon acoustic oscillations (BAO), which allows us to extend the data range up to z = 2.4, improving also the behavior at low redshifts. In Table 1 we present the above points as they appear in Zhang & Xia (2016). We note that in principle the use of data points from different data sets should be avoided; however, as was discussed in Seikel et al. (2012), the simultaneous use of the above two data sets gives the increased statistics that are necessary for the correct application of the GP. Finally, concerning the value of the Hubble parameter at present, H₀, we use the latest released local value at 2.4% precision, namely H₀ = 73.52 ± 1.62 km s⁻¹ Mpc⁻¹ (Riess et al. 2016). We mention here that since we are using direct Hubble measurements and the local value for H₀, in our analysis and hence in our f(T) reconstruction, the H₀ tension is alleviated by construction.⁹

Table 1.  Observational Data for H(z) and their Uncertainty σ_H in Units of km s⁻¹ Mpc⁻¹

z	H(z)	σH	z	H(z)	σH
0.070	69	19.6	0.4783	80.9	9
0.090	69	12	0.480	97	62
0.120	68.6	26.2	0.593	104	13
0.170	83	8	0.680	92	8
0.179	75	4	0.781	105	12
0.199	75	5	0.875	125	17
0.200	72.9	29.6	0.880	90	40
0.270	77	14	0.900	117	23
0.280	88.8	36.6	1.037	154	20
0.352	83	14	1.300	168	17
0.3802	83	13.5	1.363	160	33.6
0.400	95	17	1.4307	177	18
0.4004	77	10.2	1.530	140	14
0.4247	87.1	11.1	1.750	202	40
0.44497	92.8	12.9	1.965	186.5	50.4
0.24	79.69	2.65	0.60	87.9	6.1
0.35	84.4	7	0.73	97.3	7.0
0.43	86.45	3.68	2.30	224	8
0.44	82.6	7.8	2.34	222	7
0.57	92.4	4.5	2.36	226	8

Note. In the upper panel we present the 30 points deduced from the differential age method of cosmic chronometers, and in the lower panel we present 10 samples obtained from the radial BAO method. The data are from Zhang & Xia (2016) (see references therein for each data point).

Download table as:  ASCIITypeset image

We use the publicly available package Gaussian Processes in Python developed by Seikel et al. (2012), and we apply it to the aforementioned H(z) data. The result of this elaboration is the successful reconstruction of H(z) and ${H}^{{\prime} }(z)$ (primes denote derivative with respect to the redshift z) in a model-independent way. In Figure 1, the numerically derived mean curves and their 1σ errors are presented.¹⁰ These reconstructions shall be used to reconstruct the f(T) forms in the next section.

Zoom In Zoom Out Reset image size

In the previous section we applied the GP to the H(z) data and we obtained the reconstruction of H(z) and its derivatives in a model-independent way, namely without assuming anything about the underlying gravitational theory or cosmological scenario. In this section we will assume that the universe is governed by f(T) gravity, hence we will use the cosmological equations of Section 2, and we will use the reconstructed H(z) and its derivatives to reconstruct the form of f(T) without any other assumption. To our knowledge it is for the first time that this is performed, since up to now in the literature a specific ansatz for f(T) was always imposed by hand and then the confrontation with the data allowed the involved parameters to be constrained.

The basic feature of f(T) gravity that makes the above reconstruction procedure easy is the fact that in the FRW geometry for the torsion scalar we have relation (8), namely T = −6H², i.e., it is a simple function of H. Thus, all the involved terms and functions of f(T) cosmology can eventually be expressed in terms of H(z) and its derivatives, which have been reconstructed through the GP in the previous section. Hence, the final step of this analysis is the reconstruction of the form of f(T) itself.

We start by expressing all cosmological equations of Section 2 in terms of the redshift z. For the time derivative of a function h we have $\dot{h}=-h^{\prime} H(1+z)$, where the prime denotes derivative with respect to z, while for f_T we have

$$\begin{eqnarray}&&{f}_{T}\equiv \displaystyle \frac{{df}(T)}{{dT}}=\displaystyle \frac{{df}/{dz}}{{dT}/{dz}}=\displaystyle \frac{f^{\prime} }{T^{\prime} }.\end{eqnarray} \tag{ 15 }$$

The next step in the application of the GP is to replace f' by

$$\begin{eqnarray}&&{f}^{{\prime} }(z)\approx \displaystyle \frac{f(z+{\rm{\Delta }}z)-f(z)}{{\rm{\Delta }}z},\end{eqnarray} \tag{ 16 }$$

for small Δz, which allows us to relate the values of f at z_i+1 and z_i. In particular, it is easy to see that from Equation (6) one acquires

$$\begin{eqnarray}&&\begin{array}{l}f({z}_{i+1})-f({z}_{i})\\ \quad =\,3({z}_{i+1}-{z}_{i})\displaystyle \frac{{T}^{{\prime} }({z}_{i})}{T({z}_{i})}\left[{H}^{2}({z}_{i})-\displaystyle \frac{8\pi G}{3}{\rho }_{m}({z}_{i})+\displaystyle \frac{f({z}_{i})}{6}\right]\,,\end{array}\end{eqnarray} \tag{ 17 }$$

where $T=-6{H}^{2}$ and T' = −12HH'. Moreover, for the matter sector we adopt the EoS parameter for regular dust, which eventually gives

$$\begin{eqnarray}&&{\rho }_{m}=\displaystyle \frac{3}{8\pi G}{H}_{0}^{2}{{\rm{\Omega }}}_{m0}{\left(1+z\right)}^{3},\end{eqnarray} \tag{ 18 }$$

due to (12). Here H₀ and Ω_m0 are the Hubble parameter and the dark matter density parameter (${{\rm{\Omega }}}_{m}=8\pi G{\rho }_{m}/3{H}^{2}$) at z = 0. Inserting these into (17), we finally obtain

$$\begin{eqnarray}&&\begin{array}{l}f({z}_{i+1})-f({z}_{i})\\ \quad =\,6({z}_{i+1}-{z}_{i})\displaystyle \frac{{H}^{{\prime} }\left({z}_{i}\right)}{H\left({z}_{i}\right)}\left[{H}^{2}\left({z}_{i}\right)-{H}_{0}^{2}{{\rm{\Omega }}}_{m0}{\left(1+{z}_{i}\right)}^{3}+\displaystyle \frac{f\left({z}_{i}\right)}{6}\right].\end{array}\end{eqnarray} \tag{ 19 }$$

Equation (19) allows us to reconstruct f(z), as long as H(z) and H'(z) are known. However, in the previous sections we were able to reconstruct H(z) and H'(z) from the Hubble data. Thus, the reconstruction of f(z) is straightforward. Finally, since both f(z) and T(z) = −6H²(z) are reconstructed, we can easily reconstruct the form of f(T) in a model-independent way.

In Figure 2 we present the reconstructed f(T) resulting from our analysis, which is the main result of the present work. Let us now try to extract information on the possible forms of f(T).

Zoom In Zoom Out Reset image size

The first and clear result is that the form $f(T)\,=-2{\rm{\Lambda }}=\mathrm{const}.$, namely the cosmological constant, lies in the central part of the reconstructed region, and in particular has the value $f(T)=-6{H}_{0}^{2}(1-{{\rm{\Omega }}}_{m0})$ (the dotted line of Figure 2), which is exactly that of the cosmological constant (f(T) = −19,267 (km s⁻¹ Mpc⁻¹)²) as expected.

Nevertheless, from the reconstructed region of Figure 2 we can obtain additional information, namely the form of the mean reconstructed curve (the black curve of Figure 2). In particular, we can see that the mean curve is not a constant, but its best fit (with accuracy R² ≈ 0.94) follows a quadratic function of the form f(T) = −2Λ + αT + ξT², with $-2{\rm{\Lambda }}=-6{H}_{0}^{2}(1-{{\rm{\Omega }}}_{m0})$ (i.e., Λ is not a free parameter) and the two parameters α ≈ −0.026 ± 0.00088, ξ ≈ (−9.68 ± 0.28) × 10⁻⁸ in units of km s⁻¹ Mpc⁻¹. Since the linear term can be removed from f(T) and be absorbed by the standard linear term that exists already in (3), the above form remains with one free parameter, namely ξ, as is the case in all viable models (that is why we did not allow for a second free parameter, namely γT³, although in this case the fitting is significantly improved at R² ≈ 0.999). Hence, in summary, the mean curve of the reconstructed procedure follows the quadratic form

$$\begin{eqnarray}&&f(T)\approx -2{\rm{\Lambda }}+\xi {T}^{2},\end{eqnarray} \tag{ 20 }$$

with ξ the sole free parameter. Note that, if we use a dimensionless free parameter, then we may rewrite (20) as $f(T)\approx -2{\rm{\Lambda }}+\beta {T}^{2}/{T}_{0}^{2}$, with ${T}_{0}=-6{H}_{0}^{2}$, and hence with the dimensionless parameter $\beta =36{H}_{0}^{4}\xi$.

However, besides the mean curve of the reconstruction, in principle any curve inside the shaded region of Figure 2 is allowed to be the true f(T) form. Hence, let us confront three viable one-parameter f(T) models of the literature (Nesseris et al. 2013; Basilakos et al. 2018) with our reconstructed region. In particular, they are the power-law (f₁CDM) model

$$\begin{eqnarray}&&f(T)=\alpha {\left(-T\right)}^{b},\end{eqnarray} \tag{ 21 }$$

with $\alpha ={\left(6{H}_{0}^{2}\right)}^{1-b}\tfrac{1-{{\rm{\Omega }}}_{m0}}{2b+1}$, the square-root exponential (f₂CDM) model

$$\begin{eqnarray}&&f(T)=\alpha {T}_{0}(1-{e}^{-p\sqrt{T/{T}_{0}}}),\end{eqnarray} \tag{ 22 }$$

with $\alpha =\tfrac{1-{{\rm{\Omega }}}_{m0}}{1-(1+p){e}^{-p}}$ and ${T}_{0}=-6{H}_{0}^{2}$, and the exponential (f₃CDM) model

$$\begin{eqnarray}&&f(T)=\alpha {T}_{0}(1-{e}^{-{pT}/{T}_{0}}),\end{eqnarray} \tag{ 23 }$$

with $\alpha =\tfrac{1-{{\rm{\Omega }}}_{m0}}{1-(1+2p){e}^{-p}}$. These models coincide with ΛCDM for b = 0 (model (21)) and for $b=1/p\to {0}^{+}$ (models (22) and (23)).

As shown in Figure 3, if we expect the above f(T) forms to lie inside the reconstructed region we obtain the following constraints on their sole parameter: −0.0005 < b < 0.0004 for the power-law (f₁CDM) model, 0 < b = 1/p < 0.15 for the square-root exponential (f₂CDM) model, and 0 < b = 1/p < 0.13 for the exponential (f₃CDM) model. Interestingly enough, these constraints are improved compared to the constraints coming from usual observational analysis using Type Ia supernovae, quasi-stellar objects, the cosmic microwave background shift parameter, direct measurements of the Hubble constant with cosmic chronometers, and measurements of redshift space distortion (fσ₈), which give −0.047 < b < 0.011 for the power-law, −0.035 < b = 1/p < 0.129 for the square-root exponential, and −0.011 < b = 1/p < 0.111 for the exponential models (Basilakos et al. 2018; Anagnostopoulos et al. 2019). Note the significant improvement by two orders of magnitude in the case of the power-law model. This is one of the main results of the present work, and shows the capabilities of the reconstruction procedure using the GP.

Zoom In Zoom Out Reset image size

Finally, concerning the quadratic form of f(T) in Equation (20), we deduce that in order to lie inside the reconstructed region its free parameter ξ should be constrained as −8.0 × 10⁻⁸ < ξ < 4.5 × 10⁻⁸ (in units of (km s⁻¹ Mpc⁻¹)⁻⁴). Equivalently, using the parameter $\beta =36{H}_{0}^{4}\xi$ discussed above, we derive that −59 < β < 33. Hence, in the present work we propose the new f(T) parameterization (20), since it can also be efficient in describing the reconstructed region that was obtained from the Hubble data through the GP analysis. The efficiency of the quadratic form was actually expected, since at late times, where H and thus T are small, every function can be expanded in T-series (Iorio & Saridakis 2012; Nashed 2015; Farrugia et al. 2016; Bahamonde et al. 2019; Chen et al. 2019).

For completeness, we close this section by presenting in Figure 4 the model-independent reconstructed forms of the dark energy density parameter ${{\rm{\Omega }}}_{{\rm{DE}}}=8\pi G{\rho }_{{\rm{DE}}}/3{H}^{2}$ from (9) and the dark energy EoS parameter w_DE from (11), as they arise using the obtained reconstructions of H(z), H'(z), and f(T).

Zoom In Zoom Out Reset image size

In this work we have applied the GP analysis and Hubble function data in f(T) cosmology to reconstruct for the first time the f(T) form in a model-independent way. In particular, up to now in the literature of f(T) gravity, as well as in the majority of gravitational modifications, physicists have assumed a specific ansatz for the involved unknown function, and used observational data in order to constrain the model parameters. However, the use of the GP analysis allows one to investigate features of the form of the involved unknown functions in a model-independent way without any assumption, using only the given observational data sets.

We applied the GP analysis for Hubble function measurements, namely for H(z) data sets coming from cosmic chronometers as well as from the radial BAO method, alongside the latest released local value H₀ = 73.52 ± 1.62 km s⁻¹ Mpc⁻¹ at 2.4% precision. Application of the procedure led to the reconstruction of H(z) and its derivative without any assumption. On the other hand, f(T) cosmology has the advantage that the torsion scalar is a simple function of H, namely T = −H², and thus eventually all cosmological equations can be expressed in terms of H(z) and H'(z). Hence, having reconstructed H(z) and H'(z) allowed us to additionally reconstruct the f(T) form itself, without any assumption. To our knowledge this is the first time where a general and model-independent reconstruction for the f(T) gravity has been obtained. Additionally, we mention that since we are using the direct Hubble measurements and the local value for H₀, the H₀ tension can be alleviated by construction in our analysis and hence in our f(T) reconstruction.

A first result of our analysis is that the cosmological constant lies in the central part of the reconstructed region, as expected. However, the mean curve of the reconstructed region is not a constant, but its best fit follows a quadratic function. Hence, inspired by this, in this work we proposed a new one-parameter f(T) parameterization, namely f(T) = −2Λ + ξT², with $-2{\rm{\Lambda }}=-6{H}_{0}^{2}(1-{{\rm{\Omega }}}_{m0})$ and ξ the sole free parameter that quantifies the deviation from ΛCDM cosmology. Moreover, fitting this form into the reconstructed region, we extracted the constraints on the free parameter as −8.0 × 10⁻⁸ < ξ < 4.5 × 10⁻⁸ (km s⁻¹ Mpc⁻¹)⁻⁴.

Additionally, we have confronted three viable one-parameter models of f(T) with the reconstructed f(T) region; these are the power-law, square-root exponential, and exponential models. As shown in the main text, we obtained improved constraints for their free parameters compared to the bounds that arise from traditional observational analyses; in particular, for the case of the power-law model the improvement was more than two orders of magnitude.

In summary, by using GPs and Hubble data we obtained a model-independent reconstruction for the f(T) form, and by fitting its mean curve we proposed a new one-parameter f(T) parameterization, namely the quadratic one. Finally, confronting three viable f(T) models of the literature with our reconstructed region, we extracted improved constraints on their parameters. These features reveal the capabilities of the reconstruction procedure using GPs. Hence, it would be interesting to apply them in other theories of modified gravity too.

We thank Canmin Deng, Rafael Nunes, Xin Ren, Yuting Wang, and Gongbo Zhao for valuable comments. We also thank the anonymous referee for constructive comments. Y.F.C. is supported in part by the National Youth Thousand Talents Program of China, by NSFC (Nos. 11722327, 11653002, 11961131007, 11421303), by CAST (2016QNRC001), and by the Fundamental Research Funds for Central Universities. M.K. is supported in part by a CAS President's International Fellowship Initiative grant (No. 2018PM0054) and by NSFC (No. 11847226). E.N.S. is supported partly by USTC fellowship for international visiting professors. This article is partially based upon work from COST Action "Cosmology and Astrophysics Network for Theoretical Advances and Training Actions," supported by COST (European Cooperation in Science and Technology). All numerics were operated on the computer clusters Linda & Judy in the particle cosmology group at USTC.