Generalised nonminimally gravity-matter coupled theory

In this paper, a new generalised gravity-matter coupled theory of gravity is presented. This theory is constructed by assuming an action with an arbitrary function $f(T,B,L_m)$ which depends on the scalar torsion $T$, the boundary term $B=\nabla_{\mu}T^{\mu}$ and the matter Lagrangian $L_m$. Since the function depends on $B$ which appears in $R=-T+B$, it is possible to also reproduce curvature-matter coupled models such as $f(R,L_m)$ gravity. Additionally, the full theory also contains some interesting new teleparallel gravity-matter coupled theories of gravities such as $f(T,L_m)$ or $C_1 T+ f(B,L_m)$. The complete dynamical system for flat FLRW cosmology is presented and for some specific cases of the function, the corresponding cosmological model is studied. When it is necessary, the connection of our theory and the dynamical system of other well-known theories is discussed.


I. INTRODUCTION
Nowadays, one of the most important challenges in physics is try to understand the current acceleration of the Universe. In 1998, using observations from Supernovae type Ia, it was shown that the Universe is facing an accelerating expansion, changing the way that we understand how our Universe is evolving [1]. Later, other cosmological observations such as CMB observations [2][3][4][5], baryon acoustic oscillations [6] or galaxy clustering [7] also confirmed this behaviour of the Universe. The responsible of this late-time acceleration of the Universe is still not well understood and for that reason it was labelled as the dark energy problem. In general, there are two different approaches which try to deal with this issue. First, one can assume that General Relativity (GR) is always valid at all scales and introduce a new kind of matter which mimics this acceleration. This kind of matter known as "exotic matter" needs to violate the standard energy conditions to describe the evolution of the Universe. Up to now, this kind of matter has not been discovered in the laboratory. One can say that this approach lies on the idea of changing the right hand side of the Einstein field equations. An alternative approach to understand and study the dark energy is to assume that GR is only valid at certain scales and therefore it needs to be modified. In this approach, the left hand side of the Einstein field equations is modified and there is no need to introduce exotic matter. Different kind of modified theories of gravity have been proposed in the literature to understand the dark energy problem (see the reviews [8,9]).
One very interesting and alternative theory of gravity is the teleparallel equivalent of general relativity (TEGR) or "teleparallel gravity". In this theory, the manifold is endorsed with torsion but assumes a zero curvature. The connection which satisfies this kind of geometry is the so-called "Weitzenböck" connection, which was first introduced in 1922 [10]. It was then showed that this theory is equivalent to GR in the field equations but the geometrical interpretation of gravity is different. In TEGR, there is not geodesic equation as in GR. Instead, forces equations describe the movement of particles under the influence of gravity. Additionally, the dynamical variable is the tetrad instead of the metric as in GR. For more details about TEGR, see [11][12][13][14][15][16] and also the book [17]. Similarly as in GR, there are also modified theories starting from the teleparallel approach. The most famous teleparallel modified theory is f (T ) gravity (where T is the scalar torsion) which can describe very well the current acceleration of the Universe and also other cosmological observations (see [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32] and also the review [33]). The TEGR action contains the term T so f (T ) gravity is a straightforward generalisation of it. This theory is analogous to the well-known f (R) gravity, where instead of having the scalar curvature R in the action, a more general theory with an arbitrary function which depends on R is introduced. These two theories are analogous but mathematically they are very different. As we pointed out before, the TEGR field equations are equivalent to the Einstein field equations. However, their generalisations f (R) and f (T ) gravity have different field equations. Further, f (R) gravity is a 4th order theory and f (T ) gravity is a 2nd order theory. This characteristic can be understood using the fact that R = −T + B, where B is a boundary term. Hence, a linear combination of R or T in the action will produce the same field equations since B will not contribute to it. However, when one modifies the action as an arbitrary function f (T ) or f (R), there will be a difference in their field equations due to the fact that now the boundary term B contributes. This was fully studied in [34] where the authors introduced a new theory, the so-called f (T, B) gravity, which can recover either f (T ) gravity or f (T, B) = f (−T + B) = f (R) as special cases. Flat FLRW cosmology of this theory was studied in [35,36].
Other kinds of modified theories of gravity have been considered in the literature. Some interesting ones are theories with non-minimally coupling between matter and gravity. In standard metric approach, some alternatives models have been proposed such as f (R, T ) [37], where T is the trace of the energy-momentum tensor or non-minimally coupled theories between the curvature scalar and the matter Lagrangian f 1 (R) + f 2 (R)L m [38]. Further, another more general theory is the so-called f (R, L m ) where now an arbitrary function of R and L m is considered in the action [37]. Along the lines of those theories, modified teleparallel theories of gravity where couplings between matter and the torsion scalar have been also considered. Some important theories are for example: f (T, T ) gravity [39] and also non-minimally couplings between the torsion scalar and the matter Lagrangian theory f 1 (T ) + f 2 (T )L m [40]. Along this line, in this paper, we present a new modified teleparallel theory of gravity based on an arbitrary function f (T, B, L m ) where L m is the matter Lagrangian. In this theory, we have the possibility of for example recover f (−T + B, L m ) = f (R, L m ) or a new generalisation of [40] in the teleparallel framework with a function f (T, L m ) depending on T and L m . The later new theory is the analogous theory as f (R, L m ) gravity. We will explicitly discuss about how those models are related, with B being the main ingredient which connects both the metric and tetrad approaches.
After formulating the new f (T, B, L m ) theory, the conservation equation is obtained and exactly as in f (R, L m ), the conservation equation in f (T, B, L m ) theory is not always valid. It will be proved that for the flat FLRW case and assuming L m = −2ρ, the conservation equation is conserved exactly as happens in f (R, L m ) or in f 1 (R) + f 2 (R)L m (see [41,42]). The main aim of this paper is to also formulate the dynamical system of this new generalise theory, which is in general a 10-dimensional one. This dynamical system is a generalisation of different models such as the ones studied in [42][43][44]. After formulating the full dynamical system, different special cases are recovered. Some of them have been studied in the past, hence we only mention how our dimensionless variables are related to them and then we show that our dynamical system becomes them for the special case studied. Then, using dynamical system techniques, we will study new cases that can be constructed from our action. Similarly as in f (R, L m ) (see [42]), a power-law and a exponential kind of coupling between L m and T is studied. Additionally, another new kind of couplings between the boundary term B and L m are studied. For this theory, we study different power-law models with f (T, B, L m ) = C 1 T + C 5 B s + (C 4 + C 4 B q )L m . This model depends highly on the power-law parameters s and q. The critical points and their stability are then studied for different models. For the readers interested on dynamical systems in cosmology, see the review [45] and also see [46,47] for further applications to dynamical systems in modified teleparallel models with the boundary term B.
The notation of this paper is the following: the natural units are used so that κ = 1 and the signature of the metric is η ab = (+1, −1, −1, −1). The tetrad and the inverse of the tetrad are labelled as e a µ and E µ a respectively where Latin and Greek indices represent tangent space and space-time coordinates respectively. This paper is organized as follows: Sec. II is devoted to present a very brief review of teleparallel theories of gravity and some interesting modified theories than can be constructed from this approach. In Sec. III is presented the new generalised gravity-matter coupled theory of gravity known as f (T, B, L m ) where T, B and L m are the scalar torsion, the boundary term and the matter Lagrangian respectively. The corresponding field equations of the theory and the flat FLRW cosmological equations are also derived in this section. In Sec. IV is presented the dynamical system of the full model and for some specific theories, the corresponding dynamical analysis of them is performed. Finally, Sec. V concludes the main results of this paper.

