2024 Cosmological ﬂuids with boundary term couplings

Cosmological models can be studied eﬀectively using dynamical systems techniques. Starting from Brown’s formulation of the variational principle for relativistic ﬂuids, we introduce new types of couplings involving a perfect ﬂuid, a scalar ﬁeld, and boundary terms. We describe three diﬀerent coupling models, one of which turns out to be particularly relevant for cosmology. Its behaviour is similar to that of models in which dark matter decays into dark energy. In particular, for a constant coupling, the model mimics well-known dynamical dark energy models while the non-constant couplings oﬀer a rich dynamical structure, unseen before. We are able to achieve this richness whilst working in a two-dimensional phase space. This is a signiﬁcant advantage which allows us to provide a clear physical interpretation of the key features and draw analogies with previously studied models.


Introduction
Cosmology, the scientific study of the universe as a whole, has undergone remarkable advances in recent decades and General Relativity (GR) provides a good model to describe cosmological gravitational phenomena [1,2,3,4].On the other hand, open questions in cosmology remain, foremost among which are the dark energy and dark matter problems.The nature of dark energy, which is responsible for driving the universe's late-time accelerated expansion, is not well understood, and it often assumed to be a cosmological constant.Since the first observational evidence of an accelerated expansion [5,6] of the universe, a plethora of cosmological models to explain dark energy has emerged, for a review see [7].
The addition of a positive cosmological constant Λ to the Einstein field equations, originally introduced by Einstein [8] for his static universe, is one of the most straightforward candidates for dark energy.This paves the way for the Λ Cold Dark Matter (ΛCDM) model.However, the ΛCDM model fails to explain why the inferred value of Λ is so small compared to the vacuum energy density expected from particle physics [9].It is also unclear why its value is comparable to the matter density today.This constitutes the so-called coincidence problem [10,11].
One way to begin to address this issue is to allow for a dynamical cosmological constant [7], that is, to introduce some dynamical field able to reproduce the late-time acceleration behaviour and mimic the properties of the cosmological constant.The simplest such model is a canonical scalar field φ with flat potential V (φ), which drives the accelerated expansion of the universe.Any model of this type is referred to as quintessence [12].Scalar fields play a major role in modern cosmology as they are also able to drive inflation, the early-time epoch of accelerated expansion [13,14,15].Scalar field models have also been used as candidates for dark matter models, see [16].We are primarily interested in scalar fields as models to drive a period of accelerated expansion, at both early and late times of the universe's evolution.
Another approach is to consider dark energy as evidence for the incompleteness of GR and, hence, seek extensions or modifications of GR [17,18,19].Several models to describe the dark energy interaction with dark matter have been proposed [20,21,22,23,24,25,26,27] while some authors, e.g.[28], have emphasised a strong distinction between modified theories of gravity and dark energy models.
From a theoretical point of view, in GR, one usually restricts the Lagrangian to a linear function of the Ricci scalar, minimally coupled with matter.However, there is no reason, a priori, to assume such a restriction.So, one can modify the gravitational part of the action to allow non-linear corrections to the Lagrangian [29,30], this is the general approach followed by f (R)-theories of gravity [31,32].Some other extensions of GR increase the number of spacetime dimensions or introduce non-minimal matter couplings to boundary and topological terms [33,34,35,36,37,38,39,40,41,42].These are terms in the Lagrangian that describe how matter couples to geometrical quantities.Non-minimally coupled terms involving curvature vanish in the limit of special relativity.
We will follow the approach suggested in [43,44], and, therefore, build models of quintessence interacting with dark matter.This involves introducing couplings at the level of the action which characterise both quintessence and dark matter [45,46].In particular, we extend the Lagrangian formulation by Brown [47] which describes a perfect fluid [48].Motivated by [43,44], we introduce new couplings containing a boundary term and a pseudovector related to the boundary term.
In [40,49], a dynamical system analysis where teleparallel quintessence is non-minimally coupled to a boundary term is presented.In the same spirit, we study the background cosmology within this framework and apply dynamical systems tools to investigate the dynamics of the different models.Our ultimate goal is to examine the behaviour of these dark energy models.Dynamical systems theory has emerged as a vital tool in cosmology [50], and has been employed successfully to study modified theories of gravity in the cosmological context [51,52,53].
The paper is organised as follows.In Section 2, we present the Lagrangian formulation of our models, and, in Section 3, we focus on the cosmological field equations and introduce the cosmological variables which will be used for our analysis.Section 4 contains our analysis of the dynamical systems for a constant interaction term, and we highlight the analogies with previous models [54].In Section 5, we present the rich and novel dynamical structure in the case of a non-constant interaction term, which, for some choice of the parameters, features both early-time and late-time accelerated expansion.In Section 6, we discuss our results and suggest potential directions for future investigations.
Notation and conventions.Unless otherwise specified, we employ standard relativistic notation throughout.The signature of the metric tensor g µν is assumed to be (−, +, +, +), Greek indices are space-time indices taking values in {0, 1, 2, 3}.The coupling constant appearing in the Einstein field equations is denoted by κ = 8πG/c 4 , where c is the speed of light and G the Newton's gravitational constant.We use natural units, with c = 1 and G = 1.A dot denotes differentiation with respect to cosmological time, a prime denotes the derivative with respect to the argument, or in the case of the dynamical system equations a prime denotes a derivative with respect to the logarithm of the scale factor log(a).
2 Lagrangian formulation and field equations

