Observational Constraints on $f(T)$ Gravity from Model-Independent Data

We establish new constraints on $f(T)$ gravity models by using cosmological data. In particular, we investigate the restrictions given by the gas mass fraction measurements of galaxy clusters and transversal BAO data. Both data sets are regarded as weakly dependent on a fiducial cosmology. In addition, we also include a CMB measurement of the temperature power spectrum first peak, along with $H(z)$ values from cosmic chronometers and supernovae data from the Pantheon data set. We also perform a forecast for future constraints on the deviation of $f(T)$ models from the $\Lambda$CDM scenario by following the specifications of the J-PAS and Euclid surveys and find significant improvements on the constraints of the $b$-parameter, when compared to the results of the statistical analysis.


Introduction
The ΛCDM model has become the standard model to describe the evolution of the Universe at large scales since the discovery of its late-time cosmic acceleration [1,2]. This model is based on the existence of a fluid with negative pressure described by a cosmological constant (Λ) added to the Einstein field equations, the so-called dark energy. The model also assumes the presence of pressureless non-baryonic matter comprises most of galaxies composition denominated cold dark matter (CDM). This model is the one that best describes data from type Ia supernovae and other astrophysical objects to the cosmic microwave background (CMB) temperature power spectrum [3]. However, the increasing measurements of the latetime Universe might suggest that the framework described by general relativity (GR) is not the most general for describing gravity, motivating the idea of extensions of the standard model in a variety of ways, such as a modified theory of gravity [4,5].
One of the most discussed questions in the literature recently regarding the predictions of the standard model is the tension between data at different eras of the Universe, especially the values of H 0 (see [6] for a review) from the Planck satellite [3] and the local measurements from astrophysical objects [7][8][9]. Therefore, several modifications to GR were considered in recent years as an attempt to explain why such tensions appear and to answer other questions that the ΛCDM model does not address. One of the most studied theories involves a function of the Ricci scalar f (R) [10][11][12] in the field equations that could account for phenomena not explained by the cosmological standard model. More profound modifications from the standard picture lead to alternative options. Suppose one considers that instead of the metric, the gravitational field is described by tetrads, where the Riemann tensor, the main ingredient for the dynamics of GR, is replaced by a non-zero torsion. A change of this kind allows us to an alternative description of gravity so that the resulting field equations could lead to relevant effects in the observables. Developments towards this way led to the Teleparallel Equivalent General Relativity (TEGR), in which the torsion scalar T is related with the Ricci scalar asR = −T + B, with B being a boundary term, meaning that it is equivalent to GR at the level of field equations. A natural extension of this picture can be realized when we look at how f (R) theories are constructed, giving rise to f (T ) theories [13][14][15][16][17], with a function of T being added to the gravitational action.
A number of forms for the f (T ) function were proposed in the literature [13,18,19] where they have proved to be viable models. In Ref. [20][21][22][23][24][25][26], statistical analyses were performed to find constraints on cosmological parameters when the models are confronted with various data sets. In Ref. [27], it was shown that a power-law dependence on T can greatly alleviate the Hubble tension by increasing the associate error in the parameter, and in [24], an exponential f (T ) form was particularly preferred over the standard ΛCDM model. More recently, the impact on f (T ) constraints by different H 0 priors was investigated in [28]. In Ref. [29], a f (T ) model with exponential form was introduced as an infrared (IR) correction to GR. A complete analysis was done in Ref. [30,31], at the background and linear perturbation level, where it was found that the Hubble tension can also be alleviated. In a recent work, [26], the study of concordance with the Big Bang Nucleosynthesis helped constrain the free parameters of f (T ) models with great precision, with a deviation from the standard model at 3σ confidence level. It is also interesting to note that an extension that is being widely explored is the f (T, B) class of models [32][33][34][35][36]. Here, the boundary term B contribute to the equations of motion as being part of an arbitrary function. It makes the construction of many different models possible, and the equations of motion are generally more complicated.
In this work, we seek to constrain some of these models with data that is generally regarded as being model-independent. In particular, we want to verify the impact of including measurements of the gas mass fraction of galaxy clusters, which is the ratio between the baryonic and total mass of a given cluster. Recently, these data sets have been used to constrain cosmological scenarios [37,38], being another way of confirming cosmic acceleration. Moreover, they are a cosmology-independent way of determining cosmological parameters, especially the matter density parameter Ω m0 . In particular, we use the data set in Ref. [39], which consists of 40 points at low and intermediate redshift range of 0.078 ≤ z ≤ 1.063. Another data set that we expect to impact our analysis and produce strong constraints is the 2D baryonic acoustic oscillations (BAO) data due to its cosmological model independence. [40][41][42] Additionally, we also use H(z) values from the Cosmic Chronometers method compiled in [43,44], the position of the first peak of the CMB temperature power spectrum l 1 [45,46], and the latest Type Ia supernovae measurements of the Pantheon compilation [47].
This work is organized in the following manner. In Sec. 2, we review the f (T ) formalism, focusing on the background dynamics. In Sec. 3, we briefly review the f (T ) models that will be investigated. Sec. 4 describes the data and methodology used in the analysis. In Sec. 5, we discuss the results of the statistical analysis, while in Sec. 6 we perform a forecast on the models for future experiments. Finally, in Sec 7, we present our considerations.