II. TELEPARALLEL GRAVITY AND ITS MODIFICATIONS
Let us briefly introduce the teleparallel equivalent of general relativity (TEGR) and some important modifications under this theory. Basically, this theory is based on the idea of having a globally flat manifold (zero curvature) but with a non-trivial geometry for having a non-zero torsion tensor. Hence, the concept of paralellism is globally defined in TEGR. The dynamical variable of this theory is the tetrad which defines orthonormal vectors at each point of the manifold and they are directly related with the metric as follows where η ab is the Minkowski metric. The connection which defines a globally flat curvature with a non-vanishing torsion is the so-called Weitzenböck connection W µ a ν , which defines the torsion tensor as taking its anti-symmetric part, namely Let us clarify here that the above definition is not the most general form of the torsion tensor. The most general definition also contains the spin-connection which needs to be pure gauge in order to fulfil the condition of teleparallelism (zero curvature). In this paper is assumed that the spin-connection is identically zero. The TEGR action is defined with the so-called torsion scalar T as follows where e = det(e a µ ) and S m is the matter action. The torsion scalar is defined as the contraction of the super-potential with the torsion tensor as T = T abc S abc . Here, T µ = T λ λµ is the so-called torsion vector. The definition of T comes directly from the condition of zero-curvature where one arrives that the Ricci scalar is directly linked with it via where B refers to the boundary term which connects the Ricci scalar with the torsion scalar. From (3) and the above relationship, one can directly notice that the TEGR is equivalent to the Einstein-Hilbert action up to a boundary term. Hence, TEGR is an alternative formulation of gravity which reproduces the same field equations as GR. Although, the geometrical interpretation of these theories are different. GR lies in a manifold with a non-zero curvature (in general) with a zero torsion tensor whereas TEGR is the opposite. Moreover, geodesic equations are replaced by forces equations in TEGR (see [17] for more details about this theory). A straightforward generalisation of the action (3) is to replace T by an arbitrary function of f which depends on T , namely The former theory is the most popular modification of TEGR and it was firstly introduced in [19] with the aim to study inflation in cosmology. In some sense, this generalisation is analogous as the famous modification of GR, the socalled f (R) gravity, where instead of having R in the Einstein-Hilbert action, an arbitrary function of R is introduced in the action. The formulation described here for f (T ) gravity where the spin-connection is identically zero is not invariant under Lorentz transformations. This is due to the fact that T itself is not invariant under local Lorentz transformations so f (T ) gravity will also have this property [48,49]. In standard TEGR where T is in the action, this problem is not important since the action only differs by a boundary term with respect to the Einstein-Hilbert action so one can say that this theory is quasi-invariant under local Lorentz transformation. The problem of the loose of the Lorentz invariant produces that two different tetrads could give rise different field equations so it depends on the frame used. For example, the flat FLRW in spherical coordinates give rise to different field equations as in Cartesian coordinates. At the level of the field equations, this problem can be alleviated by choosing "good tetrads" as it was introduced in [50]. In this approach, one needs to rotate the tetrad fields and fix it accordingly depending on the geometry studied. In [51], it was proposed a new approach of teleparallel theories of gravity where a non-zero spin-connection is assumed giving rise to a covariant version of f (T ) gravity. Both approaches should arrive at the same field equations and since almost all the works based on f (T ) gravity used the approach presented above, we will continue using this approach. It is also possible to create other kind of modifications of teleparallel theories of gravity. A very interesting modification theory is given by the following action [34] where now the function also depends on the boundary term B. Under this theory, it is possible to recover either f (−T + B) = f (R) gravity or f (T ) gravity. Moreover, the theory f (T, B) = C 1 T + f 1 (B) can also be obtained from this action. From this theory one can directly see how f (R) and f (T ) are connected by this boundary term. Since R = −T + B, only if a linear combination of R and T is assumed in the action (TEGR or GR), we will have equivalent theories at the level of the field equations. It is known that f (R) gravity is a 4th order theory whereas f (T ) gravity is a 2nd order theory. Hence, f (T, B) gravity is also a 4th order theory. f (T ) and f (R) gravity have different field equation orders since the difference comes from integrating by parts twice the boundary term B.

