Interactive mixture of inhomogeneous dark fluids driven by dark energy: a dynamical systems analysis

We examine the evolution of an inhomogeneous mixture of non-relativistic pressureless cold dark matter (CDM), coupled to dark energy (DE) characterised by the equation of state parameter $w<-1/3$, with the interaction term proportional to the DE density. This coupled mixture is the source of a spherically symmetric Lema\^\ itre-Tolman-Bondi (LTB) metric admitting an asymptotic Friedman-Lema\^\ itre-Robertson-Walker (FLRW) background. Einstein's equations reduce to a 5-dimensional autonomous dynamical system involving quasi--local variables related to suitable averages of covariant scalars and their fluctuations. The phase space evolution around the critical points (past/future attractors and five saddles) is examined in detail. For all parameter values and both directions of energy flow (CDM to DE and DE to CDM) the phase space trajectories are compatible with a physically plausible early cosmic times behaviour near the past attractor. This result compares favourably with mixtures with the interaction driven by the CDM density in which conditions for a physically plausible past evolution are more restrictive. Numerical examples are provided describing the evolution of an initial profile that can be associated with idealised structure formation scenarios


Introduction
Current observational evidence supports the existence of an accelerated cosmic expansion, likely driven by an unknown form of matter-energy, generically denoted "dark energy" (DE), and usually described by suitable scalar fields or (phenomenologically) as a fluid with negative pressure [1,2,3]. Observations also point to the existence of cold dark matter (CDM) clustering around galactic halos, usually described in cosmological scales by pressure-less dust, while ordinary visible matter (baryons, electrons and neutrinos) and photons (radiation) comprise less than 5% of the total contents of cosmic mass-energy.
While both dark sources only interact with ordinary matter and radiation through gravitation, it is very reasonable to assume that there is some form of interaction between them. This assumption cannot be ruled out, given our ignorance on the arXiv:1711.09627v1 [gr-qc] 27 Nov 2017 fundamental nature of these sources. In fact, potentially useful information on the primordial physics behind dark sources may emerge by fitting various assumptions of such interactions to observational data, given the fact that interactive DE and CDM is consistent with the dynamics of galaxy clusters [4] and the integrated Sachs-Wolfe effect [5]. Several models of coupled dark sources can also be found in the literature motivated by particle physics, thermodynamics, etc. [1,6].
Observations also suggest that at sufficiently large scales the Universe is well described by linear perturbations of all sources (dark and visible) in an homogeneous Friedman-Lemaître-Robertson-Walker (FLRW) background metric, with non-linear dynamics (whether Newtonian or relativistic) needed to explain the observed local structure [7]. The interplay of local and cosmic dynamics at all scales must comply with the observed anisotropy of the Cosmic Microwave Background (CMB) [1,2,3]. Evidently, this dynamics depends on the assumptions made on DE and CDM, which leads to a model dependent power spectrum that should be contrasted with observations at large scales and in structure formation (from data and from numerical simulations). The observed data should provide interesting constraints on assumptions about the dark sources. In particular, different DE models have been considered in the linear order perturbation scheme in the literature [8,9,10,11].
Conventionally, structure formation scenarios are studied by non-linear Newtonian dynamics (analytically [12,13] and through numerical simulations [14], see review [15]), since CDM is assumed to be practically pressure-less and DE can be modelled (or approximated) by a cosmological constant. However, once we assume a fully dynamical DE source with non-trivial pressure and non-trivial interaction with CDM, it is necessary to utilise General Relativity (whether perturbatively or not) to obtain a valid description of its evolution, since Newtonian gravity can (at best) mimic sources with pressure when adiabatic conditions are assumed (see discussion in [16]). While considering a scalar field is the most common approach to dynamical DE, a phenomenological description by means of fluids with negative pressure can also be useful. Ideally, fully relativistic inhomogeneous DE and CDM interacting sources should be examined through high power numerical relativistic codes (whether assuming a continuous modelling or N-body simulations). However, since the latter codes are in their early development stages [17,18,19,20], we can resort to more idealised (yet still relativistic and non-perturbative) description by means of inhomogeneous exact solutions of Einstein's equations.
In most of the literature the term "Lemaître-Tolman-Bondi (LTB) models" is broadly understood to denote spherically symmetric exact solutions endowed with the LTB metric and associated with a pure dust source [21]. These solutions have been extensively studied (with zero and nonzero cosmological constant) and used in a wide range of astrophysical and cosmological modelling (see extensive reviews in [22,23,24,25]). In particular, a better understanding of their theoretical properties follows by describing their dynamics in terms of "quasi-local scalars" [26,27,28,29] (to be denoted henceforth as "q-scalars"), which are related to averages of standard covariant scalars and satisfy FLRW dynamical equations and scaling laws [28].
Since the deviation from a homogeneous FLRW background can be uniquely determined (in a covariant manner [28]) by fluctuations that relate the q-scalars and the standard covariant scalars, the full dynamics of Einstein's equations is equivalent to the dynamics of the q-scalars and their fluctuations (see discussion in [29]). The q-scalars and their fluctuations allow for a consistent dynamical systems study of the models (with zero cosmological constant in [30,31] and nonzero in [32]). An important theoretical connection with cosmological perturbation theory follows from the fact that the fluctuations of the q-scalars provide an exact analytic (and covariant) generalisation of gauge invariant cosmological perturbations in the isochronous gauge [31,33].
It is less known that LTB metrics admit energy momentum tensors with nonzero pressure in a comoving frame. For a perfect fluid the pressure must be uniform (zero pressure gradients), which allows to interpret the source as a mixture composed by a homogeneous DE fluid interacting with inhomogeneous dust representing CDM (see [34]). For fluids with anisotropic pressure, the latter supports non-trivial pressure gradients, leading to a similar description in terms of q-scalars and their fluctuations as in pure dust LTB models. For anisotropic fluids the anisotropy of the pressure can be related to the fluctuation of the q-scalar associated with the isotropic pressure. Since setting up fluid mixtures is possible and a wide variety of equations of state are admissible, LTB metrics with these sources have been used to model inhomogeneous mixtures of DE and CDM [35,36].
In order to extend earlier work in [35,36] and to explore a generalisation of previous work in [32] that considered DE as a cosmological constant, we studied recently [37] a dynamical systems analysis of an LTB interactive mixture of CDM (dust) and DE (fluid with p/ρ = w, w = const. < −1/3), under the assumption of an interaction driven by CDM: i.e. the interaction term J is proportional (via a dimensionless constant α) to the CDM dust density.
In the present paper we undertake a dynamical systems analysis of a similar (yet qualitatively different) configuration: we assume the same EOS for CDM and DE, but with the interaction now driven by DE, with J now proportional to the DE density. As we show along the paper, the resulting evolution is qualitatively different in both cases. While the phase space of both (CDM or DE driven) mixtures contains as critical points five saddles and past/future attractors. The attractors have very different properties: • The CDM driven mixture examined in [37]: the phase space position of the past attractor (Big Bang) depended on the parameters α and w. For α > 0 (energy transfer from DE to CDM) the past attractor describes a well behaved CDM dominated scenario. However, for α < 0 (energy transfer from CDM to DE) the DE density became negative in phase space regions around the past attractor (irrespective of initial conditions), thus signalling an unphysical past evolution that is inconsistent with all observational data. This problem was already noticed in the literature with this type of CDM-DE mixtures based on FLRW metrics [1,6]. • The DE driven mixture examined here. As opposed to the system in [37], the past attractor is now fixed and the future attractor depends on the choice of parameters α and w. Thus, regardless of the sign of α (directionality of interaction energy flow), the past evolution is now a physically plausible CDM dominated scenario (compatible with observations) while the position of the future attractor describes various plausible DE dominated scenarios that depend on the parameter choices. For α < 0 the future attractor is unphysical (CDM density becomes negative). In this case, the phase space evolution must be appropriately restricted.
Since all observations examine cosmic evolution along our past null cone, the physical plausibility of the models must be determined primarily from their past phase space evolution (the future evolution is much more open to speculation). Hence, a comparison between our results and those of [37] suggests that the DE driven mixture should be favoured, as it exhibits a wider parameter consistency with a past evolution compatible with observations. The plan of the article is summarised as follows. In section 2 we present the qscalar formalism to set up the evolution equations of for an LTB metric whose source is a mixture of CDM and DE with the coupling term proportional to the DE density. The resulting dynamical system and the characteristic features of its phase space are introduced in section 3. In section 4 we classify the critical points for different ranges of the free parameters. In section 5 we solve numerically the dynamical system for initial conditions corresponding to three different choices of the free parameters, leading to three different types of evolution. The associated phase space trajectories and radial profile evolution are displayed in detail. Finally, in section 6 we summarise our findings findings and compare our results with those of [37].