Teleparallel Gravity
In this section we briefly introduces the teleparallel formalism, the generalization for a f (T ) function, and its cosmological consequences.

Formalism
Teleparallel gravity is a way of describing gravity in which the fundamental object is the tetrad e A µ instead of the usual metric tensor g µν for GR. Gravity is then described by a nonzero torsion, while the Riemann tensor, along with the non-metricity tensor, are both zero (in the teleparallel picture). A consequence of this approach is that when deriving the field equations, the Levi-Civita connection (Γ λ µν ) is substituted by the teleparallel connection Γ λ µν [48]. The metric tensor of GR is related with the tetrad as [14,16] with capital latin letters corresponding to the tangent space, while Greek letters correspond to space-time coordinates on the manifold. The teleparallel connection is written as with E λ A being the inverse tetrad. In the teleparallel picture, the connection Γ λ µν is related to the Riemannian one as being the contortion tensor that is defined in terms of the torsion tensor T λ µν , which has an equivalent role as the Riemann tensor in GR Contraction of the torsion tensor leads to which can be related to the Levi-Civita Ricci scalarR as where B = − 2 e ∂ ρ (eT µρ µ ) is a boundary term, showing an equivalence between GR and TEGR at the level of field equations. This allows us to write a similar gravitational action as the Einstein-Hilbert one. This way, we can write the field equations that can be shown to be equivalent to GR by providing the same equations of motion, according to (2.6). Furthermore, in an analogous way to f (R) gravity, one can generalize the action by introducing a function f (T ) to the gravitational Lagrangian, so that the action becomes where L m is the matter Lagrangian, and e can be identified as e = det e a µ = √ −g. Just as in f (R) gravity, f (T ) will be responsible for deviation from GR, where for instance, if the function is taken as a constant, we reproduce the ΛCDM model.
We can vary this action concerning the tetrad to obtain the field equations as where the susbscript T denotes derivatives with respect to the torsion scalar, and T ν ρ is the energy-momentum tensor, and S µν ρ ≡ 1 2 K µν ρ + δ µ ρ T σν σ − δ ν ρ T σµ σ is a superpotential that can be used to obtain the tensor scalar as T = S µν ρ T ρ µν .

Background Dynamics
To study the cosmological implications of f (T ) gravity in the context of a homogeneous, isotropic, and spatially flat universe, characterized by e A µ = diag(1, a, a, a), we see that this corresponds to the FLRW geometry characterized by the line element with H ≡ȧ a being the Hubble parameter, and ρ, P being the energy density and pressure that come from the total energy-momentum tensor, respectively. We note that if f = 0, the f (T ) formulation is equivalent to GR, while the dynamics can be modified entirely by assuming a different f (T ) function since we in the FLRW geometry have that T = −6H 2 . Before solving the equations, we should define some quantities. Since, in general we can interpret the r.h.s. of Eq. (2.11) as corresponding to the contribution of all matter components, while the Eq. (2.12) contains the contributions from the pressures+densities of the fluids, it is possible to make the following definitions for the dark energy density and pressure, respectively: Then, the dark energy equation of state can be written as . (2.14) The cosmological fluids considered will have their evolution dictated by the conservation of the energy-momentum tensoṙ ρ m + 3Hρ m (1 + w m ) = 0 andρ r + 3Hρ r (1 + w r ) = 0, (2.15) with w m and w r being the the equation of state parameters of matter and radiation, respectively; and we can find that the defined dark energy density will also follow the same conservation equation:ρ with ρ DE and P DE defined by (2.13). Since T = −6H 2 , the normalized Hubble parameter E(z) can be written as with H 0 is the present value of the Hubble parameter, and T 0 = −6H 2 0 . Also, assuming that ρ m is pressureless dust, so w m = 0, and that radiation follows w r = 1/3, we can write the Friedmann equation (2.11) as with y(z, r) being and Ω dark0 being the dark energy density parameter today, produced by the modifying f (T ) term. Note the distortion function y(z, r) that controls the effect from the modified dynamics of teleparallel gravity, where r corresponds to the free parameters of the specific model considered. The main characteristics of this function are that GR must (preferentially) be reproduced for some limit of parameter, while at the cosmological level, the concordance model ΛCDM can also be achieved (when y = 1). Numerical analysis of the main f (T ) models indicate that deviations from the standard model are generally small [24], when the model in question can reproduce the ΛCDM one, showing that different f (T ) scenarios are concordant with the standard model, and might even compete with it.