A. General equations
Inspired by the theories described in [52] in the curvature approach and also from f (T, B) gravity, let us now consider the following gravity model where the function f depends on the scalar curvature T , the boundary term B and the matter Lagrangian L m . The energy-momentum tensor of matter T β a is defined as Now, we will assume that the matter Lagrangian depends only on the components of the tetrad (or metric) and not on its derivatives, giving us Now, by a variation of action (8) with respect to the tetrad, we obtain where we have used Eq. (10) and f T = ∂f /∂T, f B = ∂f /∂B and f L = ∂f /∂L m . Variations with respect to the torsion scalar and the boundary term are given by [34] so that by imposing δS f (T,B,Lm) = 0, we obtain the f (T, B, L m ) field equations given by The above field equations can be also written only in space-time indices by contracting it by e a λ giving us From these field equations, one can directly recover teleparallel gravity by choosing f (T, L m ) = T + L m which gives us the same action as (3). Moreover if we choose f (T, L m ) = T + f 1 (T ) + (1 + λf 2 (T ))L m we recover the non-minimal torsion-matter coupling extension of f (T ) gravity presented in [40]. Note that in our case, we have assumed that the matter Lagrangian does not depend on the derivatives of the tetrads, which according to [40] is equivalent as having Let us now study the conservation equation for this theory. First, we will use that R β is the Einstein tensor. Using this relationship, we can rewrite the field equation (16) as follows where for simplicity we have also introduced the quantity By taking covariant derivative of H λβ and after some simplifications, we find that where we have used the fact that the energy-momentum tensor is symmetric and hence S σρ λ K βσρ X λ = 0. The latter comes from the fact that field equations are symmetric, and hence the energy-momentum tensor is also symmetric. Now, we will find the condition that f needs to satisfy in order to have the standard conservation equation for the energy momentum tensor, i.e., ∇ µ T µν = 0. By taking covariant derivative in (18) and assuming ∇ µ T µν = 0, one gets that the standard conservation equation for the energy-momentum tensor is satisfied if the function f satisfy the following form which matches with the conservation equation presented in [52]. Note that in our case, we have defined the energymomentum tensor in a different way so that there is a minus sign of difference between Eq. (13) presented in [52] and the above equation. Thus, in general, f (T, B, L m ) is not covariantly conserved and depending on the metric, the model and the energy-momentum tensor, this theory may or may not be conserved. Hereafter, we will consider that the matter is described by a perfect fluid whose energy-momentum tensor is given by Here, ρ and p are the energy density and the pressure of the fluid respectively and u µ is the 4-velocity measured by a co-moving observer with the expansion so that it satisfies u µ u µ = 1. For a perfect fluid, if one assumes that in the proper frame where the particle is static, the matter Lagrangian is invariant under arbitrary rescaling of time coordinate [53]. Therefore, from (10), one gets T 00 = ρ = −(1/2)L m which is equivalent as having L m = −2ρ. This is a "natural choice" for a perfect fluid (see [40,41,53] for more details). Hence, from Eq. (21) we can directly conclude that the conservation law will be always satisfied when flat FLRW and a perfect fluid are chosen without depending on the model for the function f (T, B, L m ). This statement was also mentioned in [52], which is a special case of our theory, explicitly when