LTB spacetimes, q-scalar variables and coupled dark energy model
Following the methodology described in [35,36] (summarised in [37]), we consider an LTB metric in a comoving frame given by (we choose units with c = 1) where R = R(t, r), R = ∂R/∂r, K = K(r), with the total energy-momentum tensor given by where ρ(t, r) and P (t, r) are the total energy density and isotropic pressure, while the traceless anisotropic pressure tensor is Π a b = P(t, r) × diag[0, −2, 1, 1]. Since we are interested in describing a CDM and DE fluid mixture, we assume the following decomposition of (2) T ab = T ab (m) + T ab (e) ⇒ ρ = ρ (m) + ρ (e) , p = p (m) + p (e) , P = P (m) + P (e) , (3) where the indices "(m)" and "(e)" (whether above or below) respectively denote the CDM and DE mixture components (we adopt this convention henceforth). The total energy-momentum tensor is conserved: ∇ b T ab = 0, but the decomposition above leads to the conservation law for the mixture components where j a is the coupling current that characterises the interaction of both sources. In order to keep the symmetry of the metric, we assume that this current is a vector parallel to the 4-velocity, so that j a = Ju a and h ca j a = 0 hold. The projection along while the spatially projected conservation equation h ac ∇ b T ab = 0 holds for all the evolution.
We have now five state variables A = ρ (m) , ρ (e) , p (m) , p (e) , J which depend on (t, r). As shown in [35,36,37], we associate to each of them a q-scalar A q and a fluctuation δ A by the following rule ‡ ‡ Notice that Aq is related to the proper volume average of A with weight factor √ 1 + K. See comprehensive discussion in [28] where the lower bound of the integrals above x = 0 marks a symmetry centre such that R(t, 0) =Ṙ(t, 0) = 0, withṘ = u a ∇ a R = ∂R/∂t. In particular, it is straightforward to show (see [36]) that Other covariant scalars associated with (1) and (2) are the Hubble expansion scalar H = (1/3)∇ a u a = (R 2 R )˙/(R 2 R ) and the spatial curvature K = (1/6) 3 R = 2(KR) /(R 2 R ), where 3 R is the Ricci scalar of constant t hypersurfaces. Their respective q-scalars are given by while for the interaction term we have J = J q (1 + δ J ), whose dependence on state variables will be determined further ahead. For the mixture (3) to describe CDM and DE we will choose the following equations of state (EOS) (the same as in [37]): DE (barotropic fluid): where we have assumed that w = p (e) /ρ (e) < −1/3 is a constant. Given the EOS's (9) and (10), we will adopt the following convention Notice that (7) and (10) imply that only the DE source contributes to the anisotropic pressure: P = P (e) = δ p (e) /2. As shown in [35,36,37] Einstein's equations reduce to the following system of evolution equations § together with the algebraic constraints (14) § The system (12a)-(12f) can be used to determine the metric coefficients R and R in (1) by supplying the following two evolution equationsṘ = R Hq,Γ = ΓHqδ H with Γ = R /R, which follow from the first equation in (8)  where κ = 8π/3, (14) follows from (13) by applying the second rule of (6), δ k = (K − K q )/K q , while J q is the q-scalar associated (via (6)) to the energy density flux defined from the (or defining a) local J.
Notice that (13) is the quasi-local Hamiltonian constraint, which (from (8)) takes the functional form of an FLRW Friedman equation. Also, the evolution equations (12a-12c) for the q-scalars H q , ρ (e) q , ρ (m) q are formally equivalent to FLRW equations for a CDM-DE mixture with EOS (9) and (10). This reinforces the interpretation of the q-scalars as averaged LTB scalars that mimic at every comoving shell r = r i the corresponding scalars of an FLRW background metric. In fact, as shown in [26], an asymptotic FLRW background follows as all fluctuations δ m , δ e δ H , δ J vanish in the limit r → ∞ for all t, which is equivalent to the fact that in this limit the full system above reduces to the FLRW evolution equations (12a-12c) and the Friedman equation in (13).
The autonomous system (12a-12f) can be solved numerically for a choice of w and J q , which determines δ J through (6). In the present work we will consider an interaction term J q proportional to the DE density and Hubble q-scalars as follows: where α is an dimensionless coupling constant. Note that the interaction energy flows from DE to CDM for α > 0 and from CDM to DE for α < 0.
Interactive mixtures with the EOS (9) and (10) and the interaction energy flux term (15) were considered for an FLRW cosmology in [1,8,10]. This type of FLRW models provide background for a first order gauge invariant perturbation treatment that yields linear evolution equations for the associated perturbations of all sources, including the interaction term J, which can be considered as a phenomenological "black box" or (ideally in principle) related to some (yet unknown) early Universe physics.
As shown in [33], the dynamics of LTB metrics in the q-scalar formalism yields (through evolution equations like (12a-12f)) an exact non-linear generalisation of linear gauge invariant cosmological perturbations in the isochronous gauge (for any source compatible with the LTB metric). The advantage of using numerical solutions of (12a-12f) lies in the possibility to examine in the non-linear regime the connection between the assumptions on the CDM-DE interaction mediated by J and observations on structure formation in the galactic and galactic cluster and supercluster scales.