Gravitational and fluid action
It is well known that a total derivative term can be isolated from the Ricci scalar, yielding the Gamma squared action.This action also gives rise to the Einstein field equations when variations with respect to the metric are considered.However, the underlying Lagrangian is no longer a coordinate scalar as it differs from a true scalar by the total derivative term.We prepend the word 'pseudo' to highlight quantities which appear to be scalars or tensors but are not invariant under general coordinate transformations.As shown in [42], this allows one to write the Einstein-Hilbert action as The bulk term G is defined as and the boundary term B is given by where we have introduced the boundary pseudovector B σ given by We note that G and B are pseudoscalars.By construction, the bulk term G is quadratic in the Christoffel symbols and hence the action is called the Gamma squared action or, sometimes, the Einstein action, to distinguish it from the Einstein-Hilbert action.Recent progress was made in [42,55] on constructing modified theories of gravity based on this decomposition.These models can be linked naturally to a variety of other modified gravity models, either within the context of GR or in the metric-affine framework.The Christoffel symbols are usually interpreted as the gravitational field strengths.We can motivate this by recalling that they contain the first partial derivatives of the metric, which represent the gravitational potentials.The bulk term is thus quadratic in the field strengths, similar to other field theories, like Yang-Mills theories or elasticity theory.This analogy provides the primary motivation for splitting the Ricci scalar as in (2.1).This split naturally yields a boundary term which could be coupled to other fields present in the model.Couplings of this type are interesting as they are purely geometrical and thus have no direct links with classical physics.This is similar to Brans-Dicke theories, where a scalar field is coupled to the curvature scalar.By isolating the bulk and boundary terms, we can therefore consider more intricate couplings involving those two parts, which make up the curvature scalar.
Matter is introduced into the theory by considering a total action of the form where S matter is the matter action.This gives rise to the energy-momentum tensor Generally, a matter action would depend on matter fields, and variations with respect to those matter fields yield the equations of motion of the matter component.As we will see in the following, this is non-trivial if one wishes to model relativistic fluids using the variational approach.
Brown [47] introduced a Lagrangian formalism for relativistic perfect fluids based on the Lagrangian (density) given by where • n is the particle number density • s is entropy density per particle • ρ(n, s) is the energy density of the matter fluid, a function of n and s • J µ is the vector-density particle-number flux, which is related to n by where U µ is the fluid's 4-velocity satisfying U µ U µ = −1 • ϕ, θ, and β A are all Lagrange multipliers with A taking values in {1, 2, 3}, and the components α A are the Lagrangian coordinates of the fluid.
The independent dynamical variables of the Lagrangian (2.8) are g µν , J µ , s, ϕ, θ, β A , and α A .We note that, in this approach, the pressure of the fluid p is defined as which is consistent with the first law of thermodynamics.

Total action and interaction terms
We can now set up the total action which contains gravity, a fluid, a scalar field φ, and an interaction term.This means where L φ is the scalar field Lagrangian (density) given by with given scalar field potential V .The Lagrangian (density) L int is an interaction coupling term, which allows us to couple the fluid to the scalar field.Due to the presence of the various independent variables in Brown's approach, one can propose some types of coupling terms which do not exist in other settings.Moreover, such terms have no natural special relativistic analogue, making this potentially interesting in the context of cosmology.In previous work [43,44], one of the authors proposed interaction terms of the form f (n, s, φ) and f (n, s, φ)J µ ∂ µ φ.These models gave rise to some unexpected dynamics.In particular, as we wish to take into account boundary terms, we identified the following terms as the suitable possibilities Note that B µ was defined in (2.4).Depending on the specific interaction term chosen, one should also note that the physical dimensions of f differ in the different couplings.
Coupling (b) is very restrictive in the context of cosmology.We find that the consistency of the cosmological equations implies that f is proportional to n, thereby eliminating the scalar field from the coupling.Consequently, we find equations which largely coincide with the standard cosmological equations, and the model does not exhibit novel behaviour.
For the remainder of this paper, we will consider the interaction term (c) and L int will denote this interaction Lagrangian.This term was found to have behaviour relevant to cosmology, and gave rise to manageable cosmological equations.In principle, our analysis can be repeated for term (a), and potentially for more complicated terms.For example, f could contain an explicit dependence on G, or the Ricci scalar, or there could be higher order couplings containing terms like (B µ J µ )(B ν ∂ ν φ), etc.