B. Flat FLRW cosmology
In this section we will briefly find the corresponding modified flat FLRW cosmology of our theory. Consider a spatially flat FLRW cosmology whose metric is represented by where a(t) is the scale factor of the universe. The tetrad corresponding to this space-time in Cartesian coordinates reads e a β = diag(1, a(t), a(t), a(t)) .
For the space-time given by (23), the modified FLRW equations become where H =ȧ/a is the Hubble parameter and dots represent derivation with respect to the cosmic time. Note that the termsḟ  [42]. Note that in the latter paper, the authors used another signature notation η ab = (− + ++), so that one needs to change R → −R to match those equations. As a consequence of the conservation law holds when considering a perfect fluid as a matter content of the universe, we also know that the standard continuity equation is valid in our case. Hence, we have that the fluid satisfieṡ Let us now assume a barotropic equation of state p = wρ, so that we can directly find that the energy density of the fluid behaves as where ρ 0 is an integration constant. It is also useful to note that the scalar torsion and the boundary term in this space-time satisfy the relationship (5), namely IV. DYNAMICAL SYSTEMS A. Dynamical system for the full theory In this section we will explore the dynamical system of different theories of gravity coupled with matter. To do this, we will first study the dynamical system of the general modified FLRW by using the conservation equation given by (27) and also the first modified FLRW equation (25). By replacing the boundary term given by Eq. (29) in (25) and expanding the derivatives of f we get where we have used the conservation equation (27) to replaceL m = −2ρ = 6Hρ(1 + w). Let us now introduce the following dimensionless variables These dimensionless variables were chosen with the aim of having a similar variables as the ones presented in [42]. Further, using these variables will help us to compare both theories in the limit case where Using these variables, the Friedmann constraint given by (30) becomes Moreover, using the dimensionless variables, we can find the following useful relationṡ where we have defined N = ln a as the number of e-folding so that d/dt = Hd/dN . The effective state matter and the deceleration parameter can be written in terms of these dimensionless parameters as follows For acceleration universes, one needs thatq < 0 or equivalently w eff < −1/3. By replacing the identities (34)- (37) and the dimensionless variables defined as (32) into the second modified Friedmann equation (26), we get The conservation equation (27) can be also written in terms of the dimensionless variables, which yields where we have defined The other quantities defined above will be useful in the full dynamical system equations. Using the dimensionless variables, the identities mentioned before and the above definitions, one can find the dynamical system. This procedure is very involved since the dynamical system is a 10 dimensional one. After all of those computations, the system can be summarized with the following equations: Additionally, one can use the Friedmann constraint (33) to reduce the above system to a 9-dimensional one. In the following sections we will explore the dynamical system of different kind of matter coupled theories of gravity which can be obtained from our approach. From our action it is possible to recover a very interesting model which comes from the curvature approach. As we discussed before, this is possible due to the fact that the function also depends on the boundary term B so that it is always possible to reconstruct theories which contains the scalar curvature R. In this sense, one can construct a non-minimally coupled theory between the Lagrangian matter and the scalar curvature, explicitly by taking the function being This kind of models was first proposed in [52] where the authors suggested that in general, this theory has some extra terms in the geodesic equation. The complete analysis of the dynamical system for flat FLRW for this model was studied in [42]. From our full dynamical system, it is possible to recover the same dynamical system equations reported in the former paper. Let us recall here again that in general, the signature of the metric for curvature theories as f (R) gravity is usually taken as the opposite as in teleparallel theories of gravity and hence the paper [42] is written in another signature notation compared to our notation. The important of this notation issue is that the scalar curvature, scalar torsion and also the boundary term will have a minus sign of difference with respect to our case. Hence, in their notation Eqs. (45)-(54) will have a minus of difference in all those quantities. Therefore, to recover the same dynamical systems found in [42] we need to change R → −R (and of course T → −T and B → −B) which makes that the important derivatives appearing in the dimensionless variables become In this case, some dimensionless variables can be reduced. It is possible to connect our dimensionless variables to the dimensionless variables used in the mentioned paper by working with the variables where tildes represent the variables chosen in [42]. By replacing (56) and (32) in our dynamical system we directly find that the corresponding dynamical system becomes a 5-dimensional one explicitly given by The above equations are the same reported in [42] for f (R, L m ) gravity if one changes the variables β BBL , β BLL , α T B and β LL accordingly. In this theory, it is possible to reconstruct different interesting gravity-matter coupled models as for example standard non-minimally curvature-matter coupled models where f (R, L m ) = f 1 (R) + f 2 (R)L m which has been studied in the literature (see [38]). In [43], it was studied the dynamical system for some of those models. To find all the most important details regarding the above dynamical system for the former model and also for other more general models in f (R, L m ) gravity, see [42]. Let us now introduce a new theory of gravity based on an arbitrary function f which depends on T and L m only. As f (T ) was motivated by f (R) gravity, f (T, L m ) gravity is somehow, the teleparallel version of f (R, L m ) discussed in the previous section. Different particular cases of this theory have been studied in the past. Let us first derive the full dynamical system for the f (T, L m ) gravity and then study some particular theories. The Friedmann equation (33) for this model reads In this case, y 1 = y 2 = x 1 = x 2 = x 3 = φ ≡ 0, so one needs to be very careful with the general dynamical system (45)-(54) since some of these equation will be also identically zero. Let us clarify here the way that one needs to proceed to find the correct dynamical system. There are two ways to find out the correct dynamical system for an specific model. Let us discuss how to proceed with the model that we are interested here, i.e., where f = f (T, L m ). The first way to proceed is using the full dynamical system described by (45)- (54). If one directly replaces f = f (T, L m ) in the full dynamical system, there will be some expressions that are indeterminate or directly zero, for example terms like y 1 /y 2 or terms divided by x 2 . Hence, one first needs to replace back all the original definitions of the dimensionless variables and after doing that, one can restrict f = f (T, L m ). By doing that, several equations are directly satisfied. Indeed, one can verify that Eqs. (45), (47), (49), (50) and (52) are identically zero, as expected. Then, for all the remaining equations, one needs to introduce again the dimensionless variables needed (in this case x 4 , α and θ). A second approach is to directly assume f = f (T, L m ) in the Friedmann equations (30)- (26) and then introduce the same dimensionless variables that we defined. By doing that, we arrive of course at the same dynamical system as the first approach. We will implement the first procedure in this work. Eq. (46) gives us a constraint for the variables, namely Let us here clarify again that even though y 1 = y 2 ≡ 0, the quotient y = y 1 /y 2 = B/T = 3 +Ḣ/H 2 is clearly non-zero. If we replace the above equation and also use the Friedmann constraint (63), the remaining three Eqs. (48), (53) and (54) gives us the following set of equations, where for convenience we have introduced the following quantities In the following section we will study some interesting cases that can be constructed from this theory.

