Cosmology of cubic galileon in modified teleparallel gravity

In this present paper, we study the cosmological evolution of the cubic galileon along with modified teleparallel gravity at perturbed and non-perturbed levels. We show the dynamical equations of motion and investigate the evolution of different cosmological parameters by using the dynamical variables analysis. In addition, a detailed analysis of different cosmological evolution in the matter, radiation and de Sitter eras is presented by solving the dynamical equations numerically. In our analysis, we find that the equations of motion in the Friedmann–Robertson–Walker (FRW) background metric is characterized by a stable de Sitter era and a tracker solution in which Hφ˙\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H{\dot{\varphi }}$$\end{document} is always constant. We find also that the equation of state of dark energy associated to the proposed model in this work can deviate from − 2 at the matter era. Moreover, the conditions of avoiding ghost and Laplacian instabilities in our model are derived; then we show that the model is free of these instabilities. Furthermore we place an observational constraint on the parameters of the model through Monte Carlo numerical method using growth rate and observational Hubble data. Finally, using the best-fit values of parameters in the model we compare our growth rate of matter perturbation with the prediction of Λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varLambda $$\end{document}CDM model and the latest measurement.


Introduction
In the last two decades, comprehending the evolution of dark energy and dark matter in the context of modified gravity has received much attention [1,2]. In fact, the modifications added to General Relativity (GR) are interpreted as a dark energy that can be responsible for the accelerate expansion of the universe. The simplest way to describe dark energy is by supposing a constant Λ energy density with an equation of state (EoS) w DE = − 1, but it causes several problems, such as the fine-tuning and coincidence problems [3]. Thus, a a e-mail: boumaza14@yahoo.com (corresponding author) large family of modified gravity theories have been proposed to find the origin of dark energy [4] in which dark energy's EoS deviates from the value − 1 at late time. One can classify these theories in two classes: the first class introduces a matter field source as dark energy, like quintessence [5,6], Brans and Dicke [7], k-essence [8] and Degenrat-Higher-Order-Scalar-Tensor (DHOST) theories [9]. And the second one replace the Ricci scalar R by a general function like f (R)-gravity [10,11] and f (R, T ) [12] where T is the torsion scalar.
Horndeski gravity, proposed in 1974 [13], is the most general scalar-tensor theory with a single scalar that leads to second-order equations which allows us to avoid the Ostrogradsky instabilities [14]. Depending on the model belonging to Horndeski modified gravity Lagrangian, like the covariant galileon [15] and Dvali-Gabadzde [16] models, one can have a deviation of the effective gravitational coupling from the Newtonian gravitational constant. Therefore, the growth matter rate will also deviate from that in the ΛCDM model which might be a solution to the σ 8 -tension. However, in order to satisfy the propagation speed of gravitational waves (GWs), detected by GW170817 [17] from a binary neutron star merged together with gamma-ray burst GRB 170817A [18], Horndeski theories need to be reduced to the form L H = G 4 (ϕ)R +G 2 (ϕ, X )+G 3 (ϕ, X ) ϕ [19]. The cubic galileon [20][21][22], which is included in the Lagrangian L H , contains self-acceleration de Sitter attractors responsible for the accelerated universe, but it is in tension with the joint observation data of SN Ia [23,24], CMB [25,26], BAO [27], large-scale structure, and weak lensing. This is due to the bad behavior of the EoS corresponding to dark energy where one found that w DE = − 2 at the matter era for the tracker solution [28].
An alternative formulation of gravity, called teleparallel gravity (TG) [29], is proposed by Pellegrini and Plebanski [30] based on Moller's work [31,32]. In this theory, we use the tetrads fields e a μ rather than the metric field g μν where the gravity is considered as a gauge theory unlike, in GR,where the gravitational interaction is geometrized. In fact, TG can be unified with the fundamental forces of the Standard Model [33] because it is a gauge theory of translation.
The equations of teleparallel gravity are equivalent to those in General Relativity (GR) because of the divergence term that appears in the relation between scalar torsion T and Ricci scalar R: R = T + ∂ μ (eT μ ), where e is the determinant of the tetrads and T μ is the torsion vector. However, the field equations of these two theories will not remain the same if we consider nonlinear terms of T , which have been studied in the literature [34][35][36]. The theories described by the Lagrangian f (T )-functions [37,38] as well as a scalar field term can be considered to lead to f (T, ϕ) [39][40][41][42][43][44]. It is well known that f (T, ϕ)-gravity is a second-order theory whereas the background equations of f (R, ϕ)-gravity contain fourthand third-order derivatives of the metric. Moreover, in modified teleparallel gravity the Lorentz invariance is violated and thus introducing new extra degrees of freedom is necessary [45,46]. Indeed, unlike in GR, these new variables are involved at the perturbation level which means that the effect on the large structure in the universe is significant.
In Ref. [47], the Horndeski theory in a teleparallel gravity frame work was explored keeping the second-order equations of motion. One shows that the theory has a much richer phenomenology and there are new terms in the Lagrangian which can explain dark energy without the cosmological constant. Afterward, the theory was generalized to get c T = 1 without dropping the terms G 5 (ϕ, X ) and G 4 (ϕ, X ) in Ref. [47]. Inspired by these two papers, we study the cosmological evolution of the cubic galileon in modified teleparallel gravity at the late time of the universe. We propose to study different cosmological parameters, the conditions of ghosts and Laplacian instability, and we constrain the model with the growth rate factor (GRF) and observational Hubble data (OHD).
The present paper is organized as follows. In Sect. 2, we present a general model of cubic galileon in Teleparallel gravity frame work (CGTG) and show the equations of the model in a homogeneous and isotropic universe. In Sect. 4, we calculate the perturbation equations corresponding to the perturbed background, scalar field and matter density. In Sect. 5, we address phase-space dynamical equations and we study the analytic and numerical evolution of dynamical systems. In Sect. 6, a numerical analysis got from Monte Carlo Markov chain is presented to derive the best-fit value of each cosmological parameter. We compare the growth rate evolution prediction for our model with that in the ΛCDM model. Finally, we conclude with the obtained results in Sect. 7.

