Searching for a Cosmological Preferred Axis in complicated class of cosmological models:Case study $f(R,T)$ model

Recent astronomical observations show that the universe may be anisotropic on large scales. The Union2 SnIa data hint that the universe has a preferred direction. If such a cosmological privileged axis indeed exists, one has to consider an anisotropic expanding Universe instead of the isotropic cosmological model. In this paper, we present a detailed analysis of the dark energy dipole in $f(R,T)=f_{1}(R)+f_{2}(T)$ Cosmological Model using three types of dipole fit (DF) method which are (I)dipole + monopole fitting for distance modulus(DMFDM), (II)dipole + monopole fitting for luminosity distance(DMFLD) and (III) general dipole fitting for luminosity distance(GDFLD). We have found the maximum anisotropic deviation direction for (DMFDM) method as $(l, b)=(315^{+25}_{-25},-23^{+14}_{-15})$, for (DMFLD) as $(l, b)=(l, b)=(315^{+35}_{-37},-23^{+18}_{-18})$, and for (GDFLD) method as $(l, b)=(317^{+32}_{-32},-23^{+18}_{-18})$ which are located very close to each other. We compare our model with the $CPL$, $\Lambda CDM$ and $\omega CDM$ models. Constraints on $(l, b)$ in $f(R, T)$ model are not much different from the cases of the $CPL$, $\Lambda CDM$ and $\omega CDM$ models. Moreover, the results are consistent with other studies.


