In this work we have investigated the dynamics of a recent modification to the general theory of relativity, the energy-momentum squared gravity model f (R, T 2 ), where R represents the scalar curvature and T 2 the square of the energy-momentum tensor. By using dynamical system analysis for various types of gravity functions f (R, T 2 ), we have studied the structure of the phase space and the physical implications of the energy-momentum squared coupling. In the first case of functional where f (R, T 2 ) = f0R n (T 2 ) m , with f0 constant, we have shown that the phase space structure has a reduced complexity, with a high sensitivity to the values of the m and n parameters. Depending on the values of the m and n parameters, the model exhibits various cosmological epochs, corresponding to matter eras, solutions associated to an accelerated expansion, or decelerated periods. The second model studied corresponds to the f (R, T 2 ) = αR n + β(T 2 ) m form with α, β constant parameters.In this case it is obtained a richer phase space structure which can recover different cosmological scenarios, associated to matter eras, de-Sitter solutions, and dark energy epochs. Hence, this model represent an interesting cosmological model which can explain the current evolution of the Universe and the emergence of the accelerated expansion as a geometrical consequence.The second approach aims at modifying the geometry of spacetime, i.e. the Einstein's gravity in the general theory of relativity (GR) at large distances, specifically beyond our Solar System to produce accelerating cosmological solutions [5][6][7]. This has given rise to the concept of "modified gravity". There are numerous such theories available in literature, such as Brans-Dicke theory [8], Einstein-Cartan theory [9], Brane gravity theory [10-13], Rosen's Bimetric theory of gravity [14][15][16], Kaluza-Klein theory [17][18][19][20][21], Hořava-Lifshitz theory [22,23], etc. Extensive reviews in modified gravity theories can be found in the Refs. [24,25]. All these have their own share of merits and de-merits. Many of the theories of modified gravity aims at modifying the linear function of scalar curvature, R responsible for the Einstein tensor in the Einstein equations of GR. So it is obvious that the alterations are brought about in such a way so as to generalize the gravitational Lagrangian which takes a special form L GR = R in case of GR. An extensively studied theory in this context is the f (R) gravity where the gravitational Lagrangian L GR = R is replaced by an analytic function of R, i.e., L f (R) = f (R). Via this generalization, we can explore the non-linear effects of the scalar curvature R in the evolution of the universe by choosing a suitable function for f (R). Extensive reviews on this theory is available in the Refs. [26,27].The cosmological viability of f (R) dark energy models have been studied in [28]. In this paper the authors ruled out the f (R) theories where a power of R is dominant at large or small R. Conditions for the cosmological via...