The model
Let us assume a theory of modified gravity in which we consider the cubic galileon field and modified teleparallel gravity in the following action: with where G is the gravitational coupling and X = − 1 2 ∇ μ ϕ∇ μ ϕ is the kinetic therm of the scalar field ϕ, e = √ −g is the determinant of the tetrad field e a μ and T is the torsion scalar (for more details see Appendix A). L m and L r are the Lagrangian of matter and radiation, respectively. Our proposed model is invariant under the shift transformation ϕ −→ ϕ + const. In addition, the gravitation wave propagation speed, in our model, is constrained to the speed of light, which is proved in Ref. [48]. In other words, the model satisfies the constraint observed by the detection of gravitational waves (GW) GW170817 [17] from a binary neutron start merger together with the gamma-ray burst GRB170817A [18]. The general equations of motion corresponding to our model are obtained by varying the action (1) with respect to tetrads and scalar field. Thus, it follows that where The first term in the action (3) corresponds to the motion's equations of the tetrad field and the second corresponds to the equation of the scalar field. These equations can be written as where where i = {m, r } and ρ i , P i and u μ are energy density, pressure and the four velocity vector, of the radiation and matter present in the universe, respectively. Note that the tensor of energy-momentum of the matter fluid is conserved. The equation of conservation is given by Since the energy-momentum tensor is symmetric and invariant under local Lorentz transformation, the antisymmetric part in the tensor (7) must vanish. Thus, it leads us to the constraint The general equation of scalar field of our model is calculated by This equation is the identical to the scalar field's equation of the cubic galileon model addressed in Ref. [20]. In our model, we have the same equation of the scalar field but it not for the background.