The dynamical system
The q-scalar formalism described in the previous section allows us to define suitable dimensionless functions that transform the system (12a-12f) into an autonomous fivedimensional dynamical system that is amenable to a qualitative phase space analysis, analogous to that undertaken in [37].
For the density scalar functions we define below the following q-scalars that are analogous to the Ω factors in FLRW cosmologies Interactive mixture of inhomogeneous dark fluids driven by dark energy: a dynamical systems analysis.7 transforming the Hamiltonian constraints in (13)- (14) into the following elegant forms Next, we introduce a dimensionless coordinate ξ(t, r) that will serve as the phase space evolution parameter, so that for all comoving curves r = r i we have In terms of ξ and using the interaction term defined in (15), the system (12a-12f) is transformed into the following dynamical system ∂Ω m (20e) The system (20a-20e) can be solved numerically for a set of initial conditions for every comoving shell r = r i once we fix the free parameters of the model, w and α. We can compute afterwards H q (ξ, r i ) from Once we have computed the phase space variables Ω m q , Ω e q , δ m , δ e , δ H , all relevant quantities can be obtained: The q-scalars associated with the CDM and DE densities and the spatial curvature and its fluctuation δ k follow directly from (16a)-(16b), (17), (18) and (21), while all local quantities follow from Additionally, it is possible to recover physical time from the phase space evolution parameter ξ(t, r) from evaluating at each fixed r = r i , though it is important to bear in mind that hypersurfaces of constant ξ and t do not coincide, thus for every scalar A we have [∂A/∂r] t = [∂A/∂r] ξ (the appropriate integrability conditions are discussed in detail in [32]). Taking this into account, the LTB metric functions follow from evaluating R = exp H q dt and R = R exp − H q δ H dt for H q , δ H as functions of (t, r). Critical points Ω m q , Ω e q , δ m , δ e , δ H Eigenvalues