Introduction
Cosmological principle is one of the basic assumptions of modern cosmology. According to the cosmological principle, the Universe is homogenous and isotropic on scales larger than a few hundred Mpc, which is consistent with currently observational data sets such as the Cosmic Microwave Background (CMB) radiation data from the Wilkinson Microwave Anisotropy Probe (WMAP) [38]- [60]. However, recent observational evidence included Large Scale Velocity Flows (DarkFlow ) [29] anisotropy in the Values of the Fine Structure Constant α (α Dipole) [61,41], anisotropy in Accelerating Expansion Rate (Dark Energy Dipole) [12,19], and other effects [47,5,46] indicate that the Universe may be anisotropic on large scales. A number of authors have investigated the anisotropies of the cosmic acceleration [32,12], which was motivated in several aspects. In particular, several groups such as [51,37] have applied the hemisphere comparison method to study the anisotropy JHEP09(2016)140 of ΛCDM, ωCDM and the dark energy model with CPL parametrization. Some previous works payed attention to study the anisotropic expansion of the universe (Dark Energy Dipole) using the SNIa data and found statistically significant evidence for anisotropies.
More recently, [5] have applied the hemisphere comparison method to the standard ΛCDM model and found that the hemisphere of maximum accelerating expansion is in the direction (l, b) = (309 −23 +23 , 18 −10 +11 ) with Union2 data. [51] took use of the hemisphere comparison method to fit the ΛCDM model to the supernovae data on several pairs of opposite hemispheres, and a statistically significant preferred axis was found. [14] have investigated the anisotropic Cosmological model in the Randers space-time. They found the preferred direction as (l, b) = (306, −18). [10] have taken the deceleration parameter q 0 as the diagnostic to quantify the anisotropy level in the ωCDM model. [11] constructed a direction-dependent dark energy model based on the isotropic background described by the ΛCDM, ωCDM and CPL models and employed the Union2 dataset to constrain the anisotropy direction and strength of modulation. They found the bestfitting value of the maximum deviation direction from the isotropic background is not sensitive to the details of isotropic dark energy models. [59] have studied dipolar anisotropic expansion with cosmographic parameters. They found (l, b) = (309 • , −8.6 • ). [62] chose two simple cosmological models, ΛCDM and ωCDM for the hemisphere comparison approach, and ΛCDM for the dipole fit. In the first approach, they used the matter density and the equation of state of dark energy as the diagnostic qualities in the ΛCDM and ωCDM, respectively. In the second method, they employed distance modulus as the diagnostic quality in ΛCDM. They found a preferred direction of (l, b) = (307 • , −14 • ).
In testing for anisotropy or consistency with isotropy, we can ask which cosmological probes are most sensitive in what redshift ranges to such a hypothetical anisotropy, i.e. what constraints could be put on angular variations in the local dark energy equation of state.
We cannot make a convincing conclusion from only one dataset, model or method about the origin of the anisotropy. Anisotropy may come from systematic uncertainty, as well as the intrinsic property of the universe. If the privileged axes derived from different datasets, different methods and different cosmological models are close to each other, we can safely conclude that anisotropy is an intrinsic property of the Universe. As we mentioned above several studies payed attention to find a preferred axis of the Universe in isotropic background described by the ΛCDM, ωCDM and CPL models; however, possibility of existence a privileged axis for the Universe may enhance if modified cosmological models such as f(R,T) could predict it and produce closely result. Searching for preferred cosmological axis of the Universe, we focus on generalized gravity model f (R, T ).
Recently, the observations of high redshift type Ia supernovae, the surveys of clusters of galaxies [56,49], Sloan digital sky survey (SDSS) [1] and Chandra X-ray observatory [2] reveal the universe accelerating expansion and that the density of matter is very much less than the critical density. Also, the observations of Cosmic Microwave Background (CMB) anisotropies [7] indicate that the universe is flat and the total energy density is very close to the critical one [53]. The observations though determines basic cosmological parameters with high precisions and strongly indicates that the universe presently is dominated by a JHEP09(2016)140 smoothly distributed and slowly varying dark energy (DE) component, but at the same time they poses a serious problem about the origin of DE [57]. The most The 'cosmological constant' as the best candidate for explaining cosmic acceleration in literature faces serious problems such as fine-tuning and a huge discrepancy between theory and observations [20]- [58]. On the other hand, modification of the geometrical part of the Einstein-Hilbert action by replacing an arbitrary function of the Ricci scalar R [45] has constructed welldeveloped dark energy models. This phenomenological approach is called as the Modified Gravity. Using the Modified Gravity we can strongly explain the rotation curves of galaxies, the motion of galaxy clusters, the Bullet Cluster, and cosmological observations without the use of dark matter or Einstein's cosmological constant [44,30]. Cosmic inflation, mimic behavior of dark matter and current cosmic acceleration being compatible with the observational data are other successful predications of the f (R) theories [44]- [54].
A generalization of f (R) modified theories of gravity was proposed in [6] studies, by including in the theory an explicit coupling of an arbitrary function of the Ricci scalar R with the matter Lagrangian density L m . A specific application of the latter f (R, L m ) gravity was proposed in [48] studies, which may be considered a relativistically covariant model of interacting dark energy, based on the principle of least action. The cosmological constant in the gravitational Lagrangian is a function of the trace of the stress-energy tensor, and consequently the model was denoted "Λ(T ) gravity". It was argued that recent cosmological data favor a variable cosmological constant, which are consistent with Λ(T ) gravity, without the need to specify an exact form of the function Λ(T ) [48]. Λ(T ) gravity is more general than the Palatini f (R) gravity, and reduces to the latter when we neglect the pressure of the matter.
In this paper, a class of the Modified Gravity theories in which the gravitational action contains a general function f (R, T ), where R denotes the Ricci scalar and T is the trace of the energy-momentum tensor, has been considered. [33] introduced this type of the Modified Gravity, f (R, T ), which obtained significant outcomes: the reconstruction of cosmological solutions, where late-time acceleration was accomplished by [34] and the energy conditions was analyzed by [3]. [52] studied the thermodynamics of Friedmann-Lemaître-Robertson-Walker (FLRW) spacetimes. Moreover, the occurrence possibility of future singularities was studied by [35]. Besides these achievements, a serious shortcoming in this kind of theory has been the non-conservation of the energy-momentum tensor. To circumvent this problem, in this paper, we show that f (R, T ) functions can always be constructed in a way to be consistent with the energy-momentum tensor standard conservation. In this regard, we can assume separable algebraic functions of the form f (R, T ) = f 1 (R) + f 2 (T ) in which the function f 2 (T ) is obtained by imposing the conservation of the energy-momentum tensor. Now, in order to search for dipolar asymmetry, we construct an anisotropic dark energy model and aim to detect the maximum anisotropy direction. Furthermore, we consider the impact of redshift on the direction by using the redshift tomography method, with the Union2 data. Finally, we compare our results for the f(R,T) model with ΛCDM, ωCDM and CPL models and also some previous studies.
The paper is structured as follows: in section 2, we obtain the field equations of f (R, T ) gravity and analyze the stability of the dynamical system of the f(R,T) model. In JHEP09(2016)140 section 3, we discuss free parameters of the model in some detail and constrain these free parameters using observational data. In section 4, we investigate the scalar perturbations in f (R, T ) = f 1 (R) + f 2 (T ) type theories. In order to search for Dark Energy Dipole in the model using observational data, in section 5 we describe some important anisotropy models and method. Then we introduce and extend types of Dipole-Fitting method in order to investigate possible anisotropy from the data. We compare the results of these three types of DF method used for the f(R,T) model with each other in section 6. Moreover, in order to explore the possible redshift dependence of the anisotropy, we have implemented a redshift tomography analysis in section 7. In section 8, we have performed the anisotropy analysis for CPL, ΛCDM and ωCDM models. We have applied DMFLD method to find the anisotropy of ΛCDM and ωCDM and the dark energy model with CPL parametrization in order to make a comparison between these models and the f (R, T ) model. Finally, in section 9, we conclude, summarize and compare the results of this work with some of the recent studies of [41,61,11], searching for evidence for a preferred cosmological axis.
2 Field equations of f (R, T ) model The action of f (R, T ) gravity is of the form is an arbitrary function of the Ricci scalar and T (m) , L (m) and L (rad) are the Lagrangian of the dust matter and radiation, g is the determinant of the metric, µν is the trace of the energy-momentum tensor and we set c = 1. By varying the action (2.1), with respect to the metric tensor g µν , the field equations are obtained as where ∇ denotes the covariant derivative and Contracting of equation (2.2) yields Now, in this model, we assume the perfect fluid and the spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric

JHEP09(2016)140
where a(t) is the scale factor. Let us rewrite (2.2) as a standard form similar to GR, i.e.
Regarding the Bianchi identity, obviously in f (R, T ) gravity, the above effective energymomentum tensor is not conserved. Thus, by applying the conservation of the energymomentum tensor of the whole matter, i.e. ∇ µ T (m) µν , the following constraint must be held. That is as the Friedmann-like equation, and as the Raychaudhuri-like equation. [33] gave three classes of these models In this paper, we have focused on the second class f (R, T ) = f 1 (R) + f 2 (T ). Now, by rewriting equations (2.9) and (2.10) for this model, one can obtain For more simplicity, we introduce a few independent new variables as where the prime denotes differentiating with respect to the argument and we have used R = 6(Ḣ + 2H 2 ) for metric (2.5). By applying the conservation equation (2.8) for the minimal combination, f (R, This constraint restricts its form to a particular one, namely, where c 1 and c 2 are constants with respect to T. The conservation of the energy-momentum tensor also leads to the case in which the variable ν is a function of ξ, namely, ν = ξ 2 . Therefore, these six variables will reduce to five independent variables once the constraint equation (2.8) is applied. 20) Where N represents derivatives with respect to ln a and α ≡ where C 1R and C 2R are constants. Note that the second order nonlinear differential equations of systems was simplified to first order differential equations by introducing a few new variables. It is interesting to consider the behavior of systems around the equilibrium points ( dζ dN = 0 dη dN = 0 dϑ dN = 0 dξ dN = 0 dχ dN = 0) using jacobian stability analysis. The Jacobin stability of a dynamical system can be regarded as the robustness of the system to small perturbations of the whole trajectory. This is a very convenable way of regarding the resistance of limit cycles to small perturbation of trajectories. It gives us the possibility of study all the evolutional paths admissible for all initial conditions [24]- [27]. It

