Bouncing cosmological solutions from f(R,T) gravity

In this work we study classical bouncing solutions in the context of $f({\sf R},{\sf T})={\sf R}+h({\sf T})$ gravity in a flat {\sf FLRW} background using a perfect fluid as the only matter content. Our investigation is based on introducing an effective fluid through defining effective energy density and pressure; we call this reformulation as the"effective picture". These definitions have been already introduced to study the energy conditions in $f({\sf R},{\sf T})$ gravity. We examine various models to which different effective equations of state, corresponding to different $h({\sf T})$ functions, can be attributed. It is also discussed that one can link between an assumed $f({\sf R},{\sf T})$ model in the effective picture and the theories with generalized equation of state ({\sf EoS}). We obtain cosmological scenarios exhibiting a nonsingular bounce before and after which the Universe lives within a de-Sitter phase. We then proceed to find general solutions for matter bounce and investigate their properties. We show that the properties of bouncing solution in the effective picture of $f({\sf R},{\sf T})$ gravity are as follows: for a specific form of the $f({\sf R,T})$ function, these solutions are without any future singularities. Moreover, stability analysis of the nonsingular solutions through matter density perturbations revealed that except two of the models, the parameters of scalar-type perturbations for the other ones have a slight transient fluctuation around the bounce point and damp to zero or a finite value at late times. Hence these bouncing solutions are stable against scalar-type perturbations. It is possible that all energy conditions be respected by the real perfect fluid, however, the null and the strong energy conditions can be violated by the effective fluid near the bounce event.