Homogeneous and inhomogeneous subspaces
As in [37,32], we can split the phase space of (20a-20e) into two interrelated projection subspaces: The homogeneous subspace. It is defined by the phase space variables Ω m q and Ω e q , since they are fully determined by the evolution equations (20a)-(20b), which do not involve the other phase space variables δ m , δ e , δ H . In fact, for every trajectory (fixed r) these evolution equations are formally identical to FLRW equations for the analogous variables. The inhomogeneous subspace. It is defined by the remaining three phase space variables δ m , δ e , δ H , as these provide a measure of the departure of the local scalars from their homogeneous FLRW counterparts.
The study of the phase space will be undertaken by looking at its trajectories in terms of these two projections.

Critical points
The critical points of the system (20a-20e) and their respective eigenvalues are shown in table 1. As expected, they depend on the free parameters w and α, save for P C1.
The critical point P C1 is in fact a line parallel to the δ e axis. The eigenvalue λ 1 of P C1 is zero, corresponding to a eigenvector that is also parallel to the δ e axis, indicating that near the line there is no evolution of the space phase trajectory in that direction. For α < 0, the critical points P C4, P C5, P C6 and P C7 are non physical as their component Ω m q is negative, which means the CDM energy density should be negative. We will examine below the homogeneous subspace closely.

Homogeneous subspace.
The homogeneous subsystem for α > 0 has the following critical points (see figure 1): Interactive mixture of inhomogeneous dark fluids driven by dark energy: a dynamical systems analysis.9 • future attractor: . Both, P CA and P CR can be considered as critical points of the phase space that would result from an FLRW model, or as a projection of the P C1 − P C7 points over the [Ω m q , Ω e q ] subspace in a full five-dimensional representation. In the former case, the trajectories in the phase-space are computed for a given set of initial conditions with δ m = δ e = δ H = 0 and live completely in the homogeneous space, while in the later case the trajectories are computed with a general choice of δ m , δ e and δ H and are represented in the homogeneous subspace as projections of the five-dimensional space trajectories over the [Ω m q , Ω e q ] subspace. Additionally, in a similar way as in the interaction used in [37], we have a one-dimensional invariant subspace (a line) given in this case by where we used (20a-20b). This invariant line contains both the saddle point and the future attractor. Hence, the system can evolve from the saddle point to the future attractor (for initial conditions with Ω m plane is divided in two regions by this invariant line: the region where trajectories evolve from the Ω m q = 0 axis to the future attractor and the region where the trajectories evolve from the past attractor. The later region contains part of the attraction basin of P CA: trajectories that evolve from the past attractor to the future attractor, representing an ever expanding scenario where initially there is only CDE with CDM density increasing from the interaction with DE. This region also contains trajectories for which Ω m q , Ω e q → ∞, which correspond to comoving layers that bounce (since Ω m q , Ω e q diverge as H q → 0). We will not consider the evolution of such trajectories. For α < 0 the future attractor P CA lies in an unphysical phase space region marked by negative Ω m q . For trajectories emerging from the past attractor the physical evolution, which can only be defined up to the invariant line, describes an expanding scenario in which energy density flows from the CDM to the DE component until the CDM density vanishes on the comoving shells (at different times for different shells). However, the fact that the past evolution is not unphysical makes the coupling term (15) acceptable also when α < 0, as has been stated in the literature [38] deling with these CDM-DE mixtures in FLRW cosmologies. This stands in sharp contrast with the coupling used in [37], where α < 0 leads to grossly unphysical past evolution, which implies considering only the coupling with α > 0 (as in FLRW cosmology scenarios).

Initial conditions, scaling laws and singularities.
To specify initial conditions to integrate the dynamical system (20a)-(20e) we need to provide an initial value formulation for the LTB models under consideration. Proceeding as in [37], we specify initial conditions given at an arbitrary hypersurface t = t 0 (subindex 0 will denote henceforth evaluation at t = t 0 ). It is useful to write LTB metric (1) in the following FLRW-like form where L = L(t, r) is analogous to the FLRW scale factor. Since the LTB metric admit an arbitrary rescaling of the radial coordinate, we can always define a convenient radial coordinate by specific choices of R 0 (r). We can identify L = 0 as the locus of the Big Bang singularity, while Γ = 0 marks the locus of a shell crossing singularity [32]. From (8), (12a), (12b) and (12c) we can see how the q-scalars scale as their equivalent FLRW scalars, that is which lead to Initial conditions to integrate the system (20a-20e) follow from specifying initial profiles ρ  (6) with R = R 0 . These initial profiles must comply with the condition δ m 0 , δ e 0 , δ H 0 → 0 as R 0 (r) → ∞ to guarantee the existence of an asymptotic FLRW background (see [26]).
Whether comoving shells expand or bounce/recollapse can be determined from (27) by looking at the roots of H q throughL 2 = L −1 Q(L) where with: a = Ω m q0 + α where we used the fact that dξ = H q dt leads to ξ = ln(L). For each choice of initial conditions Ω m q0 , Ω e q0 and free parameters w and α, if Q(L) has real roots for a given comoving shell, the latter will bounce at the value of ξ where Q = 0. Conversely, if Q(L) has no real roots the layer has an ever expanding evolution. For the remaining of the paper we will only consider the phase space evolution of expanding layers.
To obtain analytic solutions we need to solve (27) and obtain Γ from the relation L /L = R 0 (Γ − 1)/R 0 , where the radial derivatives must be evaluated for constant t (see [37])). The scaling laws for the fluctuations can be found from (26). For example, using the definition of δ e , it is straightforward to show that which can (in principle) be used to evaluate Γ once we have a solution t = t(L, r) of (27). Since analytic solutions of (27) may only exist for very restricted values of α, w, scaling laws like (30) are not useful. In general, the evolution of the models needs to be determined numerically.
If Γ = 0 we have a shell crossing singularity, which means that initial conditions should be found to avoid this happening (it is not possible to provide simple guidelines for this, as in dust solutions with zero cosmological constant, see [22,23,24,27]). As shown in previous work (for example [37]) the q-scalar formalism fails at shell crossing singularities because the fluctuations δ m , δ e , δ H diverge. Hence, we will select initial conditions such that shell crossings are avoided: Γ > 0 holds throughout the full phase space evolution.