JHEP09(2016)140
Fixed points Coordinates (ζ, η, ϑ, ξ, χ) eigenvalue stability P1 (0, 0, 0, 0, 1) (1, 0, 0, 0, 0) 1, 2, 5, 7 2 , 4(α−1) Table 1. The fixed points solutions of the dynamical system problem of f (R, is especially important in cosmology where there is the problem of initial conditions. Using the dynamical systems methods one hopes to answer the question of what is the range of initial conditions and parameters of the system for which the subsequent evolution is compatible with current Evaluating the Jacobian matrix at the steady state and computing the corresponding eigenvalues of them, we can investigate stability or instability based on the real parts of the eigenvalues. Table 1 shows the property of critical points of dynamical system.

Constrain on parameters of the model
In pervious section, we investigated stability of dynamical system by introducing the dimensionless parameters {ζ, η, ϑ, ξ, χ, ν}, it is obvious that the critical points and eigenvalues dependent only on the free parameter α (see table 1). Also, we can see from equation (2.16) to (2.20) that parameter α is the only parameter which has been explicitly revealed in the set of equations and has directly affected the dynamical system; however, there are some parameters such as {C 1R , C 1T , C 2R , C 2T } which have not appeared in the set of equations (2.16) to (2.20) and have been masked by the dimensionless new variables; but they can affect on dynamics of the system. In fact, the new variables dependent on these parameters. For example, we can rewrite the variable η as We can see that this variable dependents on C 1R , C 2R and α. It is important to note that in a dynamical system with a set of equations both free parameters and initial conditions determine the dynamics of the system. Free parameters affect critical points and initial conditions affect the trajectories of variables in phase space. Here, the parameter α is the only free parameter. Although the parameters {C 1R , C 1T , C 2R , C 2T } have not appeared in the equations explicitly, they can affect value of initial conditions. In order to study the effect of these parameters and constrain them JHEP09(2016)140 with observation, we reveal these parameters in new set of equations. In this respect, we introduce some other new variables as Equation (2.9), (2.15) and (2.21) give the following constraints betweenḟ R (R, T ) and the variables x 1 to x 5 , namelẏ Now,using (3.2) and (3.3), for the autonomous equations of motions, we obtain Therefore, we have a dynamical system with four independent variables and five free parameters α, C 1R , C 2R , C 1T and C 2T . Note that this system of equations is corresponding to equations (2.16) to (2.20). We have best-fitted these parameters using SNe Ia data by χ 2 method. The likelihood for these parameters have been shown in figure 3.

type theories
Let us consider the scalar perturbations of a flat FRW metric in the longitudinal gauge and in conformal time: where Φ ≡ Φ(η, x) and Ψ ≡ Ψ(η, x) are the scalar perturbations. The components of perturbed energy-momentum tensor in this gauge are given bŷ where v denotes the potential for the velocity perturbations. The first order perturbed equations in a dust matter dominated universe, c 2 s = 0 will be obtained as [4] Φ where κ 2 = 8πG, the prime holds for the derivative with respect to η, H ≡ a /a and the subscript 0 holds for unperturbed background quantities: R 0 denotes the scalar curvature corresponding to the unperturbed metric, ρ 0 the unperturbed energy density, with

Solution of the equations using dynamical system
The complete set of equations that describes the general linear perturbations for the model have been presented in pervious section. These equations are a set of nonlinear second order differential equations with a large number of variable for which there is no analytical JHEP09(2016)140 solution except for simplest cases and only numerical analysis can be performed. Our purpose is to convert second order differential equation to first order by introducing some new variables. There are various reasons for doing this, one being that a first order system is much easier to solve numerically. Also, it allows us to investigate the behavior of the system in phase space. Phase planes are useful in visualizing the behavior of the system particularly in oscillatory systems where the phase paths can "spiral in" towards zero, "spiral out" towards infinity, or reach neutrally stable situations called centres. This is a useful method to determine whether dynamics of a system are stable or not.
The structure of phase space of the field equations is simplified by defining a few variables and parameters. These variables are generally defined as Now, for the autonomous equations of motions, we obtain

JHEP09(2016)140
CriticalP oints χ 1 After some calculation from equations (4.3)-(4.9), we can obtain the above parameters in terms of the new variables as By substituting equations (4.31)-(4.35) into equations (4.18)-(4.25), the complete set of equations that describes the behavior of the system in terms of new variables will be provided.
In general, the critical points and eigenvalues of the system will be obtained in terms of β, k. Here, we have obtained critical points of the system for β = 1, k = 0.3 (see table 2). The corresponding eigenvalues are as:

JHEP09(2016)140
Due to the fact that there are complex values in some matrix elements of the eigenvalues, the dynamical system shows attractor behavior. The attractor behavior of the system has been shown in figure 4. Note that the attractor behavior in phase space implies that the system oscillates and moves toward steady state in a critical point.  The oscillating behavior of the system has been shown in figure 5. Moreover, we are interested in behavior of the parameters δ, Ψ, Φ. Therefore, we can reconstruct them from new variables as Figure 6 shows the oscillating behavior of the parameters δ, Φ, Ψ.

Solution for
In this section, we solve the equations for the → 0. Note that this limit in f (R) theories corresponds to the large scalaron mass which defined as [23,43,28]; Applying this condition to equation (4.3) yields Ψ = Φ. Therefore, the equations (4.3)-(4.9) are simplified as follows From equation (4.41) we have Hence, the autonomous Equation of Motion for the independent variables can be obtained via Also, from equation (4.46) we obtain Applying constraint (4.53), the system reduces to a system with five independent variables.

JHEP09(2016)140
By setting x 6 = 0, behavior of the dynamical system in f (R) theory will be provided. Here, we have plotted two dimensional, three dimensional phase space and evolution of variables for f (R) theory in figure 8.
Attractor property and oscillating behavior of the dynamical system in figure 9 shows that the trajectories spirals out from an unstable focus point and moves towards a steady state point.

Searching for dark energy dipole using observational data
There are various ways to investigate possible anisotropy from the data. Generally speaking, there are three important ways.

Modification of the luminosity distance redshift relation in a specific anisotropic cosmological model
In this method, an expression is derived for the luminosity distance as a function of redshift in a specific anisotropic cosmological model. Many anisotropic cosmological models with modified luminosity distances have been proposed to match observations. Table 3 shows modified luminosity distance for some of these models as an incomplete list.
JHEP09(2016)140 Figure 9. Attractor property and oscillating behavior of the dynamical system for f (R) theory when The Bianchi I type cosmological model [12,50] and the Rinders-Finsler cosmological model [14,15] are two models which are consistent with the SNe Ia data. A scalar perturbation of the ΛCDM model may also break the spherical symmetry of the Universe such that a preferred axis arises. [13,39,59]  ) [39] [59] 2 Anisotropic d L in the Finslerian space-time [9] 4 "wind" scenario to the bulk flow (1−e 2 ) 1/6 (1−e 2 cos θ) 1/2 [12] in ellipsoidal universe Here the luminosity distance is expressed in terms of the observed redshift z, and the direction to the source, where n denotes a unit spatial vector from the observer to the source. The notation ≈ denotes an approximate equality accurate up to first-order in the potentials Φ, Ψ, and the peculiar velocities (in the conformal spacetime) v. The subscripts s and o refer to the source and the observer, respectively. The affine parameter χ is given by This relation is appropriate for investigation of the peculiar velocities. It should be mentioned that the aim of our study is to investigate Dark Energy Dipole in the model.