Nonminimal torsion-matter coupling f (T, Lm) = f1(T ) + f2(T )Lm
In this section we assume that the function takes the following form, where f 1 (T ) and f 2 (T ) are arbitrary functions of the torsion scalar. This model is an extension of f (T ) gravity, where an additional nonminimal coupling between the torsion and the matter Lagrangian is considered [40]. In this model we have that γ LL =β T LL ≡ 0. In [44], the dynamical system of this model was carefully studied. The authors used other dimensionless variables but one can verify that our dynamical system (65)-(67) give rise to the the same dynamics. In that paper, the authors used a different energy density which is related to our as ρ new = −2ρ. The dimensionless variables used in [44] are given by It is possible to show that our variables are related to them as follows which gives us the following set of equations, dX dN = − 3(w + 1)Q(X + 1)(Y + 1)(2(W + 1)X + 1) where It can be shown that the Eqs. (72)-(74) are equivalent to our equations (65)-(67) if the corresponding (71) is used properly. This for sure is a good consistency check that our equations are correct. The full study of the dynamical system (72)-(74) was carried out in [44] where 6 different kind of functions f 1 (T ) and f 2 (T ) were assumed. For some of those models, they found some critical points representing accelerating or decelerating solutions and also scaling solutions. For more details about all of this models and their dynamical analysis, see [44].

Exponential couplings for f (T, Lm) gravity
Now, let us study a new model where the function takes the following form where Λ is a positive cosmological constant. Let us take a look at this model further. In the limit where the argument is much less than one ( 1 Λ (T + L m ) ≪ 1), if one expands up to first order in the argument, the function becomes hence in that limit, one recovers the TEGR plus matter case with a cosmological constant. Therefore, the function (76) is an interesting model to take into account. An analogous model was proposed in [52], where instead of having T , the authors considered the scalar curvature R. The dynamical system of the former model was investigated in full detail in [42]. Under this theory, we directly find thatβ T LL = α T T T = 1 and by manipulating the definitions of the other quantities, β LL andβ T T L can be written as From the Friedmann constraint (63), L m = −2ρ and by using Eqs. (28) and (29) we directly find that for this model, the universe is always expanding as a De-Sitter one with a scalar factor being equal to It is interesting to see that actually this model is very different to its analogous in f (R, L m ). In the model f (R, L m ) = −Λ exp[−(R + L m )/Λ], it is not possible to directly find an unique scale factor which rules out the whole dynamic for the model. Hence, in that case, the dynamical system technique is very useful to check how the dynamics evolves on time. In our case, since the dynamics is always the same (described by the above equation) it is not important to study its dynamical properties since the solution is the standard De-Sitter universe. Therefore, the model described by an exponential coupling between L m and T as it is given by (76) mimics a De-Sitter universe.