Critical points in terms of the parameters w and α
In this section the critical points of the system are studied for different possibilities of the free parameters w and α. As both parameters are widely used in the FLRW model, we will consider a parameter range that is common in Cosmology.
The constant EOS parameter w plays a similar role as in analogous CDM-DE mixtures based on FLRW models. While observational data seems to favor the Λ-CDM model for which w = −1 holds exactly, small variations from this value are still possible [3]. We will henceforth adopt the current terminology, in the literature by referring to DE with w > −1 as "quintessence models" and w < −1 as "phantom models". The latter models present several theoretical problems, such as the violation of the second law of thermodynamics once we assign an entropy to the phantom fluid, or the presence of a negative kinetic energy of the phantom field term (when described by a scalar field) [1]. In this article we will consider quintessence and phantom models with w close to −1.
Considering the same interaction in the present article, the CDM-DE mixture in FLRW geometry examined in [40] finds that the second law of thermodynamics (based on the entropy of DE as an effective field) is violated if α < 0, while the entropy is zero for a scalar field in a pure quantum state. They also find that for positive α it is necessary to impose the bound α < 0.1 in order to reproduce the observed values of BAO and CMB anisotropy [8,9]. On the other hand, in [41] the evolution of linear perturbations in a FLRW background sharing our assumptions on CDM, DE and J q leads to the bounds −0.22 > 3α > −0.90 to comply with the constraints of CMB anisotropy. These results are specially interesting, since the dynamics of LTB solutions described by q-scalars and their fluctuations can be mapped to linear perturbation on an FLRW background [33]. While the spherical symmetry of LTB models allows for the description of a single structure, the latter can be studied exactly in full non-linear regime. In the present paper we will assume positive and negative values of α.
The critical points P C4 − P C7 share the values Ω m q = −α/w and Ω e q = 1 − α/w of the homogeneous projection, but are distinct in the inhomogeneous projection coordinates. The points P C5 and P C7 are always saddle points in the range of free parameters considered. For a choice of parameters such that 1 + w + α > 0, the critical point P C4 is a future attractor as all their eigenvalues are negative defined, while P C6 is saddle point. On the other hand, when 1 + w + α = 0, P C6 has the same components as P C4, and behaves as a non hyperbolic point as one of its eigenvalues are null while the rest are negative defined. Finally, for 1 + w + α < 0, P C6 is the future attractor while P C4 is a saddle point. Figure 1 shows the homogeneous subspace together with the critical points P CR, P CA, P CS and the invariant line for both cases: α > 0 in panel (a) (in this case we have chosen α = 0.1 and w = −1), and α < 0 in panel (b)(α = −0.1 and w = −1). Some numerically computed trajectories are shown for illustration purposes only. We have chosen to represent arctan(Ω m q ) vs. arctan Ω e q in order to deal with finite values in the plots. The same criterion is used for the rest of the homogeneous projection plots. Figure 2 shows the two inhomogeneous projections of the system (20a-20e) and some numerically computed trajectories. Panel (a), (b) and (c) represent the projection with Ω m q = −α/w and Ω e q = 1 − α/w and the critical points P C4, P C6 and P C7 for different choices of w and α > 0: panel (a) shows a set of parameters where 1 + w + α > 0 and the future attractor is P C4; panel (b) shows a set with 1 + w + α = 0 where P C4 and P C6 are the same point; and, finally, panel (c) shows a set of parameters is represented where 1 + w + α < 0 and the future attractor is P C6. The critical point P C5 is not represented as it is always a saddle point with a large δ m component, far away from the rest of points. Panel (d) shows the projection Ω m q = 1 and Ω e q = 0 for α = 0.1 and w = −1, although choosing a different value of w and α will not change the general behavior of the points or the trajectories. The projection with Ω m q = −α/w and Ω e q = 1 − α/w will be unphysical when α < 0. The projection Ω m q = 1 and Ω e q = 0, on the other hand, will be physically plausible and phenomenologically identical to that in panel 2d in the α > 0 case.

Energy density flow from DE to CDM (α > 0)
In this case all seven critical points are physical. The attractor is a different point for the different choices of w and α, as stated above. The presence of a future attractor for α > 0 allow us to find initial profiles that present inhomogeneities in the attraction basin of it, i.e the fluctuations (δ functions) evolving to constant values given by the components of the corresponding critical point. This behaviour is examined further ahead.