Variations and field equations
We begin with the variations of action (2.11) with respect to the the fields ϕ, θ, β A , and the Lagrangian coordinates α A , respectively.This yields δϕ : J µ ,µ = 0 , (2.13) δθ : (sJ µ ) ,µ = 0 , (2.14) ) These equations are independent of the gravitational action and are also independent of the interaction term.Next, variations of (2.11) with respect to the entropy density, s, give where the final term depends on the choice of f .Variations with respect to J µ yield Again, we have one term which depends on the coupling.Variations with respect to the scalar field φ yield a modified Klein Gordon equation where := ∇ µ ∇ µ .
Finally, variations with respect to the metric tensor yield the Einstein field equations where G µν is the Einstein tensor and Both are the standard forms of the energy-momentum tensors of a perfect fluid and a scalar field, respectively.The energy-momentum tensor related to the interaction term is more complicated and is given by The second line of this tensor appears due to variations of the boundary pseudovector with respect to the metric.This requires integration by parts and thus leads to the second derivative terms of the scalar field.We note that one could rewrite the partial derivatives of the metric determinants using the Christoffel symbols.However, for the purposes of this work, this would not introduce additional insights.
We motivated the introduction of L int as a new interaction term which could model an interaction between the fluid and the scalar field.However, one can adopt a different interpretation, namely to view (2.23) as an independent fluid of unusual form.

Cosmological field equations
In this section, we provide a brief overview of the necessary background material required to study the cosmological field equations of our coupled models.We do this via a dynamical systems formulation, which has proved to be a powerful tool when studying cosmological equations.
In line with current observational evidence [2,56,57], let us begin with the homogeneous, isotropic, and spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) line element where a(t) is the scale factor and N (t) is the lapse function.For all models under consideration, we will be able to set N = 1, which simplifies the cosmological equations further.In this case, t is cosmological time.In this cosmological setting, (2.20) yields the cosmological Einstein field equations given by and (2.19) leads to the modified Klein-Gordon (KG) equation Here the dot denotes differentiation with respect to cosmological time, and we remark that ρ and p are the energy-momentum density and pressure of the fluid, respectively.
A direct, but lengthy, calculation verifies that the three equations (3.2)-(3.4)imply the fluid's energy-momentum conservation equation ρ + 3H(ρ + p) = 0.This is a non-trivial result which is, perhaps, unexpected given that the coupling contains an unspecified function.Let us also note that the only dependence on the scale factor a(t) in the field equations is via the Hubble function H and its derivative.These equations feature both first and second derivatives of the scalar field, φ.However, following [54], one can introduce a new variable which depends on the first derivative of the scalar field, leading to field equations which are first order.In short, this is the key idea behind the dynamical systems formulation.
If we assume that the potential is non-negative, we can introduce the well-known dimensionless variables, first proposed in [54], We restrict to the case that H > 0, i.e. that the universe is expanding (choosing H < 0 would correspond to a contracting universe).It follows that the variables y and σ are non-negative.
In line with previous studies, we assume V has the exponential form where V 0 > 0 is a constant and λ ≥ 0 is a dimensionless parameter.We note that this form for V is invertible, which will allow us to view φ as a function of V .This potential is most convenient as the exponential form allows one to close the autonomous system of equations without the introduction of an additional variable.
When the FLRW metric is considered in (2.13) and (2.14), one immediately finds that the entropy density s = s 0 is a constant.Consequently, the coupling function is of the form f (n, φ) only.Moreover, (2.13) and (2.14) also imply that the particle number density is n = n 0 a −3 , where n 0 is a constant, which is expected.
Going back to Brown's formulation (2.8), we have that the energy density is a function of n, since the fluid's entropy s is constant, thus ρ = ρ(n).On the other hand, in standard cosmology, it is customary to assume a linear equation of state of the form p = wρ.We will now demonstrate that, given the definition of pressure in (2.10), this is equivalent to the assumption that the density is a power of the particle number density.To begin with, let us consider ρ = n w+1 , for some w, which implies For the matter dominated case, w = 0, this gives that p = 0. On the other hand, integrating (2.10) with the assumption that p = wρ implies ρ = n w+1 .Dividing equation (3.2) by 3H 2 and using the variables (3.5), one obtains We note that f must be chosen to have the same dimensions as κ −1/2 to ensure that this equation is consistent.As it is derived from the Friedmann constraint equation, we will generally refer to it as the constraint equation.This is motivated by the fact that it is an algebraic relation between all the variables, which implies that the variables are not all independent.We finish this Section by noting that for f (n, φ) = 0, we retrieve the model studied in [54], which we can view as our baseline model.When interpreting our results, we draw analogies and highlight differences with this baseline model.In that work, the constraint equation (3.8) is solved for the matter variable σ which is then eliminated from the other equations, reducing the system to two differential equations and we follow the same approach here.

Constant interaction
To begin our study, we consider perhaps the simplest non-trivial model, where the coupling function is a constant for some constant k.This model shares some similarities with [54] and is an ideal prelude to the study of more complicated models.