Hemisphere Comparison (HC) method
The HC method divides the data points into two subsets according to their position in the sky and fit the subsets to an isotropic cosmological model (e.g., ΛCDM model). Several groups such as [51,37] = g(z)(ẑ.n = g(z) cos θ) [10]  recently, [5] have applied the hemisphere comparison method to the standard ΛCDM model and found that the hemisphere of maximum accelerating expansion is in the direction of (l, b) = (309 −3 +23 , 18 −10 +11 ). [51] took use of the hemisphere comparison method to fit the ΛCDM model to the supernovas data on several pairs of opposite hemispheres, and a statistically significant preferred axis was found. Some of studies which have used this method are listed in table 4.

Dipole-Fitting (DF) method
Using the DF method, we can directly fit the data to a dipole (or dipole plus monopole) model. If the Universe is really intrinsically anisotropic and there exists a preferred direction, it should directly affect the expansion rate of the Universe, leading to the anisotropic luminosity distance and anisotropic distance modulus. In fact, this method corresponds to the fluctuation of the distance modulus. Anisotropic Dipole-fitting method has been used for searching the anisotropy of fine structure constant using quasars data on cosmological scale. [5] firstly applied this method to anisotropic study using SNe Ia dataset. [62] have applied this method to investigate dipolar asymmetry of the Universe. [17] have made a comprehensive comparison between the HC method and the DF method using the Union2 dataset.
Several studies payed attention to the fluctuation of the distance modulus in order to find the preferred axis of the Universe using the DF method (See table 4). It is worth to mention that the anisotropic property of the Universe directly affect the luminosity distance and leading to the anisotropic luminosity distance. Therefore, besides the DF method for the distance modulus, we have used this method for Luminosity distance. We have explained three types of Dipole-fitting method which are based on deviation of distance modulus and Luminosity distance from their best values in isotropic model to investigate the anisotropic expansion of the Universe.