Quintessence and cosmological constant cases
When w ≥ −1 and α > 0, the critical point P C4 acts as a future attractor and the trajectories nearby evolve to it. On the other hand, the critical point P C2 is a past attractor as all the eigenvalues of the system computed near P C2 have positive values. The rest of the critical points are saddle points with their own attraction subspace generated by the corresponding eigenvectors.
Panel 1a shows the homogeneous subspace for w = −1: it is formally identical to that of the w > −1 case, except for the position of the future attractor and the shape of the invariant line. In the panel 2a, the inhomogeneous projection Ω m q = −α/w and Ω e q = 1−α/w is shown for the w > −1 case. The attractor P C4 is also shown together with some trajectories in its vicinity that evolve to it. Also, the saddle points P C6 and P C7 are displayed (the point P C5 is not shown given that it is located far away and is always a saddle point). Finally, panel 2d displays the inhomogeneous subspace with Ω m q = 1, Ω e q = 0, plotted for the w = −1 case but, again, it is phenomenologically identical to the w > −1 case. In this projection the past attractor P C2 and the saddle points P C1, P C3 are also displayed.

Phantom dark energy case
The points P C4 and P C6 can be (respectively) a future attractor and a saddle point when 1 + w + α > 0 or (respectively) a saddle point and a future attractor when 1 + w + α < 0. Finally, when 1 + w + α = 0 both points have the same coordinates and the point is non-hyperbolic, i.e. it is an attractor in some directions and there is no evolution near it in other directions. The rest of points behave in a similar way as in the previous cases for any choice of free parameters.
The homogeneous subspace is identical to that of panel 1a except for the position of the P CR point and the slope of the invariant line. In the panel 2c, the inhomogeneous projection Ω m q = −α/w and Ω e q = 1 − α/w is displayed for a choice satisfying 1 + w + α < 0. The point P C4 is a saddle point and the point P C6 is now the attractor of the system, in contrast with the case 1 + w + α < 0 that will be similar to the quintessence and cosmological constant cases. Finally, the inhomogeneous subspace with Ω m q = 1, Ω e q = 0 is as in the previous case similar to that in 2d.

Energy density flow from CDM to DE (α < 0).
When α < 0, the energy flows from the CDM to DE. In this case only P C1, P C2 and P C3 have physical meaning while the rest of the points represent values with Ω m q < 0. In the homogeneous subsystem, the future attractor P CA is no longer physical and consequently the invariant line is also non physical. The trajectories in the homogeneous subsystem evolve from the past attractor to the Ω m q = 0 axis or to infinity. When a given shell reaches the Ω m q = 0 axis, we can assume that the CDM content of that shell has been consumed by the coupling term. Such points occur at different values of ξ (and thus different cosmic times). The evolution of the mixture is only physically meaningful up to these points.
There is no significative difference between the homogeneous space for the different possibilities of the parameter w. Although the trajectories follow a different curve for every choice of w and α, they all evolve form the critical point P CR. Panel 1b shows schematically the homogeneous subspace for α < 0. The behavior of P C1, P C2 and P C3 in the inhomogeneous subspace with Ω m q = 1, Ω e q = 0 is identical to the α > 0 case, plotted in the panel d of figure 2.
The lack of future attractor for the δ A functions (A = m, e, H, phase space variables of the inhomogeneous subspace) in this case makes it possible for some of their initial profiles to make them evolve: to infinity (shell crossing), or to some q , H to be positive). Consequently, it is not possible to keep an evolution with a δ A for more negative values than the limit −1, as this would imply negative local densities. We can argue that the evolution equations yield unphysical conditions if (somehow) δ m , δ e < −1 holds. In the next section, we explore this problem for a specific initial profile.