Introduction
Today, the standard cosmological model (SCM) or the Big-Bang cosmology has become the most acceptable model which encompasses our knowledge of the Universe as a whole. For this reason it is called also the "concordance model" [1]. This model which allows one to track the cosmological evolution of the Universe very well, has matured over the last century, consolidating its theoretical foundations with increasingly accurate observations. We can numerate a number of the successes of the SCM at the classic level. For example, it accounts for the expansion of the Universe (Hubble law), the black body nature of cosmic microwave background (CMB) within the framework of the SCM can be understood and the predictions of light-element abundances which were produced during the nucleosynthesis. It also provides a framework to study the cosmic structure formation [2]. However, though the SCM works very well in fitting many observations, it includes a number of deficiencies and weaknesses. For instance some problems which are rooted in cosmological relics such as magnetic monopoles [3], gravitons [4], moduli [5] and baryon asymmetry [6]. Despite the self-consistency and remarkable success of the SCM in describing the evolution of the Universe back to only one hundredth of a second, a number of unanswered questions remain regarding the initial state of the Universe, such as flatness and horizon problems [7]. Moreover, there are some unresolved problems related to the origin and nature of dark matter (DM) [8]. Notwithstanding the excellent agreement with the observational data there still exists a number of challenging open problems associated with the late time evolution of the Universe, namely the nature of dark energy (DE) and cosmological constant problem [9]. Though the inflation mechanism has been introduced to treat some of the mentioned issues such as, the horizon, flatness and magnetic monopole problems at early Universe [10], the SCM suffers from a more fundamental issue, i.e., the initial cosmological singularity that the existence of which has been predicted by the pioneering works of Hawking, Penrose and Geroch in 1960s, known as the singularity theorems [11] and their later extensions by Tipler in 1978 [12] and by Borde, Vilenkin and Guth in the 1990s [13] (see also [14] for a comprehensive study). According to these theorems, a cosmological singularity is unavoidable if spacetime dynamics is described by General Relativity (GR) and if matter content of the Universe obeys certain energy conditions. A singular state is an extreme situation with infinite values of physical quantities, like temperature, energy density, and the spacetime curvature from which the Universe has started its evolution at a finite past. The existence of such an uncontrollable initial state is irritating, since "a singularity can be naturally considered as a source of lawlessness" [15]. A potential solution to the issue of cosmological singularity can be provided by "non-singular bouncing cosmologies" [16]. Beside a huge interest in the solutions that do not display singular behavior, there can be more motivations to seek for non-singular cosmological models. The first reason for removing the initial singularity is rooted in the initial value problem since a consistent gravitational theory requires a well-posed Cauchy problem [17]. However, owing to the fact that the gravitational field diverges at a spacetime singularity, we could not have a well formulated Cauchy problem as we cannot set the initial values at a singular spatial hypersurface given by t = const. Another related issue is that the existence of a singularity is inconsistent with the entropy bound S/E = (2πR)/c , where S, E, R, and c being entropy, proper energy, the largest linear dimension, Planck's constant and the velocity of light, respectively [15].
During the past decades, models which describe bouncing behavior have been designed and studied as an approach to resolve the problem of initial singularity. These models suggest that the Universe existed even before the Big-Bang and underwent an accelerated contraction phase towards reaching a non-vanishing minimum radius. The transition from a preceding cosmic contraction regime to the current accelerating expansion phase (as already predicted in SCM) is the so called "Big Bounce". From this perspective, the idea that the expansion phase is preceded by a contraction phase paves a new way towards modeling the early Universe and thus, may provide a suitable setting to obviate some of the problems of the SCM without the need to an inflationary scenario. Although an acceptable model can be considered as the one being capable of explaining the issues that have been treated by inflationary mechanism, e.g., most inflationary scenarios can give the scale-invariant spectrum of the cosmological perturbations [18], problems of the SCM may find solutions in the contracting regime before the bounce occurs. The horizon problem, for example, is immediately resolved if the far separated regions of the present Universe were in causal connection during the previous contraction phase. Similarly, the homogeneity, flatness, and isotropy of the Universe may also be addressed by having a smoothing mechanism in the contraction phase, see e.g., [19] for more details. Moreover, though the fine-tuning is required to keep a stable contracting regime, the nonsingular bounce succeeds in sustaining a nearly scale-invariant power spectrum [20].
For several years, great effort has been devoted to the study of bouncing cosmologies within different frameworks. The resulted cosmological models could be obtained at a classical level or by quantum modifications. Most of the efforts in quantum gravity are devoted to reveal the nature of the initial singularity of the Universe and to better understand the origin of matter, non-gravitational fields, and the very nature of the spacetime. Non-singular bouncing solutions generically appear in loop quantum cosmology (LQC) [21], where the variables and quantization techniques of loop quantum gravity are employed to investigate the effects of quantum gravity in cosmological spacetimes [22,23,24]. The recent large amount of works done within the loop quantum cosmology (LQC) show that when the curvature of spacetime reaches the Planck scale, the Big-Bang singularity is replaced by a quantum Big-Bounce with finite density and spacetime curvature [25]. Another approach, based on the de Broglie-Bohm quantum theory, utilizes the wave function of the Universe in order to determine a quantum trajectory of the Universe through a bounce [26]. In the framework of LQC the semi-classical Friedmann equations receive corrections as [24,27] where, ρ max ≈ 0.41ρ Pl and ρ Pl = c 5 / G 2 being the Planck density [24,28]. We note that, the relative magnitude of ρ and ρ max enables one to distinguish classical and quantum regimes. By a short qualitative inspection of the above equations the general feature of the bouncing behavior will be revealed: initially, the Universe were in contracting phase at which the matter density and curvature are very low compared to the Planck scale. As the Universe contracts more, the maximum density is reached so that the quantum evolution follows the classical trajectory at low densities and curvatures but undergoes a quantum bounce at matter density ρ = ρ max , where we have H LQC = 0 and alsoḢ LQC = 4πG(ρ + p). The quantum regime then joins on to the classical trajectory that was contracting to the future. Therefore, the quantum gravity effects create a non-singular transition from contraction to expansion and thus the Big-Bang singularity is replaced by a quantum bounce. Furthermore, we see that for all matter fields which satisfy the weak energy condition (WEC) we haveḢ LQC > 0. These two results are accounted for the general conditions for the existence of a bouncing solution. Moreover, nonsingular bouncing scenarios have been also reported in nonlocal gravity model [29] where effective Friedmann equation with a slight difference to equation (1) has been proposed exhibiting a bouncing-accelerating behavior. See also [30] for probing the issue of singularity avoidance in nonlocal gravitational theories. Another type of theories are called non-singular "matter bounce" scenarios which is a cosmological model with an initial state of matter-dominated contraction and a non-singular bounce [31]. Such a model provides an alternative to inflationary cosmology for generating the observed spectrum of cosmological fluctuations [19,32,33]. In these theories some matter fields are introduced in such a way that the WEC is violated in order to makeḢ > 0 at the bounce. From equation (2), it is obvious that putting aside the correction term leads to negative values for the time derivative of the Hubble parameter for all fluids which respect WEC. Therefore, in order to obtain a bouncing cosmology it is necessary to either go beyond the GR framework, or else to introduce new forms of matter which violate the key energy conditions, i.e., the null energy condition (NEC) and the WEC. For a successful bounce, it can be shown that within the context of SCM the NEC and thus the WEC, are violated for a period of time around the bouncing point. In the context of matter bounce scenarios, many studies have been performed using quintom matter [34], Lee-Wick matter [35], ghost condensate field [36], Galileon fields [37] and phantom field [38]. Cosmological bouncing models have also been constructed via various approaches to modified gravity such as f (R) gravity [39,40,41,42,43], teleparallel f (T) grav-ity [44], brane world models [45], Einstein-Cartan theory [46], Horava-Lifshitz gravity [47], nonlocal gravity [48] and others [49]. There are also other cosmological models such as Ekpyrotic model [50] and string cosmology [51] which are alternatives to both inflation and matter bounce scenarios.
In the present work we study the existence of the bouncing solutions in the context of f (R, T) gravity. These type of theories have been firstly introduced in [52] and later, their different aspects have been carefully studied and analyzed in [53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71]. In the current work we use an effective approach that we have introduced previously in [72]. In this method one defines an " effective fluid" endowed with the effective energy density ρ (eff) and effective pressure p (eff) allowing thus, to reformulate the f (R, T) gravity field equations. In a class of the minimally coupled form, i.e., f (R, T) = R + h(T), one usually presumes h(T) functions and solves the resulted field equations. Instead, using the effective fluid description we obtain the h(T) function which corresponds to an EoS defined as p (eff) = Y(ρ (eff) ). We therefore observe that the effective fluid picture may at least be imagined as a mathematical translation of gravitational interactions between the actual matter fields and the curvature to an overall behavior attributed to a mysterious fluid with p (eff) (T) and ρ (eff) (T). In the current article we discuss three different classes of models specified by three different effective EoSs or equivalently three different h(T) functions. We shall see that in these cases we obtain where, β i , λ i , γ,ρ (eff) andp (eff) are some constants. Eliminating the trace between the effective density and pressure leaves us with an EoS as follows where the subscript "eff " has been dropped. In view of relation (5), we may conclude that in f (R, T) = R + h(T) gravity, interactions of a perfect fluid with the spacetime curvature can be mapped effectively onto the behavior of an exotic fluid obeying equation of state (5). It is quiet interesting that a reduced form of (5) has been introduced in [73] and further studied in [74]. Such a complicated EoS has been introduced to study the cosmological implication of a model with a mixture of two different fluids, i.e., effective quintessence and effective phantom. Additionally, one can find other related works investigating some "exotic" fluids which follow various subfamily of (5). These theories may be called "modified equation of state" (MEoS) models. This branch of research presumes a mysterious fluid(s) specified by an unusual EoS with the hope of dealing with some unanswered questions in the cosmological realm. For example some relevant works in the literature can be addressed as follows; in [75] the author has employed an EoS of the form p = −ρ + γρ λ in order to obtain power-law and exponential inflationary solutions. The case with λ = 1/2 has been analyzed in [77,78,74] to focus on the future expansion of the Universe.
Emergent Universe models have been studied in [76] by taking into account an exotic component with p = Aρ − Bρ 1/2 and in [79] with A = −1. Different cosmological aspects of DE with more simple form of EoS, i.e, p DE = α(ρ DE − ρ 0 ) have been investigated in [80] and the study of cosmological bouncing solutions can be found in [18]. Therefore, recasting the f (R, T) field equations into the "effective picture" may provide a bridge to the cosmological models supported by MEoS. Via this connection the problem of an exotic fluid turns into the problem of a usual fluid with exotic gravitational interactions. However, contrary to the former, in the latter case we start with a predetermined Lagrangian, i.e., f (R, T) gravitational Lagrangian. The importance of the effective picture becomes more clear when one considers the energy conditions in f (R, T) gravity. As discussed in [81], in f (R, T) gravity the energy conditions would be obtained for effective pressure and effective energy density. Therefore, it is reasonable to define a fluid as a source with effective pressure and energy density. As we shall see, the bouncing solutions in f (R, T) gravity (using only one perfect fluid in a flat FLRW background) in the framework of our effective fluid approach, exhibit nonsingular properties such that in a finite value of the bounce time t b , non of the cosmological quantities would diverge. More exactly, as t → t b we observe that the scale factor decreases to a minimum non-vanishing value, i.e., a → a b , . Therefore, non of the future singularities would appear. Also, in all cases we have W t→t b → −∞, where W being the effective EoS parameter and subscript "b" stands for the value of quantities at the time at which the bounce occurs. We then observe that if we want to describe nonsingular bouncing solutions in f (R, T) = R + h(T) gravity using a minimally coupled scalar field, a phantom field should be employed. These solutions show a violation of the NEC in addition to the strong energy condition SEC. Such a behavior is predicted in GR for a perfect fluid in FLRW metric with k = −1, 0 [82].
The current research is planned as follows. In Sec. 2 we briefly present the effective fluid picture. Sec. 3 is devoted to the bouncing solutions with asymptotic de Sitter behavior before and after the bounce. We first analyze models with constant effective pressure in subsection 3.1, they are called models of type A. We then proceed to investigate the corresponding bouncing solutions, the energy conditions, the scalar field representation and finally the stability of these type of solutions. In subsection 3.2 models which correspond to two different EoSs assuming p (eff) = Y(ρ (eff) ) will be discussed. These models are named as B, C, D and E models. An example of the matter bounce solution is considered in Sec. 4 which is labeled as model E. The connection of A-E models with MEoS theories will be presented through the effective picture. Section 5, is devoted to study of scalar-type cosmological perturbations. In Sec. 6 we give a brief review of singular models in the context of f (R, T) gravity and obtain a class of solutions exhibiting singular behavior. Finally, in Sec. 7 we summarize our conclusions.
2 Reformulation of f (R, T) field equations in terms of a conserved effective fluid In the present section we review the field equations of f (R, T) gravity theories and rewrite them in terms of a conserved effective fluid. This reformulation would allow us to better understand the properties of bouncing solutions as well as classifying them. In f (R, T) gravity, one usually chooses a mathematically suitable f (R, T) Lagrangian, then obtains the corresponding field equations, and finally tries to solve them. In [72], we introduced a novel point of view to dealing with the field equations of f (R, T) gravity. This approach is based on reconsidering the field equations in terms of a conserved effective fluid. One of the benefits of using this method is to obtain a form of f (R, T) function for a physically justified condition on the effective fluid. Therefore, instead of mathematical arbitrariness in selecting different f (R, T) functions, we have physically meaningful Lagrangian forms. Hence, in this paper we make use of those f (R, T) functions that we obtained in our previous work [72]. The action integral in f (R, T) gravity theories is given by [52] where R, T ≡ g µν T µν , L (m) denote the Ricci scalar, the trace of energy momentum tensor (EMT) and the Lagrangian of matter, respectively. The determinant of metric is shown by g, κ 2 ≡ 8πG is the gravitational coupling constant and we have set c = 1. We assume that the matter Lagrangian L (m) depends only on the metric components therefore the following form for the EMT can be defined By varying action (6) with respect to the metric components g αβ we obtain the following field equation [52] F where we have defined and The differential equations governing the dynamical evolution of the Universe can be obtained from equation (8) by choosing suitable line element. For cosmological applications one can use a spatially flat FLRW metric given as where a(t) is the scale factor of the universe and dΩ 2 is the standard line element on a unit two sphere. Applying metric (11) to field equation (8) together with using the definition for Θ µν for a perfect fluid leads to as the modified Friedmann equation, and as the modified Raychaudhuri equation, where H indicates the Hubble parameter. Note that, we have used L (m) = p for a perfect fluid in the expression (9). Applying the Bianchi identity to the field equation (8) gives following covariant equation where we have dropped the argument of F(R, T) for abbreviation. Equation (14) tells us that the conservation of EMT is not generally respected in f (R, T) gravity theories. There is Only a narrow class of solutions for which the conservation of energy is still preserved [64]. In this work we consider a specific class of models for which the Lagrangian is written as a minimal coupling between the Ricci curvature scalar and a function of trace of EMT, i.e., Also, we study bouncing solutions when the Universe contains a single perfect fluid with a barotropic EoS given as p = wρ. In the next section, we seek for non-singular cosmological solutions in f (R, T) gravity with a baratropic perfect fluid. As we shall see, a useful approach for choosing the functionality of h(T) is to rewriting equations (12) and (13) in terms of an effective fluid that respects the conservation of EMT. For the Lagrangian, (15) equations (12) and (13) are simplified to and where " " denotes derivative with respect to the argument. Also, from equation (14) we obtain As a matter of fact, one can directly solve equations (16)- (18) to obtain the scale factor or the Hubble parameter, once the h(T) function is determined. Alternatively, one can also define the pressure and energy density profiles of an effective fluid along with imposing the energy conditions in order to obtain the functionality of h(T). We then proceed in this way and rewrite equation (16) as where, With this definition the acceleration of expansion of the Universe can be obtained as followsä where, we have defined the effective pressure as Therefore, the original field equations of f (R, T) gravity can be recast into the usual Friedman form with an effective fluid. This fluid is characterized by an effective energy density and effective pressure which in turn are determined in terms of the EMT trace. Therefore, once a property/relation for energy density and pressure components is established, a first order differential equation for the function h(T) could be reached. Solving the differential equation leads to an h(T) function that conveys the specified property/relation. Hence, in this way we can obtain a minimal f (R, T) model based on the conditions on energy density and pressure profiles of the effective fluid, instead of choosing the functionality of f (R, T) based on ad hoc mathematical terms. It is straightforward to verify that equation (18) is turned to the usual conservation equation in terms of ρ (eff) and p (eff) , i.e.,ρ where the arguments are dropped for the sake of simplicity. Consequently, an effective EoS parameter can be defined for this effective fluid as