JHEP09(2016)140
(II) Calculate the angle of each supernova with respect to the dipole axis, which is determined by cos θ i = z i . n (5.2) where z i is the unit direction vector of the supernova, which can be expressed by using the Galactic coordinate system.
and n is the direction of dark energy dipole, which is the maximal expanding direction, where (l, b) is the Galactic coordinate direction of dipole axis.
(III) Define the angular distribution model with dipole and monopole where m 1 and d 1 denote the monopole and dipole magnitude, respectively,μ is the distance modulus predicted by the isotropic f (R, T ) model and µ is the true luminosity distance of the supernova.
(IV) Fit the SNIa data by minimizing the χ 2 sn value of the distance modulus. The χ 2 sn for SNIa is obtained by comparing theoretical distance modulus with observed µ obs of supernovae. We suppose the experiment error between each measurement is completely independent, so the covariance matrix can be simplified as the diagonal component, and the χ 2 sn can be written as where µ th (z i ) is the theoretical distance modulus which it will be obtain from equation (5.5) as andμ(z i ) = 5 log 10 [d L (z)] + 42.38 − 5 log 10 h also for a FRW cosmological model, one hasd dz .  In this step, we employ the Union2 dataset to constrain the anisotropic dark energy model. The directions of the SNIa that we have used here are given in [8] work, and are described in the equatorial coordinates (right ascension and declination). In order to make comparisons with other results, we convert these coordinates to the galactic coordinates (l, b) [22].
The parameters need to be constrained are (d 1 , m 1 , l, b). Using the least χ 2 sn method, we can find the best-fit parameters (d 1 , m 1 , l, b). The best-fit dipole direction is found to be towards (l, b) = (315 0 ± 25 0 , −23 0 ± 15 0 ) (5.11) Figure 10 shows the two dimensional likelihood for parameters (l, b) in f(R,T) model using DMFDM method. The distribution of Union2 SnIa Datapoints in galactic coordinates along with the dark energy dipole direction (l, b) are shown in figure 11. The magnitude of the dipole and the monopole have obtained as We can see that the magnitude of the monopole is one order of magnitude smaller than that of the dipole. This is consistent to the result of [41], who obtained the result of [17], who fitted the data with a dipole only and obtained  We have obtained the likelihood function of each parameter by performing the χ 2 analysis using 10 5 data point. The results are shown in figure 12 and figure 13.
B. Dipole+Monopole Fitting for Luminosity Distance (DMFLD). We perform a similar dipole+monopole fit using the Union2 data. Instead of ( ∆µ(z) µ(z) ) which corresponds to distance modulus deviations from its isotropic f (R, T ) value, we use the luminosity distance deviation from its best fit isotropic f (R, T ) value where,d L (z) is the luminosity distance of the supernova in isotropic background and d L (z) is the true luminosity distance or anisotropic luminosity distance of the supernova. Therfore we can use the following expression We have obtained the likelihood function of each parameter by performing the χ 2 analysis using 10 5 data point. The results are shown in figure 16 and figure 17.
We have also obtained (m 10 −4 ). Therefore, neglecting m and by considering dipole magnitude as a function of z, the general case of Luminosity Distance dipole fit will be