General properties and dynamical systems formulation
Let us start with the Klein-Gordon equation (3.4), which simplifies to where we introduce the quantity Q to match previous work on dark sector couplings [45].The energy density and pressure of the scalar field are given by ρ φ = φ2 /2 + V and p φ = φ2 /2 − V , respectively.This allows us to re-write equation (4.2) in the well-known form hence Q can be re-expressed as where q = −1 − Ḣ/H 2 is the standard deceleration parameter.It is well known from dark sector coupling models [45,46] that Q > 0 means an energy transfer from dark matter to dark energy and Q < 0 a transfer in the opposite direction.
For k > 0, equation (4.4) implies that an epoch of accelerated expansion, q < 2, gives a positive coupling, leading to energy going into the scalar field.In turn, this leads to an epoch of further acceleration and can be seen as a self-reinforcing effect.The above argument is reversed for f < 0 (i.e.k < 0).Given that the late-time universe is dark energy dominated while the early universe contains considerably more dark matter than dark energy, it is reasonable to consider f > 0 (i.e.k > 0) and it will turn out that such models indeed evolve into epochs of late-time accelerated expansion.
Next, we consider a fluid with equation of state p = wρ, which, as discussed at the end of Section 3, is equivalent to setting ρ = n 1+w .The constraint equation (3.8) then reads The quantity σ 2 is the relative energy density of matter, sometimes denoted by Ω m when discussing explicit cosmological models.For a scalar field, it is helpful to introduce the equation of state Here, w φ ∈ [−1, 1] and we get w φ = −1 when φ = 0, as is expected for dark energy.The energy density of the scalar field is given by Ω Hence, (4.5) can also be written as At this point, it is clear that one can introduce improved variables by completing the square of the x-term in (4.5).Namely, we write and now divide by the new right-hand side so that we arrive at where These variables will prove particularly useful for our subsequent qualitative analysis.Using equations (3.2)-(3.4),one can obtain the acceleration equation which can be integrated to find a(t) at any given fixed point (X 0 , Y 0 ).The right-hand side, at a fixed point, is constant.If this constant is non-zero, it is straightforward to show that the scale factor evolves as a power law in cosmological time, that is, a ∝ (t − t 0 ) µ , where µ is that power.We therefore have that µ is given by Here t 0 is an integration constant.When the right-hand side of (4.12) vanishes at some fixed point, the scale factor a(t) evolves exponentially.This corresponds to H being constant at this point, that is, a universe undergoing a de Sitter expansion.
It can be useful to define the total energy density and total pressure of the cosmological model where we set ρ int = −kH 6/κ + φ and p int = k φ/ √ 6κ, as suggested by (3.2) and (3.3).This naturally leads to the effective equation of state parameter w = p/ ρ.For power law models, this effective equation of state parameter is directly related to the power µ, and one has We note that the power µ, the effective equation of state parameter w, and the deceleration parameter q, all encode the same physical information.Similar to previously studied models, the positivity of the matter variable and equation (4.10) imply that 0 ≤ Σ ≤ 1, and hence 0 ≤ X 2 + Y 2 ≤ 1. Together with the fact that Y ≥ 0, since we are considering an expanding universe, this means that the phase space for the variables X and Y is a semicircle of radius one.
We are now ready to state the dynamical equations of the system, using the convenient variables defined in (4.11).This leads to two independent equations Here a prime denotes a derivative with respect to the logarithm of the scale factor log(a).One can now follow the standard dynamical systems approach to study this system, for a review see [50].We begin with the fixed points of (4.17)-(4.18).We note that these are two polynomial equations of degree three, meaning that one could find up to nine real distinct critical points, by Bézout's theorem.If Y = 0, the second equation is automatically satisfied, and this leads to the solutions X ∈ {−1, 0, +1}.Next, excluding Y = 0, one notes that Y appears only as Y 2 in the equations, meaning that there are up to four more solutions.Two of these are at negative values of Y , which we exclude, again because we are considering an expanding universe (Y ≥ 0).Assuming λ > 0 and −1 ≤ w ≤ 1, we obtain a total of five critical points, shown in Table 1.
Note that Point B is always located on the boundary of the phase space while Point C is generally inside the phase space, if it exists.For the special value the lower existence bound, Point C is also on the boundary.Next, one needs to study the eigenvalues of the stability matrix at each of the critical points.For more details about their classification, see Appendix A. For the first four points, O, A ± , and B, these are given in Table 2.Note that we will discuss the occurrence of possible zero eigenvalues separately to keep the discussion more straightforward.For example, one may immediately note that the choice w = 1 implies at least one zero eigenvalue for the Points O and A ± .The final critical point, C, is more difficult to study as the eigenvalues are much more involved.They are the solutions of the characteristic polynomial in ξ Solving this quadratic equation is easy, however, the explicit solutions do not offer much insight given that they contain three free parameters.For concrete parameter choices, we discuss this point in more detail below.One easy result to extract is the sum of the eigenvalues ξ 1 and ξ 2 at this point, that is, As this number is negative for w < 1, this point cannot have two positive eigenvalues and therefore will have at least one stable direction.This implies that Point C is a saddle point, stable node, or stable spiral.
There will be many parameter choices resulting in zero eigenvalues, w = 1 being the obvious one.However, the choice kλ = √ 6(1 + w) would also give a zero eigenvalue for Point O.The stability analysis of such points requires techniques beyond linear stability theory.These are well known and their applications in cosmology were discussed in [50,58,59].However, for the purposes of this work, we will assume a matter dominated universe w = 0 and employ linear stability theory.We note that our analysis can also be performed for the radiation dominated case w = 1/3 where one finds qualitatively similar results.