Power-law couplings f (T, Lm) gravity
Let us consider another interesting new model that one can consider from our approach where the function takes the following form where ǫ is a constant and M is another constant which represents a mass characteristic scale. In this case, up to first order in ǫ, the expansion of the above function becomes so that since ǫ is assumed to be very small (comparable with T and L m ), the above model could represents a small deviation of the standard TEGR plus matter case. For this model we find that Similarly as we did in the previous section, from the Friedmann constraint (63), by replacing L m = −2ρ and by using Eqs (28) and (29) we find the following equation for the scale factor, which gives us two different types of scale factors. One can directly check that if ǫ = 0, the above equation is reduced to the standard TEGR plus matter case, namely 3H 2 = ρ. For the specific case where ǫ = −1/2, we must need ρ = 0, so that this special case is not a reliable model. There is no point on going further with the dynamical system of this model since the equation can be directly solved for the scale factor. The above equation depends on the power-law parameter ǫ. For negatives values of ǫ, the only possibility is that the second bracket is zero whereas for positives values of ǫ, there will be two kind of possible scale factor. This is again different as the case f (R) = M −ǫ (R + L m ) −ǫ+1 studied in [42]. Our model seems to be simpler than the former one due to the fact that T only contains derivatives of a(t) and not second derivatives as R.
Let us know explore what kind of solutions we have from our power-law model. The first type can be obtained by assuming that the first bracket is zero, which is only valid for ǫ > 0 giving us the following scale factor, where for simplicity we have chosen that the integration constant is zero. Let us clarify here that this scale factor will rule out the dynamic only for ǫ > 0. The scale factor must be real and positive so we must ensure that the imaginary term disappears. This is possible for some values of w. If one assumes that w > −1, for the solution a + (t), the state parameter must satisfy w + = 1 6k − 1 for any positive integer number k whereas for the solution a − , the state parameter must be w − = 11 6k − 1 to ensure a positive real value of a ± (t). Moreover, for these two solutions,ȧ ± > 0 andä ± > 0 for both w ± , so this solution could describe an accelerating expanding universe for those specific values of w ± . However, only the solution a − with w − = 11/6 ≈ 1.8 represents power-law expanding accelerating universes without evoking exotic matter. Additionally, Eq. (83) can be solved by letting the second bracket equal to zero, which is valid for all ǫ = −1/2, yielding a(t) = 3ρ 0 4(2ǫ + 1) where again for simplicity we have assumed that the integration constant is zero. This solution is very similar to (84) but now the parameter ǫ plays a role in the dynamics of the universe. Let us again consider the case where w > −1 for studying this solution. For ǫ > −1/2, the scale factor and its derivatives are always positive so that the scale factor will mimic a power-law accelerating universe. For ǫ < −1/2, we need to impose that otherwise the scale factor would be negative. Moreover, all the derivatives of the scale factor would be also positive if w satisfy the above condition .Hence, only special cases of w will give rise to viable models when ǫ < −1/2. Further, all those models are in the regime −1 < w < 0 which represents exotic kind of matter. Thus, cases with ǫ < −1/2 needs exotic matter to represent accelerating expanding universes. Additionally, we can conclude that for ǫ > 1/2, the power-law f (T, L m ) will mimic power-law accelerating universes without evoking exotic matter. In this section we will study the case where the function takes the following form where C 1 is a constant and the functionf (B, L m ) depends on both the boundary term and the matter Lagrangian. The first term represents the possibility of having TEGR (or GR) in the background when we set C 1 = 1. If this term does not appear in the function, it is not possible to recover GR since one cannot construct GR fromf (B, L m ) gravity. This kind of theories have not been considered in the past, but there are some studies for the specific casẽ f (B, L m ) = f (B) + L m , which is known as f (B) gravity [35,36,54]. The full dynamical system (45)-(54) is simplified since x 1 = x 3 = x 4 = α ≡ 0 which implies that Eqs. (45), (47), (48) and (53) are also automatically zero. Hence, in our variables, this theory is a 5-dimensional dynamical system given by where for simplicity we have introduced the following quantities Let us now concentrate on a specific model based on the boundary term non-minimally coupled with the the matter Lagrangian where the function takes the following form where C 1 is a constant and f 1 (B) and f 2 (B) are functions which depends on the boundary term B. This case is analogous to the one studied in Sec. IV C 1, but the dynamical system is more complicated to deal since it is a 5 dimensional one. The aim of this section is to study some specific cases that can be constructed from the above model. Let us further study the case where the functions are a power-law type given by where C 3 , C 4 , C 5 , q and s are constants. Since we are interested on studying non-trivial couplings between B and L m we will assume that C 3 = 0. We directly find that β LL = β BBL = 0. It can be proved that for this model, the dynamical system can be reduced from 5 dimensional to a 4 dimensional one. For this case, the dynamical system is difficult to study. However, if one assumes that the exponents are related as the system becomes easier to work since it becomes a 3 dimensional dynamical system. Then, we will split the study depending on different cases which depends on the constants. Let us first study a very special case where q = 1 in (95) giving us a linear coupling between the boundary term B and the matter Lagrangian L m . This model will depend on the power-law parameter s and also on the constants C 3 , C 4 and C 5 . In this model, one can relate two dynamical dimensionless variables with the other ones making it a 3 dimensional one. In our case, we will replace θ and y 1 as follows Thus, it is possible to replace the above equations into Eqs. (90)-(91) in order to reduce the dimensionality of the dynamical system for this model. By doing that, we find that the model only has one critical point given by P : x 2 , y 2 , φ = 0, s 6(s − 1) , 0 , which depends on the power-law parameter s. The case s = 1 can be discarded since a linear combination of the boundary term does not affect the field equations. It is easily to see that the effective state parameter for this critical point is always −1, hence this critical point always represents acceleration. To find out about the stability of this point, one needs to check the eigenvalues evaluated at P . There are three different eigenvalues given by One can directly see that when 1 < s ≤ 8/7 and w > −1, the critical point P is stable.
Let us now consider the case where C 5 = 0 in (95). This model represents the case where f 1 (B) = 0. Let us also assume that q = 1 to do not have the same model as the previous case. For this case, it is possible to relate the terms β BBL and β BB with the dynamical variables as follows Moreover, the dynamical system can be reduced from 5D to 3D since some of the variables are directly related, namely By replacing (101) and (102) in the dynamical system (88)-(92) we find that this system is reduced as follows dx 2 dN = θ(4φ − 2) + 3 4(w + 2)φ 2 − (3w + 5)φ + w + 1 + x 2 (2θ + 6(3w + 5)φ − 9(w + 1)) + 6(w + 1)x 2 2 2φ , (103) This dynamical system has only one critical point given by where we have assumed that q = (1 − w)/2. This point of course depends on the parameters q and w. For this point, there is acceleration when where we have assumed that w > −1. Further, for the dust case w = 0 we can see that this point requires q > −1/2 to represent and accelerating universe. It is also possible to check that there are three Eigenvalues associated with this point. Those Eigenvalues are very long to present here but Fig. 1 represents a region plot where the point is stable. Note that besides of the values of q and w, this point is never unstable. Let us now assume the case where q = 1 − s and the constant C 5 = 0 which is a more generic model which has an additional boundary power-law contribution. As we have studied in the previous section, we can again reduce the dynamical system as a 3-dimensional one. However, the models is much more complicated than the previous two models. The dynamics of the model highly depends on the parameter s. We can relate the termsβ BBL and β BB with the dimensionless variables but now those quantities are very long for a generic s. Moreover, those terms make the dynamical system very long and difficult to treat for any s. One can also relate two dimensionless variables with the other ones. In this case, we will choose to work with the variables (y 2 , θ, φ) since the dynamical system is slightly easier to work with them. The variables x 2 and y 1 are then given by It is possible to write down the dynamical system for any generic s but it is very long a cumbersome to present it here. Moreover, the critical points highly depend on the parameter s and it is not possible to obtain all the possible critical points for any arbitrary s. Hence, we will only study some particular models. We will concentrate only on models with integer values of s. Table I represents various models with their critical points, effective state parameter and their acceleration regime. In general, for all the critical points for those models, it is possible to have acceleration for the dust case w = 0. It is also important to mention that for s ≥ 2, all the models have only one critical point with the possibility of describing acceleration depending on the state parameter w. It can be proved that the critical points in the models s = 3, 4, 5 (there is only one critical point for each model) are always saddle points. The model s = 2 can be either a saddle point or an unstable point. Hence, for all of positive models of s, the critical points cannot be stable. When negatives values of s are considered, the system becomes more complicated. For s ≤ −3, the dynamical system becomes highly complicated to analyse. For the case s = −1, the critical point is either a saddle or s Model (y2, θ, φ) weff Acceleration an unstable point so it cannot be stable. Moreover, for the dust case (w = 0), the critical point for the model s = −1 is always unstable spiral. For the case s = −2, there are two critical points P 1 and P 2 (see Table I). The critical point P 2 is either a saddle point or stable whereas the point P 1 is always a saddle point. Fig. 2 represents the regions where the point P 2 is stable. It is important to mention that only the term C = C 3 C 5 appears in the Eigenvalues so that it is possible to make 2D region plots for the model. In this figure, it was considered the case C 1 = C 4 = 1 which is equal to consider the standard General Relativity plus matter model in the background.  Table I) is stable. The point is never unstable. All the blank regions represents the regions where the point is a saddle point.