Comparison of three DF models
In the previous section, we have described three types of dipole-fitting (DF) method which has been used for statistical analysis in order to find the preferred cosmological axis of the Universe in f (R, T ) model. These three types are as follows:  method to study the anisotropy of ΛCDM, ωCDM and the dark energy model with CPL parametrization. We have applied all of these DF methods to study privilege axis of the universe in f (R, T ) model. At first, it seems that these methods have a same origin (because of the direct relation between µ and d L ). Also, the best fitted direction of preferred axis of these methods are very close to each other. In fact, DMFDM ((l, b) = (315 0 ± 25 0 , −23 0 ± 15 0 )) and DMFLD ((l, b) = (315 0 ± 37 0 , −23 0 ± 18 0 )) methods have resulted exactly the same value for the privilege axis of the universe in f (R, T ) model. However, their 1 − σ confidence level are different. As left panel of figure 16 shows, the (1 − σ) confidence region of DMFDM is smaller than DMFLD (right panel of figure 22). Moreover, they give different values of dipole magnitude which is interesting to note. The dipole magnitude obtained using DMFDM method (d 1 = (1.4 ± 0.8) × ×10 −3 ) is close to previous studies of [17,59,62]    has been mentioned in DMFDM method section. However, the dipole magnitude obtained using DMFLD method (d 2 = (0.026 ± 0.014)) is different from the value obtained using DMFDM method and also previous studies. Interestingly, the magnitude of anisotropy (d 2 = (0.026 ± 0.014)) obtained using DMFLD method is approximately equal to that of CMB dipole. The recent released Planck data show that the dipole magnitude of CMB temperature fluctuations is about A = 0.07-0.01 [15].
There are two reasons to study the dark energy dipole of the universe using the formula based on deviation on Luminosity distance (DMFLD method) instead of distance modulus (DMFDM method). The first is that if dark energy has anisotropic repulsive force, it will directly affect the expansion rate of the Universe, leading to the anisotropic luminosity distance; therefore, in formulating the dipole-fitting method it is more appropriate that the d L be revealed directly in the equation. The later reason is that most of the formulaes for modification of d L presented in table 3 with very small values of dipole magnitude  (d  1) can be simplified as which is the same as the DF equation of DMFLD and GDFLD methods.    Table 6. Constraints of the directions and amplitude of maximum anisotropy using DMFLD method for different redshift bins of the SNIa data. The error-bars quoted is 1σ error. In this section, we compare c R R α+1 + c T √ −T Gravity model with CPL parametrization, ωCDM and ΛCDM models. In the framework of a spatially flat Friedmann universe, the expansion history of the Universe is given by where H =ȧ a is the Hubble parameter, q is the deceleration parameter, Ω m0 = ρ 0 ρc is the current value of the normalized matter density, Ω x (z) is the normalized dark energy density  Table 7. Constraints of the directions and amplitude of maximum anisotropy using GDFLD method for different redshift bins of the SNIa data. The error-bars quoted is 1σ error.
as a function of redshift which evolves as Ω x (z) = Ω x0 f (z) Next, we turn to the parametrization of w(z). There are many functional forms of w(z) in the literature. In this work, we consider Chevallier-Polarski-Linder (CPL) parametrization introduced by [18,40], which invokes as barotropic factor the known expression w(z) = w 0 + w 1 z 1 + z (8.4) In this case, the equation of state becomes w(z = 0) = w 0 at present time and w(z → ∞) = w 0 + w 1 at earlier time. This simple parametrization is most useful if dark energy is important at late times and insignificant at early times. In addition to its simplicity, this CPL parametrization exhibits interesting properties. However, it cannot describe rapid variations in the equation of state. Using this functional form and equation (8.1), equation (8.3) can be written analytically as While in the case of ωCDM model, the equation of state of dark energy is parameterized by a constant ω = p ρ ; therefore, we have Using Union2 data and by χ 2 method, we have best fitted parameters Ω m0 , Ω Λ , h for ΛCDM and parameters Ω m0 , Ω CPL , h, ω 0 , ω 1 for CPL model. Figure 26 shows the confidence levels for parameters (Ω m0 , h) in both ΛCDM and CPL models.    We have also considered the ωCDM, ΛCDM and the CPL parameterized dark energy models as the isotropic background. We use the isotropic background dark energy parameters in table 8 and fit our anisotropic parameters, respectively. The results are summarized in table 9. Figure 27, figure 28 and figure 29 show the results of constraints on (l, b) which are not much different from the case of the f (R, T ) model. This means that the best-fitting value of the maximum deviation direction from the isotropic background is not sensitive to the details of isotropic dark energy models. The best fitted trajectories of the effective EoS parameter in isotropic, anisotropic c R R α+1 + c T √ −T gravity, CPL and ωCDM and ΛCDM models are shown in figure 30. Based on this, the trajectory of anisotropy is not in much difference from the case of the isotropy, also the best fitted trajectory of CPL and ΛCDM models are same at late time and different in the future and ωCDM, ΛCDM and c R R α+1 + c T √ −T gravity have same trajectory at late time and in the future.