Numerical example of idealised structure formation scenarios.
In this section we examine various numerical examples of evolving radial profiles that could lead to potentially interesting structure formation scenarios. Our intention is not to model any type of realistic configuration, but to illustrate the dynamical evolution of the mixtures through simple idealised examples. We choose an appropriate form for R 0 (r) defined for an interval r ∈ [0, r max ] and examine the evolution equations for fixed values of r specified by a partition of n elements in this interval.
To look at the numerical evolution of the initial profiles we define the dimensionless time parameter (different from ξ) given byt = H s t where H s is an arbitrary constant with time inverse dimensions (in cosmological applications it is customary to choose H s = H 0 ). This rescaling of time introduces a rescaling of the remaining variables:H q = H q /H s , κρ (m) /3 = κρ (m) /(3H 2 s ), κρ (e) /3 = κρ (e) /(3H 2 s ). For simplicity we will drop the bars on the normalised variables and will set the arbitrary scale as H s = 1, which fixes the energy density normalisation scale as κ/(3H 2 s ) = 1 (see [37]). In order for the initial profiles to define a structure formation scenario we need some of the "inner" shells (values of r around the symmetry centre r = 0) that initially expand, but at some t bounce and collapse (H q changes sign from positive to negative), whereas "outer" shells continue expanding (H q > 0 holds for all t). The bounce is defined by H q = 0, hence we can define for each r a value t = t max (r) such that H q (t max (r), r) = 0. Notice that the dynamical systems study we have undertaken does not examine phase space trajectories of shells that have bounced and then collapse (H q → 0 evolving towards H q < 0), as both coordinates [Ω m q , Ω e q ] of the homogeneous projection diverge as H q → 0. The numerical study given in this section will compensate for this deficiency.
For the mixed expanding/collapsing type of evolution described above we need the following homogeneous subspace trajectories: (i) the outer ever expanding shells must evolve from the past to the future attractor (or to the Ω e q axis when α < 0); (ii) inner shells must evolve from the past attractor to infinity Ω e q , Ω m q → ∞ as t → t max . For a bounce/collapse regime (and pending on specific initial conditions), the variables of the inhomogeneous subspace could also diverge or not evolve to the future attractor.
As mentioned before, if δ A = −1 on a given shell and A q = 0, then A = 0, which for positive definite quantities (densities) implies that the LTB dynamical yield an unphysical evolution for decreasing δ A < −1. This problem tends to occur specially in the cases with α < 0, when no physical attractor is present, but it may also occur for some configurations in the α > 0 cases, where the inhomogeneous attractor is physical but the initial profiles were set with the initial δ A functions out of its attraction basin.
We have chosen the following initial profiles to be used to probe the models for three different sets of free parameters α, w K 0 = k 10 + k 11 − k 10 1 + tan 2 (r) , k 10 = −4.10, k 11 = 7.50; together with the coordinate choice R 0 (r) = tan(r). Hence, we consider a partition of n = 20 elements for r going from 0 to π/2.  Figure 3a shows the homogeneous projection of the configuration with the invariant line of the system in grey, and the initial conditions for every shell as red points. Figure 3b shows the inhomogeneous projection of the trajectories, and the initial conditions as red points. Figure 4 displays the radial profiles of the local scalars A = H, ρ (m) , ρ (e) and J at different instants of time in the plane arctan(A) vs. arctan(R 0 (r)). For t = 0.50 no shell has collapsed yet, for t = 0.70 shells 1, 2 have collapsed, for t = 1.00 the shells 3, 4 have collapsed, for t = 2.00 shell 5 has collapsed and the last inner shell is about to collapse. Finally, for t = 4.00 all inner shells have already collapsed and outer shells evolve into a homogeneous profile. For expanding trajectories the q-scalars and the fluctuations δ A for the outer shells tend to their attractor values, while the profiles of local scalars A tend to a constant profile, as δ A → 0 (A = m, e, H) for the attractor P C4. It takes a long time for the functions to evolve into their attractor values. The evolution of the these profiles to a constant profile can be appreciated in panel 4a for the local scalar H at t = 4.00, 6.00, 9.00 and in the panel 4b for the local scalar ρ (m) at the instants t = 4.00, 6.00, 9.00. The local scalar ρ (e) , which was initially constant, needs an even longer time to evolve into a constant profile, but the line representing it at t = 9.00 is clearly more homogeneous at the outer shells than in previous instants.

5.2.
Positive α and 1 + w + α < 0 We choose α = 0.1 and w = −1.15. Only the shells r 1−3 collapse, while the rest evolve towards the attractor P C6. As for the other choice of parameters, the initial δ A functions are in the attraction basin of P C6. In this case, the inhomogeneous projection of the attractor is not zero as δ m = δ e = 0.05, δ H = 0.05/2. Consequently, from (6), the profiles of ρ (m) and ρ (e) do not evolve to a constant profile as in the previous case, but to a profile whose r-dependence is given by ρ (m) = ρ (e) = R 0.15 (1.05), while the local H is given by H = R 0.07 (1 + 0.05/2). Figure 5a shows the homogeneous projection of the configuration with the invariant line of the system in grey, and the initial conditions for every shell as red points. Figure 5b shows the inhomogeneous projection of the trajectories evolving to the future attractor, and the initial conditions as red points. Figure 6 displays profiles of local scalars at different instants of time. At the instant t = 0.80 no inner shell has collapsed yet, at the instant t = 1.00 shells i = 1, 2 have already collapsed and at the instant t = 2.00 all the inner shells have collapsed. Since the q-scalars and fluctuations δ A for the outer shells tend to their attractor values while trajectories expand, the profiles of local scalars tend to the terminal profile shown in panel 6a for H at t = 4.00, 6.00 and in the panel 6c for ρ (e) at instants t = 4.00, 6.00.

5.2.1.
Negative α For any choice of α < 0, the future attractor is no longer physical.
To illustrate this case, we chose α = −0.1 and w = −1. The shells r 1−4 collapse while outer shells, r 5−20 , keep their expanding evolution up to a point where δ m = −1 and the LTB evolution is no longer physical. In particular for this profile and this choice of parameters the function δ m tends to −1 very rapidly for the outer shells. Figure 7a shows the homogeneous projection of the configuration with the invariant line of the system in grey, and initial conditions for every shell as red points. All the outer shells evolve to the Ω m q = 0 axis. Figure 7b displays the plot δ m vs. ξ for the outer ever expanding shells i = 5 − 19. The first shell to reach the value δ m = −1 is i = 20, which is not represented in the figure as δ m → −1 occurs immediately for this shell. The value ξ i for which δ m = −1 occurs (i.e. δ m (r = r i , ξ i ) = −1) puts an upper limit to the range of ξ for which we can use eqs. (20a-20e) to obtain all the local scalars. Figure 8 shows local profiles of various scalars at different instants of time. The profile of the local scalar ρ (m) tends to zero for every shell at the instant mentioned