f (T ) Models
We present the f (T ) models investigated in this work. The three selected functions are well studied in the literature, and previous numerical analyses have shown that they are among the best ones preferred by data when compared to the ΛCDM model. We will see how different data can affect the predictions for each model while verifying the consistency with previous works.
• Power-law model Currently, one of the most favored by data f (T ) models is the power-law form given by [13] where α and b are the two free parameters that can be related through . The distortion factor becomes simply and then the Friedmann equation is We can easily see that b = 0 reproduces the ΛCDM cosmology. This model gives a de-Sitter limit for z = −1, and deviations from the standard model are more evident for higher |b|. However, these deviations are generally small, as verified by numerical analyses performed in past years [20,21,24]. Also, the power-law model is capable of alleviating the Hubble tension [27]. The parameter b is anti-correlated with H 0 , meaning that a larger H 0 is achieved for small b, a feature that does not happen with the other models investigated due to a strong degeneracy between parameters.
• Exponential model Another model investigated is inspired by f (R) gravity, where an exponential dependence exists, and the f (T ) function takes the form [20] where, again, α and p are dimensionless parameters that can be related though the Friedmann equation as so the distortion term is Consequently, the Friedmann equation for this model becomes where we can define b ≡ 1/p, so the ΛCDM model is recovered for b → 0 + , while the GR limit is achieved for b → +∞.
• The square-root exponential model The last f (T ) model we consider here is the exponential form studied in [18], with functional form where the α and b parameters are related as and the distortion factor becomes . (3.11) with p = 1/b. In a similar manner to the f 2 model, one can see that the limit b → 0 + reproduces the ΛCDM model, while b → +∞ corresponds to the pure GR limit.

Observational data and methods
To check the viability of these models, we will perform a statistical analysis using the Monte Carlo Markov Chain (MCMC) method, where we compare the predictions with different data sets of the cosmological observables.

Data Sets
• Gas mass fraction data The first data used in this work is the cluster mass gas fraction f gas ≡ Mgas M total [49-54] (one can check [39] for further references). Since these clusters can be assumed as containing a good part of the total content of non-relativistic matter in the Universe, we can initially, approximate f gas as where Ω b is the total fraction of baryonic matter, while b gas (z) is some function that expresses how different the cluster mass gas fraction is from the cosmic one. As f gas ∝ d L d