Conclusions
We have investigated the cosmological solutions of f (R, T ) gravity, in isotropic and anisotropic space-time. In both isotropic and anisotropic cases, our studies are based on the phase-space analysis (the dynamical system approach). In this approach, we convert a set of second order differential equations to a new set of first order ones by defining some -- Table 9. The preferred direction of the universe in ΛCDM, CPL and ωCDM models.  dimensionless variables and parameters. There are various reasons for doing this: a first order system is much easier to solve numerically, and also phase planes are useful in visualizing the behavior of dynamical systems, especially in oscillatory systems where the phase paths can "spiral in" towards zero, and "spiral out" towards infinity. Moreover, it gives us useful information about (in)stability of the system and critical points of the system. At first, we have obtained the field equations of f (R, T ) gravity in isotropy case and have analyzed the stability of the dynamical system for f (R, T ) = f 1 (R) + f 2 (T ). Then, JHEP09(2016)140 we have studied the evolution of scalar cosmological perturbations in the metric formalism. The main purpose of scalar perturbations is to find explicit expressions for the parameter Φ, Ψ and δ in the framework of nonlinear f (R, T ) model. Unfortunately, the system of equations for scalar perturbations is very complicated in the case of nonlinearity. It is hardly possible to solve it directly. Therefore, we have used phase space approach to simplify the nonlinear equations of the f (R, T ) = f 1 (R) + f 2 (T ) model. We have also reconstructed the parameters Φ, Ψ and δ from new variables. In the model, the evolution of matter density perturbations for different cases has been studied and the corresponding results have been shown in figure 4, 5, 6, 7, 8 and 9. The attractor property (spiral in and out) of the system leads to an oscillating behavior of the matter perturbations and other variables and parameters. This behavior is predictable for the critical points of the dynamical systems whose eigenvalues are complex.
In section 5, we have investigated Dark Energy Dipole in the f (R, T ) model using Dipole Fitting method. There is a range of independent cosmological observations which indicate the existence of anisotropy axes. This appears to be one of the most likely direc-JHEP09(2016)140  tions which may lead to new fundamental physics in the coming years. These cosmological observations along with their preferred directions and the corresponding references are summarized in table 11.
In this paper, we present a detailed analysis of the dark energy dipoles in f (R, T ) = f 1 (R) + f 2 (T ) cosmological model using three types of dipole-fitting (DF) method which are (I)dipole + monopole fitting for distance modulus (DMFDM), (II)dipole + monopole  fitting for luminosity distance (DMFLD) and (III) general dipole fitting for luminosity distance (GDFLD).
Several groups have applied DMFDM method to study the anisotropy of ΛCDM, ωCDM and the dark energy model with CPL parametrization. [11] have applied GDFLD method to study the anisotropy of ΛCDM, ωCDM and the dark energy model with CPL parametrization. We have applied all of these DF methods to study privilege axis of the universe in f (R, T ) model. At first, it seems that these methods have a same origin (because of the direct relation between µ and d L ). Also, the best fitted direction of preferred axis of these methods are very close to each other. In fact, DMFDM ((l, b) = (315 0 ± 25 0 , −23 0 ± 15 0 )) and DMFLD ((l, b) = (315 0 ± 37 0 , −23 0 ± 18 0 )) methods have resulted exactly the same value for the privilege axis of the universe in f (R, T ) model. However, their 1 − σ confidence level are different ( figure 22). The (1 − σ) confidence region of DMFDM is smaller than DMFLD. Moreover, they give different values of dipole magnitude which is interesting to note. The dipole magnitude obtained using DMFDM method (d 1 = (1.4 ± 0.8) × ×10 −3 ) is close to previous studies of [17,59,62] as it has been mentioned in DMFDM method section. However, the dipole magnitude obtained using DMFLD method (d 2 = (0.026 ± 0.014)) is different from the value obtained using DMFDM method and also previous studies. Interestingly, the magnitude of anisotropy (d 2 = (0.026 ± 0.014)) obtained using DMFLD method is approximately equal to that of CMB dipole. The recent released Planck data show that the dipole magnitude of CMB temperature fluctuations is about A = 0.07-0.01 [15]. Also, it is close to the result of [15] which have obtained the magnitude of dipolar asymmetry as |D| = 0.044±0.018, using modified luminosity distance in anisotropic cosmological model in the Finsler-Randers spacetime (formula 2 of table 3). ). The results for preferred direction in other models are presented for contrast.Point × denotes the result of [59], point denotes the result of [62],point denotes the result of [11], point denotes the result of [16], point denotes the result of [16], point denotes the result of [14], point denotes the result of [61], point • denotes the result of [14], point denotes the result of [41], and point + denotes the result of [11]. The light green represents the 1-σ errors on the Dark Energy dipole direction, which includes the results for preferred direction in other models.