V. CONCLUSIONS
In this work we have presented a new modified theory of gravity based on an arbitrary function f which depends on the scalar torsion T , the boundary term B and the matter Lagrangian of matter L m . Different kind of modified theories of gravity can be recovered from this theory. The incorporation of B in this function is with the aim to have the possibility to recover and connect standard metric theories based on the curvature scalar. This is possible since R = −T + B, so that it is possible to recover the generalised curvature-matter Lagrangian coupled theory f (R, L m ). Fig. 3 shows the most important theories that can be constructed from our action. The graph is divided into three main parts. The left part of the figure represents the scalar-curvature or standard metric theories coupled with the matter Lagrangian. Different interesting cases can be recovered from this branch, such as a generalised f (R, L m ) theory or a non-minimally scalar curvature-matter coupled gravity f 1 (B) + f 2 (B)L m or just standard f (R) gravity. The entries at the middle of the figure represent all the theories based on the boundary term B and the matter Lagrangian L m . In this branch, new kind of theories are presented based on a general new theory where the term C 1 T is added in the model to have TEGR (or GR) in the background. The right part of the figure is related to teleparallel theories constructed by the torsion scalar and the matter Lagrangian. Under these models, a new general theory f = f (T, L m ) is highlighted in box, allowing to have new kind of theories with new possible couplings between T and L m . As example, in this paper we have considered theories with exponential or power-law couplings between T and L m . Under special limits, these theories can represent a small deviation of standard TEGR with matter with or without a cosmological constant. As special case, this theory can also become a non-minimally torsion-matter coupled gravity theory f = f 1 (T ) + f 2 (T )L m , presented previously in [52]. Thus, different gravity curvature-matter or torsion-matter coupled theories can be constructed. Some of them have been considered and studied in the past but others are new. The relationship between all of those well-known theories have not been established yet. From the figure, one can directly see the connection between modified teleparallel theories and standard modified theories. The quantity B connects the right and left part of the figure. Hence, the connection between the teleparallel and standard theories is directly related to this boundary term B. Therefore, one can directly see that the mother of all of those gravity theories coupled with the matter Lagrangian is the one presented in this work, the so-called f (T, B, L m ). In this work, we have also studied flat FLRW cosmology for the general f (T, B, L m ) theory of gravity. Explicitly, we have focused our study on the dynamical systems of the full theory. In general, the theory is very complicated to work since it becomes a 10 dimensional dynamical system. This is somehow expected since the theory is very general and complicated. Using the full dynamical system found for the full theory, we then study different special interesting theories of gravity. For the case f = f (−T + B, L m ) = f (R, L m ), it was proved that our full dynamical system becomes a 5-dimensional one. Moreover, we have proved how one can relate our dimensionless variables with the ones used in [42] giving us a possibility of checking our calculations. We have found that the dynamics of this model is the same as it was described in [42].
The case f = f (T, L m ) is also studied, where in general the dynamical system can be reduced to be a 3-dimensional one. This theory is analogous to f (R, L m ) but mathematically speaking, it is different. It is easier to solve analytically the flat modified FLRW for a specific model for the f (T, L m ) than f (R, L m ). Further, for the later theory, for the exponential/power-law curvature-matter couplings one needs to study the dynamical system to understand the dynamics. For the f (T, L m ) case, the exponential/power-law torsion-matter couplings are directly integrated, giving us a scale factor of the universe directly from the modified FLRW equations. Hence, one does not need dynamical system technique to analyse the dynamics of those two examples. Another special interesting case studied was f (T, L m ) = f 1 (T ) + f 2 (T )L m . The dynamical system for this case is reduced as a 2-dimensional one. We have proved that our dimensionless variables can be directly connected to the ones introduced in [44]. This also gives us a good consistency check that our full 10-dimensional dynamical system is correct mathematically, at least for those special cases. Thus, the dynamics of those models are consistent with the study made in [44].
Finally, we have also studied the dynamics of modified FLRW for C 1 T + f (B, L m ) gravity using dynamical system. The dynamical system for this case becomes a 5-dimensional one, exactly as the f (R, L m ) case. The dynamics for this model is more complicated than f (T, L m ). This is somehow expected since B contains second derivatives of the scale factor and T only contains first derivatives of the scale factor (see Eq. 5). Further, R also contains second derivatives of the scale factor, exactly as B, so it is not so strange to see that the dimensionality of the dynamical system of f (R, L m ) is the same as C 1 T + f (B, L m ). Under the boundary-matter coupled model, we have studied a specific case where the matter Lagrangian is non-minimally coupled with B as f 1 (B) + f 2 (B)L m . By assuming some power-law boundary functions f 1 (B) = C 5 B s and f 2 (B) = (C 4 + C 3 B q ), we analysed the dynamics using dynamical system techniques. In general, the dynamical system for this power-law couplings are 4 dimensional but for the specific case where q = 1 − s, becomes a 3 dimensional one. Thus, we have analysed this model depending on three different limit cases: (i) q = 1, (ii) C 5 = 0 and lastly the case (iii) C 5 = 0, q = 1 − s. In general, the dynamics of all of these models are similar. As we have seen, mainly only one critical point is obtained for mainly all of them. The stability of those points were also studied, showing the regions where the critical points become stable.
As a future work, it might be interesting to study further other models that can be constructed from the full theory. In principle, one can use the same 10 dynamical system that we constructed here, and then simplify it by assuming other new kind of couplings between T ,B or L m . In addition, one can also use the reconstruction technique to find out which model could represent better current cosmological observations. Further, we can also incorporate the teleparallel Gauss-Bonnet terms T G and B G to have a more general theory f (T, B, L m , T G , B G ) (see [55]) or even a more general new classes of theories based on the squares of the decomposition of torsion T ax , T vec and T ten (see [56]). Then, one can study the dynamics of the modified FLRW for this general theory. By doing all of this, it will give a powerful tool to determine which models are better describing the current acceleration of the Universe, or other cosmological important questions.

ACKNOWLEDGMENTS
The author would like to thank Christian Böhmer for his invaluable feedback and for helping to improve the manuscript. The author is supported by the Comisión Nacional de Investigación Científica y Tecnológica (Becas Chile Grant No. 72150066).