The matter dominated case
In what follows, we set w = 0. Points O, A ± and B are independent of w and all results discussed above apply.The location of Point C, if it exists, depends on w and so do its corresponding eigenvalues.We now outline some physical properties of the critical points of the system, with the values of effective equation of state parameter w and the deceleration parameter shown in Table 3.
In order to analyse the stability of the fixed points, we look at the different regions in the (k, λ)plane, see Fig. 1, and we recall that λ > 0. First, we remark that the fixed Points O and A ± exist for all values of λ and k.Moreover, there are four distinct regions of values of k and λ, which yield different stability properties of the critical points, and, hence, different cosmological phenomena.We discuss these four different cases and comment on their suitability as a cosmological model.2.

I II III IV
Region I.For values within Region I, there are only three critical points, namely O and A ± .In particular, O is a stable node, A − is an unstable node, and A + is a saddle.Since O is the only attractor of the system, all trajectories will eventually approach it.The Point A − can be thought of as the past-time attractor, in the sense that all trajectories would start at A − .Lastly, some trajectories are attracted towards A + , but are eventually repelled and move towards O.This case is not of physical interest.We do not show a phase space diagram.
Region II.In Region II, Point B does not exist.We therefore have four critical points: the unstable node at Point A − , the stable node at Point C, and the saddle Points O and A + .We note that here, Point C represents the scaling solution [60] as the effective equation of state parameter matches the matter one ( w = w = 0).Hence, the universe expands as if it was completely matter dominated despite the scalar field's influence, according to (4.13).We note that this is not an accelerated expansion but this solution is of physical relevance for the coincidence problem.The type of dynamics is illustrated by Fig. 2, where we set k = 1 and λ = 2.We remark the analogy with one of the cases discussed in [54], however we also point out that, in our example, no acceleration region is present.Region III.In Region III, there are five critical points in the phase space.Points A − , O, and C, still behave as an unstable node, a saddle point, and a stable node, respectively.Point A + is now an unstable node.Point B exists and is a saddle point.This is shown in Fig. 3, where k = 1/2 and λ = 3/2.Point C always lies outside the acceleration region, so it does not represent a late-time inflationary solution.This is, again, a scaling solution.We highlight the analogy with another case discussed in [54].All trajectories connect Points A ± to Point C, with the exception of the orbits along the boundary.
Region IV.In Region IV, there are again four fixed points since Point C lies outside the physical space.Here, Point A − is always an unstable node and can be seen as the past attractor.Similar to Region III, A + is an unstable node.Point O is a saddle point, whereas Point B is a stable node and therefore the late-time attractor.We note that here Point B lies within the region of accelerated expansion, hence we are in the presence of a cosmological solution with accelerated expansion.This is illustrated in Fig. 4. Once again, we emphasise the analogy with one of the cases discussed in [54].

A non-constant interaction model 5.1 Equations of the model
We are now considering a model with a non-constant interaction term.As we wish to exploit dynamical systems techniques without increasing the number of independent variables, we consider an interaction of the form [61,62] where α is a fixed power.Let us make the following observations to motivate this particular choice for f .In cosmology, s = s 0 and n = n 0 a −3 .Moreover, if we consider the linear equation of state 2) The dynamics of the system depends on the parameters w, λ, k, and α.We note that one could consider the limit α → 0 and recover equation (4.5).
It is clear that larger (integer) values of α can make the study of this system difficult, since the constraint (5.2) would become a polynomial equation of high order.At the same time, even the values α = ±1 introduce challenges as one has to deal with cubic equations.In fact, the two simplest cases that can be studied explicitly, without introducing further complications, are α = ±2.In what follows, we consider α = 2 and w = 0, and note that the radiation dominated case (w = 1/3) leads to broadly similar results.
When α = 2, the constraint (5.2) can be written as allowing us to eliminate σ from the equations, so that the dynamical system remains two-dimensional.Moreover, as σ 2 ≥ 0, we have This gives rise to the physical regions of the phase space 1 − x 2 − y 2 ≥ 0 and y 2 − kx > 0 , (5.5) 1 − x 2 − y 2 ≤ 0 and y 2 − kx < 0 . (5.6) Since the non-constant coupling leads to a considerably more complicated dynamical system, we restrict our study to the matter dominated case w = 0.The dynamical equation for x is given by where the functions A and B are defined by (5.9) Similarly, the y equation reads where (5.12) Here a prime denotes a derivative with respect to the logarithm of the scale factor log(a).This way of writing the dynamical equations, namely isolating the terms in powers of k, is useful as it allows us to consider the limit k → 0 easily.In that case these equations reduce to those of [54].
We remark that the acceleration equation, which follows from (3.2) and (3.4), is where the functions E and F are given by (5.15)