1/2
A and following the Ref. [39,55], we can use the cosmic distance duality relation d A = d L /(1 + z) 2 to write f gas in terms of the angular diameter distance where D A (z) is the angular diameter distance for a given model, normalized by a fiducial model that is taken as a ΛCDM one with Ω m0 = 0.3 and H 0 = 70 km/s/Mpc for the data we are using. The A(z) factor is the angular correction between two models, which is usually close to one, but can be modeled as where η is estimated as η = 0.442 ± 0.035 [39]. K(z) and γ(z) are respectively the calibration bias and the depletion factor, where the former takes into account instrumental inaccuracies as well as astrophysical effects in the cluster mass, while the latter measures the depletion of hot gas in the cluster relative to the baryon cosmic fraction. Some works in the literature have investigated the possible variation of these quantities with redshift (in partcular γ) [56,57], but in this analysis, we take them as constants, as estimated by hydrodynamical simulations; therefore, we use the values γ = 0.848±0.085 [58], and K = 0.96±0.09±0.09 [59]. By using these three parameters (η, K, γ as discussed, the χ 2 function for gas mass fraction measurements is given by In this expression, f gas (z i ) represent the theoretical predictions given by Eqs. (4.2-4.3), f obs gas are the observational values; the uncertainties σ f gas,i have the effective form where σ 2 obs,i are the uncertainties associated with the data. We use the following data set for the analysis. In [39], the fraction was derived for 40 cluster measurements at the radius r 2500 1 improving the previous work done in Ref. [55]. These points cover the redshift interval of 0.078 ≤ z ≤ 1.063. A recent application of these data points in constraining cosmological parameters is described in [38] for the ΛCDM and wCDM models, in a way that we can compare our results with theirs, especially in the determination of the Hubble parameter H 0 . It is good to mention that other measurements of f gas are available in the literature [60][61][62][63], from lower to higher redshifts (0.0473 ≤ z ≤ 1.235), but measured in the radius r 500 , which will not be used in the present work.
• BAO 2D data As in the case of gas mass fraction data, the BAO data we will use in this analysis is computed in a way that can be regarded as almost model-independent. The method, presented in Refs. [40][41][42], involves the 2-point correlation function for a distribution of galaxies, where only the angular separation is considered in redshift shells of the order δz 10 −2 . This allows one to obtain information on the BAO transversal signal without the effect of a fiducial cosmology 2 , and can be used to test different cosmological scenarios. The expression gives the BAO angular scale θ BAO where r s is the comoving sound horizon, obtained as 3 where Ω b0 is the present baryon density parameter, and Ω γ0 is the present photon density parameter. The redshift at the drag epoch z d is estimated by the fitting formula [64]  The total BAO χ 2 function (χ 2 BAO ) is then (4.10) These data points have been used previously in different investigations. For instance, in [65], the ΛCDM and CPL models were analyzed along with Planck data. It was shown that a dynamical dark energy scenario in this context provides a value for H 0 that is compatible with local measurements. Following this work, the same data was used with H0LiCOW data [66] to obtain constraints on the H 0 − r d plane independently of CMB data and investigate the impact on spatial curvature. In [67] cosmological constraints were obtained by imposing observational and thermodynamics limits on interacting dynamical dark energy models.
• Type Ia Supernovae (SNe) data We also use in this analysis the SNIa Pantheon compilation [47]. In particular, we consider the binned version where the 1048 points are compacted to 40, which span the redshift interval 0.01 ≤ z ≤ 1.6. The χ 2 function is given as where C −1 SN e corresponds to the inverse covariance matrix of the data, and ∆µ = µ i − µ i,th is a vector with the difference between the observational and theoretical distance modulus. The distance modulus is defined as µ = m B − M, where µ B is the observed apparent magnitude at a given redshift, while M is the absolute magnitude which is treated as a nuisance parameter in the statistical analysis. This is compared with the theoretical form calculated via where D L (z) = (1 + z) z 0 dz H(z ) is the luminosity distance.
• H(z) data We use measurements of the Hubble parameter obtained from the differential age method, also known as cosmic chronometer (CC) data. This method of measuring the differential age of galaxies allows us to determine the Hubble parameter at a certain redshift without assuming a specific model. Here, we will consider 31 points cataloged in [43], and compiled in Table 1 of [44] spanning the redhsift range of 0.07 ≤ z ≤ 2.
The χ 2 function is constructed as (4.13) • CMB data The last data set used in this work is the Planck CMB data encoded on the first peak of the temperature power spectrum, indicated by l 1 , expressed as [45] (4.14) with r dec being the acoustic sound scale, and all quantities are evaluated at the decoupling redshift z dec [68]. The measured value of the peak we use is l 1 = 220.0 ± 0.5 [46] Therefore, to analyse the impact of the f gas and BAO data, we consider four total χ 2 functions: χ 2 = χ 2 Base +χ 2 fgas , χ 2 = χ 2 Base +χ 2 BAO , χ 2 = χ 2 Base +χ 2 fgas +χ 2 BAO along with all data sets combined, and where χ 2 Base ≡ χ 2 SN e + χ 2 CC . We assume uniform priors on H 0 , w, b and Ω m0 and a Gaussian prior on the baryon parameter density of Ω b h 2 ≡ ω b = 0.0226 ± 0.00034 [69]. To perform the MCMC analysis, we use the emcee sampler [70], and the GetDist [71] Python module to plot the results.

Model Selection
After determining the parameters posteriors distributions for each model, we must use a way to compare them, which will help us to determine which model is more favored by the data used. The most robust estimator used in cosmology for statistical comparison is the Bayes factor, the ratio between the Bayesian evidence of a model of interest and a reference model. We also compute the value of the Akaike Information Criteria (AIC) [72], which, under the assumption of at least near gaussianity of the posterior distribution, it is given as [73] AIC ≡ −2 ln L max + 2k(k + 1) In (4.15), L max is the value of the maximum likelihood for a given model.  Table 1. Cosmological constraints for all models investigated with 1σ uncertainties. We divide the results into Base+f gas , Base+BAO 2D , Base+f gas +BAO 2D and Base+f gas +BAO 2D +CMB data sets. same classification as [24], where ∆IC ≤ 2 corresponds to statistical compatibility between models, 2 < ∆IC < 6 represents a tension between them, while ∆IC ≥ 10 represents robust evidence against the model we want to compare with the reference one. We also use the Bayes' factor as an evidence-based statistical estimator for model selection. This quantity considers not only the best-fit point (the minimum χ 2 parameters values) but also the entire probability distribution. The definition of the Bayes' factor, B 01 , is the evidence ratio between two models: being E i the evidence in the Bayes' theorem for the i−model and E 0 the evidence for a reference one. As in the AIC criterion, the comparison is performed with a reference model and a qualitative inference is interpreted by the Jefreys' logarithmic scale [74]. In this scale, the logarithm of the Bayes' factor determines the preference for a model with the higher Bayesian evidence. The characteristic values of the scale are: | ln B 0i | < 1, 1 < | ln B 0i | < 2.5, 2.5 < | ln B 0i | < 5, and | ln B 0i | > 5 for inconclusive, weak, moderate and strong evidence, respectively.

Results
The results of the statistical analysis for all models considered are displayed in Figures 1, 2,  3 and 4 and Tables 1 and 2, where the values for the statistical criteria in Table 2 correspond to the combination of all data sets. For the standard ΛCDM model, we realize the following: When using the combination Base+f gas we have the lowest value for H 0 , of 68.962 +1.688 −1.741 , compatible with the value of Planck [3] of 67.36 ± 0.54 Km/s/Mpc at 1σ confidence. This value is also considerably higher than the one obtained in [38]. The present matter density, on the other hand, is the highest, of Ω m0 = 0.301 +0.012 −0.011 , as suggested by the anti-correlation of the parameters in Figure 1 (grey contour). When we consider Base+BAO 2D , we note a considerable increase in H 0 and a decrease in Ω m0 , where H 0 goes more towards the R19 value [8], and there is an improvement in the error bars, when considering BAO 2D data. For the third combination, Base+f gas +BAO 2D , we note a slight increase in H 0 , but with the uncertainties essentially preserved. However, there is now a significant improvement in the uncertainties of Ω m0 . Finally, for all data sets combined, there is a great decrease in the error bars for H 0 , with similar results as in [24], but in our case, the values of parameters remain almost the same.
We note similar behaviors in the constraints of the f (T ) models, where the lowest value of H 0 is always obtained when we consider the combination Base+f gas , while the highest is achieved for Base+f gas +BAO 2D . For the f 1 power-law model, we have H 0 = 69.602 +1.618 (for the Base+BAO 2D combination) and H 0 = 70.990 +1.114 −1.217 (for the Base+f gas +BAO 2D combination), being quite different from the results in Ref. [28], where the 'Base' data set was used, but with another BAO data, and closer to Ref. [24], where Base+f σ 8 was used. A good improvement in the H 0 constraint was obtained by adding CMB data. As for the parameter b, the combination Base+f gas +BAO 2D leads to 1σ concordance with the ΛCDM scenario. Such a conclusion was also obtained from previous analyses in the literature. One interesting feature of the power-law model is the ability to greatly alleviate the H 0 tension due to its anti-correlation with b. In our results (Figure 2), this anti-correlation is preserved, while there is an inversion of correlations in the H 0 −Ω m0 plane, when we consider Base+f gas and Base+BAO 2D data sets. For all data sets combined, we note that this specific correlation is not as evident, so the value of Ω m0 is better determined. In addition, the value of b agrees even more with the standard model, where this time, we obtain a small positive value.
As shown in Figures 3 and 4 , the results are also similar for both exponential f 2 and f 3 models. For all data sets combined, H 0 is better constrained when compared with the power-law model, while the parameter b does not have the ΛCDM limit at 1σ level. These general results are also present in the literature and are part of the explanation as to why the power-law model can solve the H 0 tension, while the f 2 and f 3 models cannot [27]. We then see that using these data for constraining f (T ) models leads to results consistent with recent previous studies. We also achieve a similar level of restriction as other data sets available, despite the larger associated error bars, as it is in the case of f gas and BAO 2D data.
To conclude, we discuss the statistical results from the Bayes factor and AIC. The AIC criterion indicates the ΛCDM model as the best one, followed by f 1 , f 3 and f 2 models, which is consistent with the previous results in the literature, where the power-law model is the best one. From the scale described in the previous section, all models are in mild tension with ΛCDM for these data. Also, as pointed out in Ref. [24], a slight difference of ∆IC between models makes it difficult to establish which of the competing models is the best. Hence, we also have a statistical equivalence between the f (T ) models. As for the Bayes'  Figure 1. 1σ and 2σ confidence contours and posterior distributions for the ΛCDM model. The grey contours represent the analysis with SNe+CC+f gas data, the purple contours correspond to SNe+CC+BAO 2D data, the green contours correspond to SNe+CC+f gas +BAO 2D , and the blue contours are referent to all data combined, SNe+CC+BAO 2D +f gas +CMB. factor, we have different results. Although the ΛCDM is still preferred, we note a significant preference for the f 3 model, when we look only at the f (T ) models; the f 1 models seem to be performing the worst in the light of this criterion.

Forecast
The next generation of surveys mapping the Large Scale Structure of the Universe will obtain tighter constraints on cosmological models. These surveys will be important to distinguish between modified gravity theories and dark energy models by considering precise data in a wide range of cosmic history. In this section we follow the specifications of the J-PAS and Euclid surveys to quantify the future constraints on f (T )-gravity that will be obtained Base + fgas Base + BAO Base + fgas + BAO Base + fgas + BAO + CMB Figure 2. Same as Figure 1, but for the f 1 power-law teleparallel model. from the radial BAO signal [75-80] 4 . For this purpose, we use the expected H(z) relative errors to simulate Hubble parameter data considering the fiducial cosmology given by the results presented in Table 1 (Base + f gas + BAO 2D + CMB). We replace the BAO 2D real data with the simulated H(z) data in the statistical analysis to avoid double counting of the same observable and maintain the rest of the data sets. It is worth mentioning that this approach considers the constraints obtained from the distribution of galaxies and their effect in conjunction with the other observables used in this work in their current state.
In Table 3, we present the results of the f (T )-gravity parameter constraints using the J-PAS-like and Euclid-like H(z) estimates. As is shown, the most essential improvements on the b-parameter constraints occur for the f 1 and f 3 models in comparison with the results presented in Table 1. Such results are particularly relevant for the f 3 model because it would be possible to measure deviance of the ΛCDM model in ∼ 2σ. Another important point is that the constraints for these two surveys are similar, being the Euclid H(z) estimates are more precise while the J-PAS estimates cover a wider redshift range.

Conclusions
This work has explored the statistical viability of some f (T ) gravity models with the essential cosmological datasets. The analysis considered the H(z) from Cosmic Chronometers, Type Ia supernovae from the Pantheon set, and the first peak of the CMB temperature power spectrum, with the addition of other data regarded as being model-independent, in our case, the gas mass fraction of galaxy clusters and radial BAO data. Our main goal was to study the impact of these data in constraining the cosmological parameters, especially  Table 3. Results of the future constraints for each f (T )-gravity parameter model using the J-PAS-like and Euclid-like H(z) estimates and real data. Base + fgas Base + BAO Base + fgas + BAO Base + fgas + BAO + CMB Figure 3. Same as Figure 1, but for the f 2 exponential teleparallel model. the b parameter present in the f (T ) functions that control the deviation from the standard ΛCDM scenario. We have found that the free parameters can be constrained with similar accuracy as previous works. Although the associated error bars in the f gas and BAO 2D data are considerably large, it would be interesting then to observe the impact of using such data in the constraints of other modified gravity models. From a statistical point of view, we have seen that the AIC/Bayes' factor criteria prefer the standard scenario, so the ΛCDM remains the best model. However, among the f (T ) models, the AIC indicates a statistical equivalence, especially between the f 1 and f 3 models, while the Bayes' factor shows the f 3 model as the best one. Finally, we have performed a forecast for the statistical analysis using two next-generation galaxy surveys: J-PAS and Euclid, to predict the improvement in the measurements of the b-parameter of the three f (T ) models. Compared with the results of Table 1, we have found a significant improvement in the error bars of b, especially for the f 1 model. We also note Base + fgas Base + BAO Base + fgas + BAO Base + fgas + BAO + CMB Figure 4. Same as Figure 1, but for the f 3 square-root exponential teleparallel model. that for both J-PAS and Euclid, similar constraints with the simulated H(z) measurements are found, and in the context of the f 3 model, a ∼ 2σ deviance from the standard ΛCDM model is observed.