Cosmological equations
Now, we consider a flat, homogeneous and isotropic spacetime where the metric is written as where a(t) is the scale factor. This metric can be expressed by the following tetrad: Note that, since we are working in the framework of the teleparallel equivalent to general relativity (TEGR), we can define the torsion tensor in (126) by a vanishing spin connection (ω a bμ = 0). In addition, as is shown in Ref. [38], even with an arbitrary tetrad in an arbitrary coordinate system along with the corresponding spin connection, the field equations in covariant f (T ) gravity are always the same. In our paper, we will use the diagonal tetrad form to study the cosmological evolution of our model. In this case, the scalar torsion is obtained by T = 6(ȧ/a) 2 and the background equation of motion in FRLW spacetime is obtained by replacing the metric (12) in (5) and assuming that the scalar field depends only on the cosmic time t. In this case, Eqs. (5) and (11) are reduced to where H =ȧ a is the Hubble parameter and the dot represents a derivative with respect to t. Equations (14) and (15) can be written in terms of dark energy densities and pressure as follows: where the density ρ DE and the pressure P DE are given, respectively, by Moreover, the equation of state of dark energy can be defined as Finally, substituting the metric (12) in (11), we geṫ where P m = 0 ( for matter) and P r = 1 3 ρ r (for radiation).

Perturbed equations
In this section, we present the perturbation equations of motion for the model we proposed in Sect. 2. However, we start by deriving the equations of scalar perturbation by considering the scalar part in the scalar-vector-tensor (SVT) decomposition of δe a μ which is given by where Φ, Ψ , F, B and α are scalar quantities and depend on the spacetime. We observe the appearance of new extra degrees of freedom in the scalar teleparallel quantities, which is not familiar for the GR metric. In fact, the components B and α play a crucial role in the cosmological perturbation theory, since the Lorentz symmetry is violated in F(T )-gravity.
In addition, we chose this particular form in order to introduce the Newton gauge and flat gauge in the following. The metric corresponding to the tetrad (23) is given by In order to simplify the equations of the perturbation, we fix the gauge by setting B = 0 and F = 0 (Newtonian gauge). In addition, the scalar field, the matter density and four velocity of matter are decomposed into homogeneous and perturbed parts: x, y, z) and u μ −→ δu μ (t, x, y, z), respectively. Substituting the expressions of tetrads (23) in the Eq. (5), the components of the perturbed equations are obtained as follows.
The energy density (0, 0): The energy flux and momentum density (0, i): The momentum flux (i = j): The anisotropic parts (i = j): This equation shows that the extra degree α implies Φ = Ψ in the case of f T T = 0. Using Eq. (11), we also derive the first order of the scalar field equation: Δδϕ Note that this equation is the same equation as derived for the cubic galileon [20]. In teleparallel gravity, the constraint (10) is satisfied at background level but in the linear perturbation becomes Finally, the equations of conservation corresponding to matter at linear perturbation is written as