Critical points and stability
To find the critical points, we need to solve the equations x ′ = 0 and y ′ = 0 simultaneously, which is a non-trivial task since both numerators are polynomials of degrees seven, giving up to 49 roots.Many of those will lie outside the physical phase space, while others will come in complex conjugate pairs which also have no physical significance.At this point, it is not clear how many physical critical points this system will have for arbitrary λ and k and hence one has to investigate the system carefully to extract them.One way to find the critical points is to draw inspiration from the previous model.For example, setting y = 0 in (5.10) leads to y ′ = 0, while setting y = 0 in (5.7) means that x ′ = 0 simplifies to This yields the first set of critical points (−1, 0), (0, 0), (1, 0), and ( 3/2 /λ, 0).Secondly, we investigate critical points on the unit circle.By substituting x = cos θ and y = sin θ into (5.7) and (5.10), at a critical point we obtain This gives another critical point at (λ/ √ 6, 1 − λ 2 /6).Lastly, one can verify that setting x 0 = √ 3 ( √ 2λ) in the dynamical equations gives four additional solutions, other than y 0 = 0, which are We are not able to find other critical points in the physical phase space, either analytically or numerically.The critical points discussed above are summarised in Table 4, together with their corresponding value of the effective equation of state parameter and the value of the deceleration parameter.Note that there will be parameter regions where the critical points with y-coordinate y ± , called D ± , may not exist or where only one of these exits, see Fig. 5.
Table 4: Critical points of the dynamical system (5.7) and (5.10), for which an explicit expression could be found.We are now ready to investigate the stability of the critical points.This is straightforward for the Points O, A ± , B, and C. The result are collected in Table 5.

Point Eigenvalues Classification
Table 5: Stability properties of the critical points assuming λ > 0 for the non-constant interaction model.
For the Points D ± , the closed form expressions for the eigenvalues are very long and do not offer physical insight.However, when presenting specific cases, we give numerical values for the eigenvalues and discuss the various critical points in more detail.
Table 4 suggests that Point C is of particular interest to us.This point has an effective equation of state parameter w < −1/3 if λ < 15/2, and it is not located on the boundary of the phase space.For such a choice of λ, we note 15/2 > 3/ √ 2, which means this point will be an attractor of the dynamical system.In turn, such a model will naturally give rise to a period of late-time accelerated expansion.Moreover, since 15/2 > √ 6, Point B will not exists in this case.The physical phase space for these models is delineated by the upper semicircle of unit radius centred at the origin and the parabola x 2 = y/k, as described by (5.5) and (5.6).Subsequent figures will make clear which regions form the physical phase space.For all k, the semicircle and the parabola will intersect, creating two regions which meet at a point, and which trajectories can traverse.Note that the intersection point is, in general, not a critical point.

Phase space diagrams and physical interpretation
Different choices of λ and k result in rather different cosmological models, since the number of critical points and their location vary significantly.Below we consider several cases which illustrate the diversity of the dynamical behaviour exhibited by the model.We select values of λ and k systematically, but do not necessarily include every possible scenario which could arise in these models.
Case (i).We begin with λ = 3/2 and k = 8, as shown in Fig. 6.In this case, Point C does not exist, however, the other six critical points do exist.The phase space contains a region of accelerated expansion and we note that only Point O is in this region.Point O is an early-time attractor of the phase space, hence, this point could correspond to an early-time universe undergoing accelerated expansion.We note that this is negative but not less than −1/3, and therefore, not accelerating.Depending on the initial conditions chosen, some trajectories will approach Point D − with w = 0, which is matter dominated.On the other hand, trajectories starting at A + will either also terminate at B, or reach D + .This latter point is a stable spiral with w = 0, and hence corresponds to a matter dominated universe.It is interesting to note that trajectories in this case can briefly go through a region of accelerated expansion before reaching D + .While these parameter values yield an interesting phase space with a rich structure, this specific model has limited applicability for modern cosmology, since the stable fixed points do not lie within the accelerated region of the phase space.
Case (ii).Next, we consider the case λ = √ 3 and k = √ 8, see Fig. 7.In this particular case, Point D − coincides with Point B, this is true for all k.As in the previous case, Point C does not exist in the physical phase plane.Point O, as before, is an unstable node and acts as an early-time attractor.Point A − is a saddle and corresponds to the to other possible early-time attractor.According to Table 5, we note that Point B has eigenvalues 0 and −3/2, which means that we are dealing with a non-hyperbolic point.This point is a centre and one can verify that it is unstable.While this can be shown rigorously, it essentially follows from the fact that trajectories near B move towards the attractor D + , which is Case (iii).We will briefly comment on the case where λ = 3/ √ 2 and k = 2/ √ 3.For this particular choice, Points C and D − do not exist, while the Point D + is located at the intersection of the two regions of the phase space.This case is mathematically quite interesting, however, less so from a physical point of view.Of mathematical interest are the following facts: D + is a critical point of the system, however, both the numerators and the denominators of (5.7) and (5.10) vanish while giving a finite limit.The stability matrix is singular at this point, meaning that linear stability theory cannot be used.We will not discuss this case further.
Case (iv).Next, we consider the case λ = k = 1, which turns out to be physically interesting as a cosmological model, see Fig. 8.There are two unstable nodes, Points O and C, which act as earlytime attractors.As in the previous cases, Point O corresponds to an early-time universe undergoing accelerated expansion.Trajectories starting near O will eventually leave the acceleration region and be partially attracted to the saddle Point A − , after which they will reach the late-time attractor, the stable node at B. This point is also in the accelerated region, which means that this model not only allows for early-time acceleration (inflation) but also for late-time accelerated expansion.The effective equation of state parameter at B is w = −2/3, as can be seen from Table 4.All trajectories starting out in the right part of the phase space will also be attracted to B, making this the global attractor of the system.We remark that, in this sense, the dynamical behaviour is similar to that shown in Fig. 4.