Conclusions.
We have undertaken a full study of the phase space evolution of expanding and interactive CDM-DE mixtures, under the assumption that the interactive term (see (15)) is proportional to the DE energy density. These mixtures are the source of an inhomogeneous and spherically symmetric exact solution of Einstein's equations characterised by an LTB metric. The present article, together with a recent article [37], generalise previous work [32] for a CDM source with DE modelled by a cosmological constant. As in [32,37], we examined the dynamics of these LTB solutions by means of q-scalars and their fluctuations, which transform in a natural way Einstein's equations into a dynamical system evolving in a 5-dimensional phase space. The latter was studied in terms of two interrelated subspace projections: the 2-dimensional homogeneous subspace whose variables are q-scalars (Ω m q , Ω e q ) analogous to covariant scalars of an FLRW model with same type of CDM-DE mixture, and a 3-dimensional subspace involving the fluctuations of the q-scalars (δ m , δ e , δ H ) that control the inhomogeneity of the models (deviation from FLRW).
The critical points associated with the phase space are listed in Table 1: a past attractor, a future attractor and 5 saddle points. All of them (save the past attractor) depend on the two constant free parameters of the solutions: the EOS parameter w and the proportionality between the interaction term and DE density, α, whose sign determines the directionality of the interaction energy transfer (CDM to DE for α > 0 and DE to CDM for α < 0). The phase space evolution was examined for "quintessence" models (−1 < w < −1/3) and "phantom" models (w < −1), keeping in either case w close to the value -1 that is favoured by observations.
It is important to compare our results with those found in our recent study in [37] involving a similar CDM-DE mixture, but with the interaction term proportional to the CDM density through the dimension-less constant α. The main difference between this assumption and that of the present work (interaction proportional to DE density) is in the parameter dependence of the past attractor, which in both mixtures can be associated with the initial Big Bang singularity. For the mixture of [37] the phase space position of this attractor depended on α, w, while in the present mixture it does not, which means that it is a fixed point in the phase space.
The above mentioned difference in the phase space position of the past attractor has very important consequences, as all available cosmological observations survey our past light cone. For α < 0, and for every set of initial conditions, the past attractor in [37] was located in an unphysical phase space region (DE density becomes negative), hence α > 0 (energy flows from CDM to DE) was the only physically plausible choice that can be (in principle) compatible with observational constraints for all choices of EOS parameter w. As a contrast, in the mixture examined here we found that regardless of the sign of α the past attractor is fixed, taking physically plausible values: an Einstein de Sitter state with zero DE density and positive CDM density with unity Omega factor (Ω e q = 0, Ω m q = 1). Hence, both directions of the interaction energy transfer are (in principle) compatible with observations.
Since the future attractors correspond to times much beyond the present cosmic time, they cannot be contrasted with observational data, and thus are more amenable to speculation. For the case examined in [37] this attractor simply marked a fixed de Sitter state (Ω m q = 0, Ω e q = 1). However, in the present study the phase space position of the future attractor depends on the choice of w, α, and for α < 0 this position is not in a physically meaningful phase space region (CDM density becomes negative). While this can be problematic (as all trajectories must terminate in this attractor), it can still be acceptable provided we only consider the evolution of the models up to phase space points where the CDM density vanishes. Such points correspond to values for which δ m = −1 that mark different cosmic times for different comoving shells. On the other hand, for α > 0, the shells with initial conditions in the attraction basin evolve to a point where both CDM and DE reach a terminal energy density. In this case a choice of parameters with 1 + w + α > 0 lead to a final homogeneous state, as the phase space variables of the inhomogeneous subspace (the fluctuations) tend to zero. For 1 + w + α < 0 the comoving shells evolve towards a nonzero density profile, since the fluctuations tend to nonzero values that depend on the radial coordinate.
The comparison between the results of the present work and those of [37] provide what can be, perhaps, the most interesting conclusion that follows from the present article: if CDM and DE are assumed to interact, then mixtures in which the interaction term is proportional to the DE density (as in this paper) offer more possibilities for educated speculation, and thus could be preferable over those in which it is proportional to the CDM density (as in [37]). The reason for this is, as explained before, that the past attractor (which determines the past evolution surveyed by observations) is always physical for any reasonable choice of parameters, in contrast with the coupling used in [37] where the choice of α < 0 was not possible.
We also examined three structure formation scenarios, given (for example) by the profile in (32). These initial conditions correspond to inner expanding shells that are not in the attraction basin while the rest of the expanding shells (outer shells) that evolve to the attractor or to the Ω m q = 0 axis. The inner shells will evolve to a point where they bounce and start a collapse. The outer shells expand forever and evolve into a profile determined by the choice of 1 + w + α, as mentioned before. Some examples of those scenarios can be found in section 5, where we have computed the local scalar functions over physical time t.
It is also important to compare our results (and those of [37]) with those obtained for DE models or similar CDM-DE mixtures in FLRW cosmologies [1,6,8,9,10,11,40,41], since the dynamics of LTB solutions described by q-scalars and their fluctuations can be mapped to linear perturbation on an FLRW background [33]. While the spherical symmetry of the LTB models allows for the description of a single structure, the evolution of the latter can be studied exactly throughout the full non-linear regime.
We believe that our results can provide interesting clues to test theoretical assumptions on dark sources (DE and CDM) in terms of observations in the scales of structure formation and in the non-linear regime. In fact, it is straightforward to generalise LTB models to non-spherical Szekeres models, which are endowed with more degrees of freedom and thus allow for a fully relativistic description and modelling of multiple structures [42,43,44]. We are currently undertaking further efforts to extend the present work to probe non-linear observational effects of theoretical assumptions on CDM and DE sources. In particular, we aim at considering "spherical collapse models", as well as less idealised non-spherical models, whose source is the type of CDM-DE mixtures we have examined here, but now attempting to fit more realistic observational constraints of structure formation.