Sub-horizon approximation
In order to show the evolution of the growth rate and confront our model with the observations of large-scale structures, one must derive the expression of the effective gravitational coupling. To do so, we study the cosmological perturbation of our model on sub-horizon scale k a H and we define the gauge-invariant contrast as where in deep Hubble radius δ m ≈ δρ m ρ m . Thus, taking the derivative of Eq. (31) and using Eqs. (32) and (33), one gets in Fourier spacë where k is a comoving wavenumber. On a sub-Hubble approximation scale, Eqs. (25) and (29) are reduced to Solving Eqs. (28), (30), (35) and (36) with respect to Φ, Ψ , δϕ and α, we take the limit k −→ ∞, we get Hk 2 (8c 2 Hφ f T + c 2 2φ and the effective gravitational potential, By taking into account only the first-order approximation of the Eq. (37), we can define the effective gravitational coupling as In the case c 2 = 0 and f T = 1, we recover General Relativity (G eff = 1) which is also the case of quintessence and kessence models [8]. However, if we consider c 2 = 0, the effective gravitational coupling takes the form This result agrees with G eff found in Ref. [44]. Finally, if we set f T = 1, we obtain The latter result corresponds to the gravitational coupling of the cubic galileon studied extensively in Ref. [22].

Ghost and Laplacian instability
Now, we derive the condition for the absence of ghost and Laplacian instability in the small-scale limit by taking the flat gauge in consideration where we set B = 0 and Ψ = 0. So, after expanding the action (1) up to second order and integrating by parts, the action is written as where the explicit form of A 1,...,5 and Σ is given by For the matter part, we used the Schutz-Sorkin action [49], which has been studied in detail by expanding up to second order as [50] where c 2 i is the sound speed of the matter fluid. Varying the action (44) with respect to Φ, F and α, we obtain the constraints Note we assumed that c 2 m = 0. And the equation of perturbed matter and radiation in flat gauge are obtained by varying the action (44) with respect to v and δρ m . Thus, the equations are written as follows: By solving Eqs. (48)(49)(50)(51) for the quantities α, Φ, F and v, then plugging the results in the action (44), we can write the order action in Fourier space as where K , N , M and B are 3 × 3 matrices and After integrating the terms containing the k 2 /a 2 in B and having absorbed them into N , we get the non-vanishing matrix components, in the small-scale limit (k −→ ∞), Providing that there no off-diagonal components in the matrix K and N the matter and the scalar field are decoupled.
To avoid the scalar ghosts, the matrix K must be definite positive. Proving that K 22 and K 33 are positive, one requires that K 11 > 0. The dispersion relation corresponding to the action (53) takes the form The solutions of this equation correspond to the matter speed squared c 2 m and c 2 r and the speed squared of the perturbed scalar field δϕ defined by c 2 S = N 11 /K 11 . Since c 2 m > 0 and c 2 r > 0, the Laplacian instability is absent only if c 2 S is positive.

Tensor perturbation
The derivation of the speed of gravitational waves requires one to use the perturbed tensor h i j as component of the tetrads, which obey the traceless and divergence conditions (h ii = 0 and ∂ i h i j = 0) and it given by Hence, the perturbative components of the metric become Expanding the action (5) up to second order and after integrating by part, we obtain The variation of this equation with respect to h i j = 0 yields the following equation: where the dimensionless parameter is defined as We can see that the speed of GWs is equal to the speed of light. Thus, we showed that the model is not disfavored by the experimental observation constraint of GW170817. However, unlike the GR theory the parameter β T deviates from zero and grows in time when the scalar and Hubble factor evolves. The parameter β T will affect the amplitude of the GR waveform given by [51] where h G R is the amplitude of gravitational waves emitted by a binary system of masses orbiting one another.

Dynamical analysis and numerical solution
In this section, we will give the form of the function f (T ); then we use the expressions derived in the last sections. In addition, we will solve the equations of motion numerically in terms of the dynamical variable. Thus, we will see the cosmological effect of the terms in our model.

Dynamical equations
Now, we present the dynamical analysis of our model by choosing a particular function f = T + κ T m , where κ and m are real constants. Thus, the equations of motion (14) and (17) are reduced to If we consider that the first order of the Hubble parameter and scalar field are constant at the de Sitter epoch, we can derive the constant c 2 and κ. By imposingφ = 0 andḢ = 0 at the de Sitter era and by solving Eqs. (62) and (64) for c 2 and κ we find where the de Sitter point is represented by the "ds" subscript. It is convenient to study the FRLW equation in the framework of dynamical variables defined as follows: By using the variables (65), we can rewrite Eq. (62) as where Equations (64) and (16) can also be written as By solving these two equations for H and ϕ , we obtain In terms of dynamical variables and by using Eqs. (70) and (71), the expressions of w DE , G eff , c 2 S and β T read G eff = 6(2m − 1) 3x 2 1 (mx 3 + 6) − 2x 1 (2mx 3 + 9x 2 + 12) where By solving Eqs. (68) and (69) for H and ϕ and substituting the results in the derivatives of x 1 , x 2 , x 3 and Ω r with respect to E-floding time N = ln(a(t)), we find where the prime represents the derivative with respect to N . Note that we used the constraint Ω DE + Ω m + Ω r = 1 to eliminate Ω m . In order to avoid divergences during the cosmological evolution, we require that The critical points of the system are derived by solving the equations x 1 = 0, x 2 = 0, x 3 = 0 and Ω r = 0 [52,53] and we summarized their values and the stability properties in Table 1. In order to study the stability of this point we perturb the dynamical equations of x 1 , x 2 , x 3 and Ω r around each fixed point obtained: 3 , Ω r = Ω c r + δΩ r . Thus, we get the following linear differential system of equations: where H is Hessian matrix and − → x is given by We see there is a stable and an unstable fixed point, corresponding to the de Sitter solution. We see also there are two critical saddle points (O 2 and O 6 ), with w eff = 0, which correspond to the matter era solution. Moreover, the points O 7 and O 4 , with w eff = 1/3 (radiation domination), are saddle but this is not the case for O 4 when m 1. In order to understand the role of the fixed points we will estimate the evolution of the cosmological parameters analytically in the three eras. Then we solve the differential equations of (79)-(82) numerically, in the next section. However, since the point O 5 has the strange value w eff = 1/5 and the point O 3 is saddle, they will not be taken into consideration.

Tracker solution
It is obvious that Eq. (79) is characterized by the solution x 1 = 1, which is called a tracker solution [20]. In fact, along the cosmological evolution of this quantityφ H remains constant. In this case the dynamical equations are reduced to From Eqs. (87) and (88), we can deduce the relation And from (87) and (88), we can also deduce Integrating these last two equations with respect to N leads to Since (N = 0). So, if we want to keep the dark energy density sub-dominated at the radiation epoch, we must impose the condition m < 1. Furthermore, it is easy to have an equation corresponding to Hubble parameters by substituting the expression of Ω r , Ω m , x 2 and x 3 in the constraint (66). Thus, we obtain Along the tracker the equation of state of dark energy and the effective total equation of state, defined by w eff = −1− 2 3 H , evolve as In the matter and radiation epoch the effective total equation of state (| x 2 | 0 and | x 3 | 0) is reduced to w eff = Ω r /3. As regards the critical point, this can be seen as a transition from the point O 6 to the point O 7 . However, the evolution of EoS of dark energy, at the early epoch, can take two regimes: The first regime is characterized by | x 3 | | x 2 | 0 in which we have This behavior is identical to that of the the cubic galileon and is not favored by the experimental observations [21]. But in the second regime, where we the condition | x 2 | | x 3 | 0 is satisfied in the matter and radiation epochs, we have Here, since the EoS takes the value m − 1 and −1 + (4m)/3 at the matter and radiation areas, respectively, one can have w DE > −2 in the matter area for m < −1. The quantities associated to c 2 S , Q S and β T for the tracker solution are given by Table 1 The critical points with the eigenvalues found by using the Hessian matrix. In the last column we give the value of w eff for each critical point (0, 0, 0, 1) ( 5/2, 1, 1/2, 4(m − 1)) Saddle for m < 1 and unstable for m 1 1 / 3 with During the cosmological evolution of the universe in the matter and radiation eras (| x 2 | 0 and | x 3 | 0), the quantities c 2 S , Q S and β T evolve as c 2 S = (Ω r +4)/6, Q S = 3x 2 and β T 0, respectively. Hence, for x 2 > 0, the ghost and the Laplacian instability are automatically avoided in the matter and radiation epochs. Finally, expanding the effective gravitational coupling about x 2 = 0 and x 3 = 0 gives For x 2 x 3 , we observe that G eff depends on m where G eff can be smaller than 1 but this is not the case when x 3 x 2 .

Small regime
We can study this regime by considering a linear perturbation of x 1 , x 2 and x 3 around zero. Thus, it follows that We see that the points O 4 and O 2 correspond to critical points of Eqs. (102)-(105) and thus a transition O 4 to O 2 may occur. By integrating these equations, we obtain the analytical solution as follows: with where C 1 , C 2 , C 3 and C 4 are constants of integration. From Eq. (107), we can characterize the radiation era by a 1/C 1 and the matter era by a 1/C 1 . Thus, the dynamical variables x 1 and x 2 increase in the two eras but the evolution of x 3 depends on the parameter m. The dynamical variable x 3 increases if m < 1, remains constant for m = 1, and decreases in the case m > 1. In this regime, the EoS of dark energy given by Eq. (72) becomes For the case | x 3 | | x 2 |, we have w DE = (3 − Ω r ) /12 which is the same in the cubic galileon's EoS in the small regime [20]. In addition, if | x 2 | | x 3 |, the equation of state is reduced to w DE = (mΩ r + 3m − 3) /3. In fact, the EoS evolves as 1/6 (radiation area) −→ 1/4 (matter area) in the first case, whether it is evolves as −1 + (4m)/3 (radiation area) −→ (m − 1) (matter area) in the second one like we found in the tracker. By taking the limit x 1 −→ 0, x 2 −→ 0 and x 3 −→ 0, the propagation speed of scalar field perturbations and the parameter β T are given by When x 2 > 0, we observe that the ghost and Laplacian instability are absent where c 2 S = 1/4 and c 2 S = 5/24 at the matter and the radiations epochs, respectively. In addition, in these epochs the sign of the parameter β T depends on the value of m and the sign of x 3 . In the small regime, G eff is expanding to We have the same result as found in (101) for | x 2 | | x 3 | but we have a small change in the case | x 3 | | x 2 |.

Numerical analysis
In order to show the cosmological evolution of the dynamical variables in more detail, we solve the system of equations Hence, the cosmological evolution of the dynamical variables depends on the numerical values of x 1 and m. If we assume that x (0) 1 = 1, we get the tracker solution discussed above because the variable x 1 is equal to one along the three cosmological areas in this condition. Furthermore, as depicted in the right of Fig. 1, the value of x 2 at the de Sitter epoch depends on the initial conditions. We note that the dashed line represents the points where the denominator vanishes.
In Fig. 1, the evolution of Ω DE , Ω m , Ω r and w eff are plotted versus N for m = 0.25 with the initial conditions x Since in the two cases 0 < c 2 S < 1 along the cosmic evolution, the Laplacian instability of the scalar perturbation is absent. However, a small Laplacian instability appears in the small regime, which is due to the initial conditions we chose. This does not seem to be catastrophic because the duration is short and one can arrange to avoid this instability by changing the values the parameters in the model. We also showed the evolution of β T as a function of N for three values of x (0) 3 . We observe that in three cases β T vanishes in all epochs but it deviates from zero in the late universe by changing the value of x (0) 3 . This means that if our model describes the evolution of our universe, the amplitudes of the gravitational waves have changed in the past.

Observation constraint
In order to better see the deviation of our model from General Relativity and confront with the measurement of observational data, we study the density perturbation of the matter by introducing the growth rate of the matter density where σ 8 (a) = σ 8 δ m (a)/δ m (1) is the rms fluctuation of the linear density field inside a radius of 8h −1 Mpc, and σ 8 represent its present value. The equation of the contrast δ m , in the sub-horizon regime, is given by where Ω m = 1 − Ω r − Ω DE . In the following, we will only focus on the tracker solution. In order to place an observational constraint on the cosmological parameters of the model, we perform a Markov chain Monte Carlo integration implemented with the Metropolis algorithm [54] by using the growth rate factor (GRF) and observational Hubble data (OHD) listed in Table 3. So, the combined likelihood function of the data reads where χ 2 GRF and χ 2 OHD are defined as follows: where It is obvious that χ 2 OHD defined by Eqs. (121) is minimized when H 0 = −B/A. In this case, we write Note that in the following we will use this result to constrain our model instead of that in Eqs.(121). Since our analysis is independent from H 0 , the numerical integration is done over the parameters Ω m 0 , x (0) 3 , m and σ 8 . In Table 2, we show the best-fit parameters with their 1σ confident level for our model and ΛCDM model as well as the Bayes factor. The minimum value for χ 2 of the model studied in this paper is slightly smaller than that in the ΛCDM model where the difference in two models is δχ = 0.9352, which is consistent within 1σ level. However, since the Bayes factor is 1.1 < ln(B 12 ) < 3, the evidence against our model is not strong.
In Fig. 5 we plot the 1σ and 2σ confidence contours in the planes (Ω m 0 , m), (m, σ 8 ), (x  Table 2. We notice that our best-fit values of σ 8 and Ω m 0 are close to the results obtained in the ΛCDM model.  The variation of G eff and w DE is plotted in Fig. 6 for the best-fit values. We observe that effective gravitational coupling starts to evolve at late universe; then it becomes constant in the future. Indeed, G eff starts to decrease at low redshift until it reaches the value 0.74 at N = − 0.30 before going to the value 1.251 in the future of the universe. Moreover, the actual value of G eff and its first derivative appears to be consistent with the constraint obtained from the Solar System [55] and nucleosynthesis constraint [56] observations, which are given by Finally, we compare the growth rate of the two models using the best value presented above with the experiment data Fig. 6 The evolution of G eff and w DE versus N for the best-fit values addressed in Table 2 Fig. 7 The evolution of growth rate of matter fluctuation for our model (Black curve) and ΛCDM model (dashed curve) using the best-fit values addressed in Table 2. We plot also the RSD observations data taken from Table 3 of the Table 3, in Fig. 7. A difference appears between the models during the cosmological evolution where the growth rate in the ΛCDM model is slightly lower than ours from the redshift z = 0.5 to z = 0.

Conclusion
In the present work, we have investigated the cubic galileon in modified teleparalellel gravity frameworks where the Lagrangian is constituted by a general function f (T ), a kinetic therm X and a cubic galileon term X ∇ μ ∇ μ ϕ. Then we calculated the equations of motion in FRW-metric and then the equations corresponding to scalar and tensor perturbations. We found that the expression of the effective gravitational coupling G eff in our model can be a generalization to the cubic galileon and teleparallel gravity. We found also that the gravitation wave speed is equal to the propagation speed of light but the amplitude can evolve in time. The scalar propagation speed sound has been derived for the model proposed in Sect. 2.
In this model, the dynamical analysis showed the existence of a stable point (de Sitter solution) characterized bẏ ϕ = const and H = const. In addition, two interesting behaviors have been found, where the first is as a tracker in which the value of Hφ stays constant with the variation of the cosmic time and the second one is characterized by the conditions x 1 1, x 3 1 and x 2 1 during the matter and radiation dominance. In the tracker solution, the equation of state of dark energy evolves as −1 + (4m)/3 (radiation domination) → m − 1 (matter domination) → −1 (dark energy domination), for x 2 x 3 1 and −7/2 (radiation domination) → −2 (matter domination) → −1 (dark energy domination), for x 3 x 2 1 where the one mentioned last is disfavored by the observational data. However, in the small regime (x 1 1, x 3 1 and x 2 1), the EoS of dark energy behaves like w DE found in the tracker solution for x 2 x 3 1 but we got the same evolution as found in the small regime of the cubic galileon for x 3 x 2 1. In these two regimes, the derived expression of the effective gravitational coupling deviates smoothly from GR in the radiation and matter epochs; then it becomes constant at the de Sitter epoch.
The numerical integration of dynamic equations allowed us to check the evolution of different quantities as a function of N for different initial conditions and values of m. In fact, we showed that the EoS can avoid the value −2 at the matter era when m > 1 and x (0) 3 > 0. Furthermore, in our numerical estimation, it is not only the ghost and Laplacian instability that is absent for the values chosen in Fig. 4, but a variation of the amplitude of GW occurs at late universe.
In Sect. 6, we studied the growth of a matter perturbation by solving the equation corresponding to the contrast δ m numerically in sub-horizon scale. Afterward, we implemented a MCMC numerical integration using the Metropolis algorithm to get the best-fit values of Ω m 0 , x (0) 3 , m and σ 8 . The best-fit values Ω m 0 and σ 8 obtained are consistent with values of Ω m 0 and σ 8 in the ΛCDM model. We also compared the growth rate with experimental data and the prediction of the growth rate associated to the ΛCDM model. We noticed that our model and the ΛCDM model show a smooth difference from each other but the ΛCDM model is not preferred over our model because of the lack of experiment data. Therefore, a combined data analysis of SNIa, CMB, and BAO on our model would be interesting but we leave this to future work.

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