Case (v)
. We now present the case λ = 1 and k = 20, see Fig. 9.Here all the fixed points we found analytically exist within the physical phase space.This model not only has a rich dynamical structure, but is also of physical relevance.We have two early-time attractors, Point O in the acceleration region, similar to the previous models, and Point C. Notice that Point C always satisfies w > 1, and so is not of physical interest.We are therefore most interested in trajectories starting from Point O.These will initially move towards A − , before leaving the left part of the phase space.By doing so, they will enter the acceleration region and move towards B where w = −2/3.Other than the various complications introduced by the other critical points, and the more complicated phase space structure, the physical situation is again somewhat similar to those shown in Fig. 4 and Fig. 8. Let us also mention that Point D + represents scaling solutions as the effective equation of state parameter is zero, and the universe evolves as if it were only matter dominated while also containing the scalar field.
Case (vi).To complete this section, we consider a case where the sign of the coupling is negative.We set k = −1/4 < 0 and λ = √ 5, and the phase plane is shown in Fig. 10.This model displays significantly different features than the cases where the coupling is positive.Point O can still be seen as an early-time attractor in the acceleration region.However, depending on the chosen initial condition, trajectories will either terminate at Point D + or Point C. Such trajectories can come close to the saddle Point D − , where w = 1/2.However, the effective equation of state parameter is also quite large at the other two points, meaning that one cannot have a model with a late-time behaviour close to a matter dominated universe.None of the late-time attractors appear in the acceleration region.Moreover, the left-hand side of the phase plane indicates the presence of a critical point at infinity,