Asymptotic de-Sitter bouncing solutions in f (R, T) gravity
In this section we study different bouncing cosmological solutions of f (R, T) gravity. We extract those solutions which correspond to some properties of the effective fluid. Such an approach may help us to understand how these solutions can emerge in f (R, T) gravity. Furthermore, there can be obtained more bouncing solutions, however, to show that f (R, T) gravity theories are capable of describing a non-singular pre-Big Bang era, we are restrict ourselves to study only few examples. We set κ 2 = 1 in the rest of the work.

Type
A models: Solutions which correspond to a constant effective pressure, p (eff) (T) = P Let us begin with bouncing solutions which are obtained by assuming an effective fluid with constant pressure. We show that these type of models lead to a de-Sitter era at late times [72]. From definition (22) for a constant effective pressure we obtain Substituting the above function into the modified conservation equation (18) together with using T = (3w − 1)ρ A for a perfect fluid we get where ρ 0 is an integration constant. Substituting solutions (25) and (26) back into the Friedman equation (19) leaves us with the following differential equation for the scale factor For P = 0, a particular solution to the above equation is found as a A (t) ∝ t 2/3 . Moreover, in case in which the middle term is absent, the solution represents a de-Sitter phase for P < 0. From the Friedman equation (27), we observe that the effective pressure P plays the role of a cosmological constant. Other forms for the Friedman equation similar to equation (27) have been found in the literature. For example in [80] authors have worked on a perfect fluid, dubbed as DE with a linear EoS, in a flat FLRW background. They used where α and ρ 0 being some constants. The DE pressure is then decomposed to a dynamical part which corresponds to p α as well as a constant part which corresponds to p Λ . These assumptions lead to an equation that will reduce to (27), when only p Λ is considered. In the upcoming subsections when we study the model C, we find that in the context of f (R, T) gravity one can obtain the same results as those given in [80]. Equation (27) admits three different general solutions which we will consider only one of them. A general solution of equation (27) is obtained as Note that once we set the integration constant in equation (27) such that a A (0) = 0, we get B = 0 and furthermore, if we consider the case w = 0, we obtain A = Pρ 1/3 0 . Thus the solution (28) reduces to the familiar form a(t) = (−ρ 0 /P) 1/3 sinh 2/3 ( −3P/2t). However, the integration constant in solution (28) is fixed such that we have a(t = 0) = a 0 . The most important feature of solution (28) is the appearance of cosine hyperbolic function which allows us to have a non-singular behavior for P < 0. This solution describes the de-Sitter expansion of the Universe which is initially dominated by dust [83]. In Figure 1 the thick black curve shows the scale factor solution (28). This solution behaves exponentially in the far past and far future from the bounce and reaches a nonzero minimum value at the bounce time, given by The Hubble parameter is obtained as Figure 1 shows that the Hubble parameter tends to constant values before and after the bounce and vanishes at the bounce. On the other hand, the Hubble parameter and its time derivatives diverge at an imaginary time t s = iπ/ √ −3P−t b , where t b is the (real) time at which the bounce occurs. Therefore, we observe that the Hubble parameter behaves in a well-defined way (without any future singularity) so that it respects the bounce conditions. We can also obtain the effective energy density and the effective EoS parameter. To this aim, we sub-stitute (25) into equations (20) and (24), respectively. Thus we obtain As can be seen from Figure 1, the effective density reduces from a constant The black middle sized dashed line shows the Hubble parameter, the solid blue curve indicates the effective energy density, the black dot-dashed curve shows the pressure, The long-dashed blue curve shows the effective EoS and the black dotted curve shows the effective energy density. We see that W A < −1 always, signaling that the scalar representation of model A should be defined by a phantom field. We have set w = 0.6, ρ 0 = 1, P = −2.8 and a 0 = 1.12. The horizontal axis represents the cosmic time.
value and tends to zero near the bounce. From equation (19) we see that the vanishing of Hubble parameter at the bounce demands that the effective density becomes zero. Also for the same reason, the effective EoS diverges at the bounce. Such behaviors are common for all bouncing models that we shall present in the framework of minimal f (R, T) gravity. the matter energy density itself increases from small values to a maximum value near the bounce. Based on the exchange of energy between gravitational field and matter constituents (that the mechanism of which is explained in [55]) one may explain the process of bouncing behavior; the interaction of the real fluid with curvature leads to some transformations of energy from gravitational field to matter before the bounce where the spacetime curvature is dominant in comparison to matter energy density. Such a transmutation, that the start of which is triggered at the time far past the bounce, gives rise to an increase in the energy density as the bounce event is approached. At the bounce time the energy density of matter grows to a maximum value after which the process of transmutation is reversed until the density falls back to zero (the post-bounce regime). Note that the effective energy density remains constant in the de-Sitter era. However, some physics is needed in order to explain the process of matter production from curvature component which disturbs the stability of de-Sitter era to enter the bounce event.
Informations from Figure 1 can help us to discuss the energy conditions. In GR the well-known energy conditions are the NEC, WEC, SEC and the dominant energy condition (DEC). In a modified gravity theory with defined effective energy density and pressure, these conditions can be written as [81] WEC ⇔ ρ (eff) ≥ 0, and For the model A, we obtain the following results It is obvious that the fulfillment of NEC and thus WEC requires ( (31) we see that the effective energy density tends to −P at late times and vanishes at the time of bounce. Therefore, we obtain 3P ≤ ρ (eff) + 3p (eff) ≤ 2P, (remember that P < 0). Thus, since only negative values are valid for the effective pressure P, the SEC is always violated. However, the validity of NEC depends upon the value of w. We plot the diagrams for w = 0.6 in Figure 1. This figure shows that the NEC and SEC are violated in this case. Our studies show that the bouncing behavior is achieved from solution (28) for w > 1/3. Note that, as we have mentioned before, solution (28) is only one of the three possible solutions of equation (27). Investigating other solutions may validate the cases with w < 1/3 from energy conditions point of view.
One may be tempted to reinterpret the source of matter as that of a scalar field. Such a representation is also used in similar works [75,79,18,84]. In the case of constant effective pressure, the particular solution (28) (which is valid for 1/3 < w < 1) corresponds to W A < −1 which can be realized from relation (32), see also the long-dashed blue curve in Figure 1. Therefore, if we want to translate mutual interactions of perfect fluid and curvature as the behavior of an effective scalar field, we should employ a phantom scalar field. Therefore, we can define where, the subscript Ph denotes the phantom field and ";" indicates covariant differentiation. We then geṫ A straightforward calculation reveals that where we have defined We can also obtain φ (eff)A in terms of the scale factor by solving solution (28) for time t, which gives We have plotted φ (eff)A and V (eff)A in Figure 2 for the same parameters of Figure 1. Another important issue that needs to be treated is to investigate the stability properties of solution of equation (27). Substituting solution (31) in (19) and also (26) in (17) and taking H and ρ as dynamical variables, we arrive at the following dynamical systeṁ where β = (1 − w 2 )/2(3w − 1) and we have dropped subscripted A . Note that the validity of solution (28) requires that β > 0. The System (46)-(48) has two critical points far from the bounce; P (±) = (± −P/3, 0) which corresponds to the eigenvalues λ (±) = (∓ −P/3, 0). We see that the stability properties of solutions are independent of w. Due to the appearance of the zero eigenvalues, one may not decide about the stability properties of these fixed points, however, by inspecting equations (46) and (47) it is possible to figure out the nature of the fixed points. For fixed point P (−) equation (47) becomeṡ Therefore, in the vicinity of P (−) within the phase space, by indicating the values of ρ and H on vertical and horizontal axises, respectively, we haveḢ > 0 &ρ > 0 for all points with ρ > 0 andḢ < 0 &ρ < 0 for points with ρ < 0. These show that when ρ → ρ + δρ for t → t + δt, the solution at P (−) will not stay stationary and hence it is a repulsive fixed point. On the other hand, for P (+) we havė Therefore, in this case we haveḢ > 0 &ρ < 0 for all points with ρ > 0 anḋ H < 0 &ρ > 0 for points with ρ < 0 in the vicinity of P (+) . Hence it is a stable fixed point. The bounce corresponds to the point P (b) = (0, −P/β) for which we haveḢ = −P/2 > 0. Therefore, at this point the tangent vectors on the phase space trajectories are directed toward the right side. We plot a typical trajectory in the phase space in Figure 3.

Solutions which correspond to a general effective EoS,
These class of models can be constructed by imposing a particular condition on the effective profiles. This approach can be viewed as a sort of classification Generally, one can obtain a class of h(T) functions for a determined property which is specified by an effective EoS. In the following sections we consider two subclasses based on conditions on the effective densities. We find that each class of h(T) solutions that exhibit bouncing behavior correspond to an effective EoS which is already introduced or obtained for an exotic fluid in the literature [18,79,80,84,85].

Type B Models: Solutions which follow the relation dρ
Applying this condition together with using the definitions for the effective energy density and pressure, we arrive at a differential equation for h(T) which can be solved as follows [64] h where, Γ B and Λ B are integration constants and n is an arbitrary constant. Substituting the relation dρ eff /dT = [n/(1 + w)T](ρ eff + p eff ) into equation (18) gives Next we proceed to find a non-singular bouncing solution by solving the modified Friedmann equation (16) for solutions (51) and (52). We first try to obtain solutions of the form a B (t) = R(cosh[(t − t 0 )/R] − S) where R and S are constants. This type of solution has been discussed in [18] under assumption of an MEoS and has the following form of the Friedmann equation Applying (51) and (52) in equation (16) gives the Friedmann equation for arbitrary constants w and n. We can check that, there are only two cases which correspond to equation (53) and thus to the scale factor a B as the solution; when w = −1/5, n = 12/5 and w = −1/5, n = 6/5. However, the latter leads to similar physics to the former. The physical quantities constructed out of the bouncing solution for w = −1/5, n = 12/5 are given as follows The behavior of the above quantities are depicted in Figure 4. In order that the ansatz a B (t) satisfies equation (53) we must have Solution (54) shows that the Universe shrinks from an infinite size to a minimum radius which is equal to R(1 − S). The size of Universe at the time of bounce is controlled by the constant Λ B as well as the coupling constant α. As can be seen, in equation (53) and the subsequent solutions, these two constants appear as a multiplied form. This means that if either of these constants become zero, the bouncing solution would disappear. Also, expression (15) together with solution (51) indicate that the bouncing solution disappears in model B when α = 0.
The expressions for effective energy density and pressure, i.e. (56) and (57), can be rewritten as a sum of three densities and pressures given by From these expressions, we could suppose that our effective fluid consists of a combination of three perfect fluids with EoSs w i = ρ i /p i . In this view, there is a DE component which corresponds to a cosmological constant with w 1 = −1, a quintessence with w 2 = −2/3 and a fluid which drives an expanding Universe with zero acceleration, with w 2 = −1/3. Such a description has been presented in [18]. In the framework of f (R, T) gravity this decomposition may be translated as follows; the effects of unusual interactions of matter (here a perfect fluid with w = −1/5) with curvature, could produce the same behavior as the case of GR with three different fluids. Eliminating the scale factor from solutions (56) and (57) leads to an effective EoS of the form In the cosmological applications, a general type of EoS p = βρ+γf (ρ) is ascribed to some exotic or dark fluid which can determine the evolution of the Universe. Different choices for function f (ρ) are studied in the literature. Such an equation has been already discussed in [85] where the authors considered cosmological consequences of an EoS of the form p = −ρ − ρ q + 1. Substituting expressions (56) and (57) in the NEC and SEC conditions (34) and (35) leads to the following conditions However, expressions (63) and (64) are never satisfied; far from the bounce where a → ∞, the second term of (63) which always has negative sign dominates and in the limit a → a b = R(1−S) the expression for NEC becomes −2/R 2 (1−S) which is also negative because S < 0. The same line of reasons can be used to prove the violation of SEC. Thus, in model B the NEC and SEC are never satisfied. Note that, the NEC is violated only near the bounce, because expression (63) tends to zero as the scale factor gets large values. By considering the above discussion for energy conditions, we find that again, a phantom scalar field can be used for modeling the behavior of bouncing solu-tion. In the case of model B we havė Our studies show that the behavior of the above solutions are similar to those of model A which in Figure 2 a typical example is demonstrated. By choosing H B and ρ B as dynamical variables one can rewrite the Friedmann equation as followsḢ where we have dropped the subscript B. This system admits two critical points P (±) B = (± 1/R, 0) far from the bounce. In these situations we have ρ → 0, hence, the dynamics of equation (67) is determined by the second term. Also, near the point of the bounce as specified by H b = 0, ρ b = 3ρ 2 0 /10R(R + 1), we always haveḢ > 0. Therefore, the stability properties of the system is similar to the model A.
The Friedmann equation (16) for general function (51) takes the form Seeking for a general solution demands that one substitutes the solution (52) into (70) (using the fact that T = (3w−1)ρ) and solves for the resulting differential equation to find the scale factor. However, the resulting equation cannot be solved analytically for arbitrary values of w and n. Nevertheless, for particular values of these parameters a non-singular solution can be obtained as given in expressions (54)- (58). But, we are still able to find more general solutions.
As the third type of solutions named as C, we work on a bouncing solution for which the governing differential equation is given by where we have set w = 1. Choosing a relation between Γ C and ρ 0 , a 0 , n and Λ C as the following where Q is an arbitrary constant and for the Hubble parameter we have By substituting (52) for w = 1 within the definitions (20), (22) and (24) along with using relation (72) we get the effective quantities as follows 1 We have plotted the cosmological parameters (73)-(77) in Figure 5. Unlike 1 We note that as we are concerned with the solutions that respect the conservation equation (18), expression (52) is used for these class of solutions.  (75) and (76) gives the effective EoS which can be viewed as the characteristic equation for the model C. We therefore get Some of the cosmological properties of model C have been investigated in [80].
The authors have considered a model of DE for which an EoS of the form p DE = γ(ρ DE − θ) 2 is assumed for a perfect fluid. We therefore observe that if we apply ρ (eff) → 3ρ DE within equation (19), we obtain the same Friedman equation as the one given in [80]. Also, by redefining the parameters as n → 1/γ and α(1 + n)Λ C /2 → θ in (73) we will obtain the corresponding solution for the scale factor. These considerations show that the problem of dark fluid with an unusual EoS (which may not clearly correspond to a definite Lagrangian) can be explained in the framework of f (R, T) gravity.

Type D models: solutions which are consistent with the relation dρ (eff) /dT = m
Applying this condition on the definition of effective energy density, i.e., definition (20), leads to a second order differential equation for h(T) function; the solution then reads where, m is an arbitrary constant and Γ D , Λ D are integration constants. Substituting (79) into the conservation equation (18), we obtain a first order differential equation for the matter energy density in terms of the scale factor. However, since the mentioned equation cannot be solved for an exact general solution for arbitrary values of w and m, we proceed with particular cases. Note that, further investigations may give other exact solutions or even numerical simulations can be utilized to study other solutions. At the present, we work on a particular case w = 1. The conservation equation (18) then yields where we have again used T = (3w − 1)ρ. In the limit a → ∞ one obtains Γ 2 D α 2 /8m 2 from solution (80). As can be seen, in the model D the matter energy density evolves from a non-vanishing initial value far from the bounce event. This property cannot be seen in the previous models. To obtain the solution for the scale factor, we substitute (79) in the Friedmann equation (19) for w = 1, which gives where From solution (81), we can obtain the time at which the bounce occurs as well as the radius of the Universe at the moment of bounce, as and a (b) Differentiating solution (81) with respect to time gives the Hubble parameter as follows and from definitions of the effective quantities we find By inspection of expressions (87) and (88), we observe that model D corresponds to the following EoS In view of what we discussed in the paragraph right after equation (62), model D corresponds to a generalized EoS with β = 1. This may be interesting since the most bouncing models have been obtained for β = −1. Note that, the behavior of obtained cosmological quantities for two values of ζ as given in (82), are similar to those of model C. By choosing the model parameters so that the term including sinh in (81) disappears, we arrive at a different model. In this case, the behavior of cosmological quantities is the same as the case for which Υ = 0. However, the evolution of matter energy density and the effective pressure are different. We typically plot these quantities for both situations in Figure 6. The thick curves belong to the solution (81) and the thin ones show the solution with Υ = 0. The energy condition considerations show that models of type D lead to the violation of NEC near the bounce event. Also the phase space and the scalar field representation of this model are the same as model A.

Matter bounce solutions in f (R, T) gravity
In this section, we deal with the well-known matter bounce scenario which can be established through the models that obey (51) and (52). We specify these class of solutions as type E models. A branch of the matter bounce scenarios have been discussed with the characteristic scale factor where Q, Z and M are positive constants. Note that Q = uρ max , with u = 2/3, 3/4, 4/3, Z = 1 and M = 1/3 [44,41,42] and M = 1/4 [86] has been used in the literature. The scale factor (91) gives the following expression for the Hubble parameter Note that in model E the Hubble parameter and its time derivatives never diverge since all of them are proportional to negative powers of Qt 2 + Z. By eliminating the time parameter between (91) and (92), we get the Hubble parameter in terms of the scale factor and thus the Friedmann equation can be obtained as From another side, substituting (51) into (16) together with using (52), the Friedmann equation in terms of the scale factor reads 3H 2 (a) = 2αΓ E (w + 1)nρ where we have used subscript E for integration constant Γ E and we have set Λ E = 0. Comparing equations (94) and (93) we find two different type of solutions which only one of them can be accepted. A consistent solution is valid for Eliminating ρ 0 between Z and Q leads to Valid solutions for which {Q, Z, M} > 0 holds, are given by As can be seen, in the context of f (R, T) gravity, there can be found a matter bounce solution for every values of w (except for w = 1/3, 1). The value w = −1 can be accessed for large values of M (see relations (95)). From (52) and (95) we see that for models of type E we have ρ E = ρ 0 a −1/M E . The effective quantities are then obtained as Model E corresponds to an effective EoS as Estimating effective pressure (101) in the limiting times t → 0 and t → ±∞, indicates that the effective EoS follows the expression p − (eff)E far from the bounce and obeys p + (eff)E near the bounce. In t → 0 pressure (101) gives p (eff)E = −4MQ/Z which is consistent with p − (eff)E and in t → ±∞ we have p (eff)E = 0 which can be explained only with p + (eff)E . Figure 7 shows the behavior of different quantities. As is seen in the left panel, the scale factor decreases till reaching a minimum non-zero value at t = t b where the Hubble parameter vanishes. The Universe experiences four phases during its evolution from pre-bounce to post bounce. Before the bounce occurs, the Universe has been within an accelerated contracting regime till the first inflection point (t 1inf < t b ) is reached at which the accelerations vanishes. At this point the Hubble parameter maximizes in negative direction and correspondingly the effective energy density gains a peak value. The Universe then enters a decelerating contracting regime so that its velocity decreases in negative direction. The collapse of Universe halts at the bounce time after which the Universe goes into an accelerating expanding phase where bothä E > 0 and H E > 0. Once the second inflection point (t 2inf > t b ) is reached, the Universe enters a decelerating expanding regime so that the speed of expansion decreases at later times.
Let us now check the behavior of NEC (eff) and SEC (eff) which for model E take the following form as It is obvious that in the expression (104), the sign of Qt 2 − 1 determines the validity of NEC. We therefore find that near the bounce, NEC is violated within the range −1/ √ Q < t < 1/ √ Q, and out of this range it is preserved. The SEC is violated within a larger time interval, i.e., in −1/ Q(1 − 2M) < t < 1/ Q(1 − 2M) for 0 < M < 1/2 and is always violated for M > 1/2. In figure 8 we have plotted the expression of NEC and SEC for two different values of M, i.e., for M = 0.3 and M = 0.7. We see again that, in the background of f (R, T) gravity a bouncing behavior corresponds to violating the NEC. We then may conclude that scalar representation type E models can be constructed by a phantom field within the time interval where NEC is violated (notice that from solution (102) we see that for the same intervals in which the NEC is violated, we have W E < −1). In this cases, the kinetic energy of scalar field is negative. Thus, for a phantom scalar field in model E we obtain In Figure 9 we have presented typical diagrams for the effective phantom field and its corresponding potential. To see the stability properties of model E, we  The effective phantom scalar field and the corresponding potential for the model E. As before, the horizontal axis shows the time values. In the model E, the NEC is violated only near the bounce (within the period −1/ √ Q < t < 1/ √ Q) in contrast to the other discussed models. A phantom scalar field representation can be used to equivalently describe the model which can be valid when the NEC is violated.
rewrite the field equations as follows The analysis of these equations is similar to those of the previous models. Firstly, note that only for 1/6 < M < 1/2, equation of (110) leads to the standard Friedmann equation when the correction term is absent. The system (108) and (109) have two fixed points with coordinates P (±) E = (±0, 0) which correspond to the limit ρ → 0. The bounce event occurs at ρ E = ρ 0 at which we haveḢ E = −3(2M − 1)/6M − 1 > 0. The stability of the system can be analyzed similar to the way used in the previous sections. A simple study shows the stability properties are analogous to the other models; the evolution of the Universe begins from an unstable state, then passing through an unstable bounce phase and finally reaches a stable state.

Stability of the bouncing models
In this section we verify the wholesomeness of the bouncing models which has been introduced in the previous sections through considering the possibility of occurring serious instabilities. We therefore examine the evolution of scalartype perturbations in the discussed models within the metric formalism. Since f (R, T) gravity introduces unusual coupling of matter to curvature part of its action, the evolution of matter density perturbations (specially, as the effect of bounce event on the evolution of matter perturbations is not obvious) can be problematic. In order to study such type of perturbations we consider the matter density perturbations in f (R, T) = R + ακ 2 h(T) models for a fat FLRW metric in the longitudinal gauge where the metric scalar perturbations Φ and Ψ are functions of four coordinates (t, x, y, z), generally. In the current work we shall obtain necessary equations for models including a barotropic perfect fluid with equation of state p = wρ and a general h(T) function. In this respect, the authors of [53] have already considered the matter perturbations in a narrow class of f (R, T) models 3 for a pressure-less perfect fluid. The perturbations of EMT in the longitudinal gauge are given by [87] δT t t = −δρ m , where, v is a covariant velocity perturbation [88]. Using the background equations (16) and (17) we obtain the following equations for the scalar perturbations in Fourier space 2 k 2 a 2 Ψ + 6H(HΦ +Ψ) = δΣ t t + .,E and "eff" are dropped for abbreviation.
as the ADM energy constraint (G t t component of the field equation, if we rewrite (8) as G ν µ = Σ ν µ ), as the ADM momentum constraint (G t i component). Moreover, we have as the ADM propagation equation (G j i − 1/3δ j i G l l component), as the perturbed version of Raychaudhuri equation as the trace equation (G µ µ = Σ µ µ ), as the time component of perturbed EMT conservation and finallẏ v + 3Hσv − Φ + λδ = 0, as the spatial component of perturbed EMT conservation 4 . Note that equations (116)-(122) are the most general equations describing scalar perturbations in minimal f (R, T) gravity for condition F = 1 when a barotropic perfect fluid is included. These equations are not independent, so that it is possible to obtain one equation from another one; for example (120) follows from multiplying (116) by 2 then adding it to (119) and using (127). In the above equations and relations we have used the following definitions for the source terms, which appear in the right hand sides of field equation (8), its trace and in equation (14) when is written as ∇ β T β α = Σ α , respectively The gauge-invariant density contrast in the longitudinal gauge is defined as and the perturbed Ricci scalar for the FRLW metric can be obtained as Also, using definitions in (123)-(125), the coefficients in (121) and (122) can be obtained as follows where and and prime denotes derivative with respect to the trace. Note that wherever is needed we have usedρ which can be obtained from the EMT conservation equation (14). Now, using equations (121) and (122) we arrive at the evolutionary equation for the matter perturbation, as follows, From equations (120) and (127) we obtain a dynamical equation for the perturbed potential Φ, as where we have used solution (118) and defined the following coefficients Hence, we have two differential equations (133) and (134) along with relation (116) to be solved for δ and Ψ. Obviously, the coefficients of these equations are some complicated functions of model properties, H, a, ρ, T, F, F , F , and thus it may not be possible to obtain exact solutions. However, one can resort to numerical methods or even in the case of stiff equations, obtain approximated solutions. We have plotted the evolution of δ and Φ in Figure 10 for models A, B and C. Numerical simulations show that the behavior of perturbations are typically similar in A, B and C models for the same initial values. We have sketched two set of plots in Figure 10. These models show a zero value for the fluctuations when the initial values δ(0) = 0, δ (0) = 1, Φ(0) = 0 and Φ (0) = 1 are assumed. In this case, the fluctuations tend to zero (left panel) and constant values (right panel) in the regime of large times, see the black curves in Figure 10. As another case, the fluctuations increase from zero to a maximum finite value in the period of bounce if the initial data are set as δ(0) = 1, δ (0) = 0, Φ(0) = 1 and Φ (0) = 0, see the gray curves. As can be seen from the evolution of scalar perturbations across the bounce, though small temporary fluctuations, no instability occurs during the bounce, nor does it happen in the limit of large times for these models. Unfortunately, the system of differential equations (133) and (134) become stiff for models D and E so as it is not possible to plot reasonable diagrams for δ and Φ. In this case, we proceed to obtain approximate solutions. For models E and D (in case in which Υ = 0), equations (133) and (134) at times near the bounce will take the following forms for which the solutions up to second order can be found as where the superscript "0" denotes the values of quantities in the limit t → 0 and the subscript i shows the initial values required for integrations. As we see, the solutions are stable in the period of bounce. Note that other coefficients except those which are shown in equations (138) and (139) vanish in the limit t → 0. In the limit of large times, we have numerically plotted the evolution of the matter contrast δ and the potential Φ in Figure 11. As can be seen, from (140) and (141), it is obvious that there happens no instability at the period of bounce in models D and E, however, far away from the bounce point, both δ and Φ increase dramatically in model D, see Figure 11. For model E, in the limit of large times, t → ∞, we geẗ with the solutions Therefore, depending on the initial values, the perturbations in model E can grow before and after the bounce. Therefore we conclude that, though the scalar-type perturbations in models D and E behave regularly at the bounce point, these solutions are unstable asymptotically.

Singular solutions
In this section we seek for possible solutions that exhibit singular behavior, specially the big-bang singularity. Some singular solutions have been already studied in the literature which we suffice to give a short discussion for them. Presuming the EMT conservation (i.e., supposing T = T 0 a −3(1+w) ), equation (18) can be solved to give the following solution where Q 1 and Q 2 are some constants. In [64] the authors have analytically shown that for the case of pressure-less perfect fluid, i.e., for w = 0, one finds which shows a singular behavior. This singular solution is the only one which respects the EMT conservation. Relaxing such a constraining condition, one can obtain other kind of solutions. Note that function (146) is the only form that respects the EMT conservation, that is, to obtain another solution one should assume some suitable form for h(T). For example, in [64] for f (R, T) = R+ακ 2 T the authors have obtained ρ(a) = ρ 0 a 6(α−1) a II SIN (t) = 3 2 2−3α 6(1−α) for a pressure-less matter. In this case singular solution can be found for some valid range of values of the model constant α. Also, the following solution has been obtained for models of form f (R, T) = R + ακ 2 T −1/2 ρ(a) = 2 −2/3 α + 2ρ 0 a III SIN (t) = 27 256 As a new class of (big-bang) singular models, we examine the models which admit the following solution where being a real positive number and t 0 being the time at which the scale factor gets the present value a 0 . Any solution has to satisfy two of the three equations (16)- (18), for a presumed scale factor function. In this case, two unknown functions T(a) and h(T) should be obtained via solving the resulted equations. To proceed further we make use of the following anzats where constants C i 's and µ should be fixed according to the following considerations. Substituting (152) and (153) in (18) and solving for T(a) we get where an integration constant has been set so that T(a = a 0 ) = T 0 . To ensure that three functions (152), (153) and (154) provide an analytic solution, they must satisfy at least one of the equations (16) and (17). Substituting these functions in (16) shows that we have a solution only for the case of stiff fluid, i.e., for w = 1. Therefore, a singular solution can be obtained provided that Thus, for conditions (155) we obtain (157) Therefore, besides the non-singular solutions obtained in the present paper, one can still find a set of singular solutions. In this brief section in addition to addressing some previous results we obtained a new singular model, though a coherent study can be performed in order to deal with possible conditions for which big-bang singularity would occur; however working on this issue is beyond the scope of the present paper and comprehensive studies on this subject will be reported elsewhere. It is worthwhile to mention that some studies have been already made to consider other forms of singular solutions, e.g., [90]. Beside the above results, Bianchi type I cosmological model with magnetized strange quark matter in the framework of f (R, T) gravity have been investigated and it is found that the model begins with big-bang and ends with big rip [91]. Using Lie point symmetry analysis method, the authors of [92] have shown that for a Bianchi type I spacetime, both singular (big-bang) and nonsingular solutions could exist subject to the type of specified symmetry.
Recently, the authors of [93], have considered some cosmological features of f (T ) gravity (where here T denotes torsion scalar) using the dynamical system approach both generally and for some specific forms of f (T ) functions. The core of their studies is taking the advantage of this fact that the torsion scalar can be used interchangeably with the Hubble parameter (i.e., T = −6H 2 ). Thus, the field equations reduce to a single equation (in the case of pressure-less matter) in the form ofḢ = F(H), since the matter density can also be rewritten as a function of the Hubble parameter. Briefly, they have shown that in f (T ) gravity a single equation (which can be interpreted as a simple one dimensional dynamical system) can govern the dynamics of field equations. Benefiting this useful result they investigated phase space portraits of various cosmological evolutions such as, singular and non-singular solutions. Likewise, one may be motivated to utilize such an approach in order to investigate the cosmological solutions of f (R, T) gravity (especially, in the case of present work, i.e., the function given in (15)) through phase portrait diagrams. However, looking at equations (17) and (18) one finds that it is impossible to obtain an equation likė H = F(H) so that it reflects full information of the field equations. In this case for an assumed function h(T), we have a two dimensional dynamical system without any further reduction. Thus, the procedure proposed in [93] would be generally failed in f (R, T) gravity.

Concluding remarks
In the present work we studied classical bouncing behavior of the Universe in the framework of f (R, T) = R + h(T) gravity theories. We assumed a single perfect fluid in a spatially flat, homogeneous and isotropic FLRW background. Having obtained the resulted field equations, we employed the concept of effective fluid (which is firstly introduced in [72]) via defining an effective energy density and pressure and also reformulating the field equations in terms of these fluid components. In this picture, one could recast the field equations of f (R, T) gravity for a real perfect fluid into GR field equations for an effective fluid. It is also shown that in a modified gravity model the energy conditions are usually obtained by using the effective EMT, not the one for real fluids. In f (R, T) gravity, the definitions for effective energy density and pressure have already been used to obtain the energy conditions [81]. The effective fluid has an EoS of the form, p (eff) = Y(ρ (eff) ), which corresponds to an h(T) function. In this method one firstly specifies an effective EoS or a condition on the effective components and then obtains the corresponding h(T) function and other cosmological quantities.
It is also possible to make a link between f (R, T) gravity in effective picture and models which use some exotic or dark component with unusual EoS. These models which have been widely discussed in the literature (to deal with some cosmological issues) are also called theories with "generalized EoS". The mathematical representation of effective components provides a setting within which unusual interactions of a real perfect fluid with gravitational field can be translated as the presence of an exotic fluid which admits the EoS of the form p (eff) = Y(ρ (eff) ). In this paper we have shown that it is possible to recover generalized EoS models which have been previously studied in the literature ( see e.g., [18,79,80,84,85]), in the framework of f (R, T) gravity. Therefore, the problem of exotic fluid in the context of generalized EoS models which are mostly without a determined Lagrangian may be discussed in a Lagrangian based theory of gravity like f (R, T) gravity.
In the current research, we discussed four different bouncing models in f (R, T) gravity. We labeled them as the models A, B, C, D and E and briefly mentioned their main properties in Table 1. Each model can be specified either by an h(T) function or by an effective EoS. Models A-D mimic an asymptotic de Sitter expansion in the far past and future of the bounce. The model A corresponds to a constant effective pressure, p (eff) = P; for the model B we have p (eff)B = −ρ (eff)B /3 + b B ρ (eff)B + d B + e B , the model C is specified by p (eff)C = j C ρ (eff)C + e C , the model D corresponds to p (eff)D = ρ (eff)D + b D ρ (eff)D + d D + e D and finally the model E obeys the EoS, p (eff)E = a E ρ (eff)E + b E ρ (eff)E + d E + e E , where the constants b, d, e and j are written in terms of model parameters. In all models the matter density grows to a maximum value at the bounce which corresponds to a minimum for the scale factor. The effective density varies from zero at the bounce to a positive value in the far past and future of the bounce. The effective pressure varies between negative values; in model A, it is a constant, in the model B it increases at the bounce, in the models C and E it decreases and the model D admits both behaviors. The effective EoS has the property −∞ < W < −1 when the bounce point is approached. The Hubble parameter satisfies H(t) = 0 and dH/dt > 0 at the event of bounce and also all its time derivatives have regular behavior for all models. Therefore, these bouncing solutions do not exhibit future singularities which are classified in the literature of cosmological solutions. We can consider the inherent exoticism hidden behind f (R, T) gravity in another way. As already we mentioned, this issue can be described as an unusual interaction between gravitational field and normal matter or introducing an effective fluid. From the point of view of the energy conditions, in all discussed models the SEC and NEC are violated (note that for a normal fluid NEC is not violated [15]) near the bounce and the effective density gets minimized to zero. Such a result has been previously predicted in GR [82]. As discussed in [15], the exoticness can be understood as a minimization in the effective pressure. In the other words, a minimum in the effective energy density corresponds to a minimum in the scale factor. Such a behavior is permitted provided that W < −1. Note that, for a normal matter, a minimum (maximum) compression leads to a minimum (maximum) energy density. Thus, in f (R, T) gravity an abnormal or effective fluid which leads to an uncommon balance in the density and pressure can be responsible for the bouncing behavior. An interesting feature of the bouncing solution in f (R, T) gravity is that one can construct solutions in which the SEC is respected by the real perfect fluid. Such solutions cannot be found in GR [82]. Also note that the real perfect fluid with w > −1 never violates the NEC. Therefore, we have solutions without the future singularities and all energy conditions can be respected by a real perfect fluid. By this discussion, one may use the definition of an (effective) phantom scalar field if one asks for the matter source to be reinterpreted as that of a scalar matter field. We obtained the equivalent scalar field φ (eff) (t) and its corresponding potential V (eff) (t) in each case. Moreover, we have studied the dynamical system representation of these models. We found that the evolution of the Universe can be displayed by trajectories which initially start from an unstable state, passing through an unstable fixed point (the bounce event) and finally are absorbed by a stable point. The initial and final states are de-Sitter era in models A, B, C and D and the decelerated expanding Universe in model E. Another important issue discussed in this work is related to the study of stability of bouncing solutions through scalar-type cosmological matter perturbations in the bouncing universe. Our numerical analysis of density perturbations for models A, B and C revealed that, though a slight jump (depending on the initial conditions) at the bounce point, the amplitude of matter density perturbation (δ) and perturbed potential (Φ) behave regularly throughout the bounce phase. Therefore, since the time interval during which the fluctuations that occur within density contrast and perturbed field is short, the instabilities do not have enough time to grow to a significant magnitude. However, this case does not happen for the two remaining models.
As the final remarks we should emphasize that our models were obtained by indicating different conditions on the effective density and pressure which led to different h(T) functions. This means that the models A, B, C, D and E are not the only possible models for the bouncing behavior. It is obvious that one can still choose other h(T) functions or consider other assumptions on the effective density and pressure to obtain new bouncing solutions (with even new features). Our aim was to show the existence of varieties of bouncing solutions in f (R, T) gravity and study their properties. Especially, our study was confined to the Lagrangians of type f (R, T) = R + h(T) though other forms of Lagrangians can be investigated. The other issue is that our study was performed in the effective picture. In case such an approach is not taken seriously, one can think of it as only an alternative mathematical method. One can still investigate a nonsingular cosmological scenario without employing the equations which are written in terms of the effective quantities. In this case it is enough to assume a Lagrangian and solve the field equations to inspect for a bouncing solution. However, cosmological solutions for the f (R, T) gravity model presented here are not singularity free and as we observed under certain conditions, a class singular solutions could be obtained.