Conclusions and discussions
The entire field of cosmology has seen remarkable progress in recent decades, including the use of dynamical system techniques to study the background behaviour of cosmological models.These techniques offer a systematic approach to understanding the underlying dynamics, which allows us to investigate the suitability of such models as realistic approximations of the universe.Our analysis involves mapping the cosmological equations onto a phase space, a step which relies heavily on the choice of suitable variables.This is rather non-trivial as various different variables could be employed and there is no particular reason to prefer one set of variables over any other.We therefore work with those variables that are known to be well suited for our task, see the review [50].
One of the main motivations of this work was to study models derived from a variational principle, in particular, we used Brown's approach for the formulation of the perfect fluid Langrangian for the cosmological matter.This approach allowed us to introduce new coupling terms, including boundary term couplings which have not been studied before in this context using fluids.So far, we have considered an algebraic vector coupling of the form f (n, s, φ)B µ J µ and noted that, in cosmology, one obtains the highly restrictive condition that f is proportional to n.Therefore, such a coupling does not yield interesting phenomena.We therefore focused on a derivative type coupling, f (n, s, φ)B µ ∂ µ φ, motivated by previous work [44].In that previous work, the coupling f (n, s, φ)J µ ∂ µ φ led to a minor change of the phase space compared to [54]: the critical point at the origin moves along the x-axis, depending on the choice of parameters.Our new coupling allows for a significantly different dynamical behaviour with features unseen before.Of particular interest to cosmology are situations where the model evolves through two periods of accelerated expansion, which replicate inflation and dark energy in a single model.
To gain an initial understanding of the resulting cosmological model, we began by studying the constant interaction model.This displayed similar behaviour to the well-known exponential potential quintessence model [54], where no interacting term is present.In fact, through a carefully chosen change of variables, we were able to arrive at a phase space which mirrors the one studied by those authors.These results demonstrate that our model can be seen as a natural extension of previous work.
We then proceeded to consider an interaction of the form n α/(2(1+w)) V −α/2 .This choice was motivated by the fact that couplings of this form will not increase the dimension of the phase space: it will remain two-dimensional.The key advantage of this assumption is that one can directly compare results with many previously studied models.In particular, we focused on the matter dominated case, w = 0, and chose α = 2.This choice leads to a rich dynamical structure with several distinct scenarios that can be of physical relevance.We were able to obtain solutions with early-time inflationary attractors, as well as late-time acceleration.Our models also included scaling solutions, which have received recent attention, see [63], as they may help to resolve the Hubble tension, that is, the discrepancy between the value of the Hubble constant inferred from measurements of the early universe and those derived from more recent observations [64,65].
Our approach to constructing coupled models lends itself to a significant amount of further study.First of all, one can study the constant coupling model in the radiation dominated universe.Our preliminary work suggests that the results are qualitatively similar to the matter dominated case, which is why we did not include them here, for the constant coupling model.One could attempt to present a comprehensive study for all w, however, this would not be without challenge as the convoluted equations would make analysing the stability at fixed points difficult.
Regarding the non-constant coupling model we proposed, there are three obvious extensions to our work, namely the cases α ∈ {−2, −1, 1}.For integer values of α, the phase will remain two-dimensional.However, one encounters other challenges which can be seen in equation (5.2).For example, when setting α = −2 one should eliminate y from the equation instead of the matter variable σ.This is unusual and has rarely been considered in the past.As a starting point, one would have to go back to the baseline model [54] and study it using a different choice of variables.In this way, comparisons could be drawn.For large values of α, on the other hand, one is dealing with a polynomial of a high degree, which is difficult to handle.In such cases, the best way forward would be to eliminate the variable x.This is equally unusual, and has also not been considered in the past.Given the complexities of these models, we are not able to predict the qualitative features of the resulting systems.The shape of the phase space alone changes significantly when varying the parameter α.
On top of all of this, one could, of course, drop the assumption of working with an exponential potential.It would be interesting to study our coupled models for power-law potentials or others.
Most of what has been done here will have to be re-investigated from scratch.For example, it is not even clear which types of couplings will admit a two-dimensional phase space.One would expect the constant coupling models to be similar to the uncoupled models, however, we refrain from speculating about results beyond this most basic of statements.
To determine the behaviour of trajectories near those fixed points, we can linearise the system around its critical point, by using a Taylor expansion for f in the neighbourhood of the fixed point.The dynamics of the linearised system are qualitatively equivalent to the original system.The eigenvalues of the matrix ∇f (x * ), known as the Jacobian matrix or stability matrix, contain the information about the local behaviour of f near x * .One generally speaks of stability or instability: • if all eigenvalues have positive real parts, we have an unstable fixed point or repeller • if all eigenvalues have negative real parts, we have a stable fixed point or attractor • if at least two eigenvalues have real parts with opposite signs, the corresponding fixed point is called a saddle point • if an eigenvalue is zero and at least one other eigenvalue has positive real parts, we have an unstable point • if an eigenvalue is zero and all other eigenvalues have negative real parts, linear stability theory does not suffice.
For more details, see for example [59].

Figure 1 :
Figure 1: Existence and stability regions in (k, λ)-plane.The plotted curves follow from the stability criteria shown in Table2.

Figure 2 :
Figure 2: Phase space with k = 1 and λ = 2. B is a stable node, that is, the only attractor describing a scaling solution with w = w = 0.No acceleration region present.

Figure 3 :Figure 4 :
Figure3: Phase space with k = 1/2 and λ = 3/2.Here the only attractor is Point C where the universe expands as if it is completely matter dominated (scaling solution), while Point B is a saddle point.The shaded region represents the area of the phase space where there is accelerated expansion.

Figure 5 :
Figure 5: Point D + exists in the whole shaded region.Point D − exists only in the dark grey region.

Figure 6 :
Figure 6: The parameter values are λ = 3/2 and k = 8.The eigenvalues of D + are −0.75 ± 1.3713i and the eigenvalues of D − are −2.1570 and 0.6570.The shaded area represents the part of the phase space where there is accelerated expansion.

Figure 7 :
Figure 7: The parameter values are λ = √ 3 and k = √ 8. B is an unstable centre and D + is a stable spiral.The shaded region represents the part where the phase space is accelerating.a stable spiral.It has eigenvalues −3/4 ± i √ 15/4.Similar to the previous case, Point O is an earlytime attractor corresponding to an early-time universe undergoing accelerated expansion.There is no late-time attractor within the acceleration region.

Figure 8 :
Figure 8: The parameter values are λ = 1 and k = 1.The shaded region represents the part where the phase space is accelerating.

Figure 9 :Figure 10 :
Figure 9: The parameter values are λ = 1 and k = 20.The eigenvalues of D + are −0.75 ± 1.11193i and the eigenvalues of D − are −2.34057 and 0.840574.The shaded area represents the part of the phase space where there is accelerated expansion.

Table 2 :
Stability of the critical Points O, A ± and B for system (4.17)-(4.18).The classification assumes that the eigenvalues are non-zero.

Table 3 :
Physical properties of the fixed points for the matter dominated case, for system (4.17)-(4.18).