Mimetic Tensor-Vector-Scalar Cosmology: Incorporating Dark Matter, Dark Energy and Stiff Matter

Phenomenological implications of the Mimetic Tensor-Vector-Scalar theory (MiTeVeS) are studied. The theory is an extension of the vector field model of mimetic dark matter, where a scalar field is also incorporated, and it is known to be free from ghost instability. In the absence of interactions between the scalar field and the vector field, the obtained cosmological solution corresponds to the General theory of Relativity (GR) with a minimally-coupled scalar field. However, including an interaction term between the scalar field and the vector field yields interesting dynamics. There is a shift symmetry for the scalar field with a flat potential, and the conserved Noether current, which is associated with the symmetry, behaves as a dark matter component. Consequently, the solution contains a cosmological constant, dark matter and a stiff matter fluid. Breaking the shift symmetry with a non-flat potential gives a natural interaction between dark energy and dark matter.


I. INTRODUCTION
The Universe is found from many experiments and observations to include 25% of matter and 68% of dark energy [1][2][3][4].There have been many models for dark matter [5][6][7][8][9].The description of dark energy, however, suffers from theoretical problems [10][11][12][13]), what gives motivation for modified theories of gravity [14][15][16], as the General Relativity (GR) could not describe the data through the whole range of scales.One interesting class of modified theories of gravity are theories, where a scalar field mimics the dark matter [17] using a Lagrange multiplier [18,19] (see also [20] for a review and [21][22][23][24] for some developments of the mimetic theory).In its original formulation, the mimetic theory of gravity can be obtained starting from the GR with a Lagrange multiplier that forces the kinetic part of the scalar field to behave as a constant in time [25].Consequently, the contribution of the mimetic field in the gravitational field equations is a pressureless perfect fluid.On a Friedmann-Robertson-Walker (FRWM) metric this fluid behaves precisely as a dust component on the cosmological scales.A minimal generalization of the mimetic dark matter was proposed by including a potential for the mimetic scalar field [19].The cosmological dynamics of the theory was further studied in [26].The theory contains dynamic dark energy along with mimetic dark matter, and it was stressed that the early-time solution is not dominated by stiff matter unlike in simpler scalar-field models of dark energy, like in quintessence.
The Modified Newtonian Dynamics (MOND) [27,28] is also an alternative successful explanation for the flat rotation curves of galaxies.By violating the Newton's second law at low accelerations, the Tully-Fisher relation is recovered, without introducing an additional dark matter [29][30][31].The Tensor-Vector-Scalar theory (TeVeS) is a covariant theory, which produces MOND in the low energy limit [32].TeVeS is a stable theory and free from ghosts [33].A new relativistic theory for MOND has been proposed [34] in order to reproduce the key cosmological observations, while still reproducing all the galactic and lensing phenomenology in the same way as TeVeS does.
Cosmology with vector fields has been considered widely in the literature [35][36][37][38][39][40].General consideration of scalar-vector-tensor theories and their cosmological solutions is given in [41].A vector field model of mimetic matter was proposed in [25].The consistency of those theories with respect to the absence of ghost instability was studied in [42], including a proposal of a new Tensor-Vector-Scalar gravity (MiTeVeS), which generalizes the previous models.The vector field was introduced to enable the possibility of rotating flows of mimetic matter and to avoid caustics of the geodesic flow [25].The scalar field in MiTeVeS, with its coupling to the vector field, has a significant impact on cosmology.In this paper we discover the homogeneous solution for this theory considering couple of possible interactions for the vector and scalar fields.
First in section III we confirm that the theory with the added interactions does not suffer from pathologies like ghosts, gradient instabilities or tachyons.The homogeneous solution is obtained in sections IV and V. Cosmological constraints on the solution are obtained in section VI.Interacting dark energy and dark matter are discussed in section VII.

A. Lagrangian
Gravitational theories with mimetic matter involve two metrics which are related by a conformal transformation.The physical metric Φ 2 g µν is coupled minimally to all matter fields, where Φ is the conformal factor field that relates the physical metric to the metric g µν .The Mimetic TeVeS theory includes a scalar field φ and a vector field u µ .The conformal factor is related to the scalar and vector fields by including a constraint with a Lagrange multiplier λ.The Lagrangian is written as where R is the Ricci scalar determined by the metric g µν , the kinetic term of the vector field is defined in terms of and the Lagrangian density L m,int (Φ 2 g µν , χ, ∂χ, u µ , φ, φ ,µ ) represents matter fields (denoted collectively as χ) and all interactions among the gravitational fields and the matter fields.Matter couples universally to the physical metric Φ 2 g µν .
The concept of scale invariance is an attractive possibility for a fundamental symmetry of nature.Dimensionless coupling constants for example appear to be related to good renormalizability properties.In its most naive realizations, scale invariance is not a viable symmetry, since nature seems to have chosen some typical scales.The scale invariance can nevertheless be incorporated into realistic generally covariant field theories.Scale invariance has to be discussed in a more general framework than that of standard generally relativistic theories.In the MiTeVeS theory, the action incorporates the scale invariance under the transformation where L = L( x, t) is an arbitrary function of space and time.The fields φ and u µ are the invariant of the theory: This indicates that the theory is capable to have inflationary solutions naturally.We can simplify the derivation of the field equations by fixing the conformal symmetry of the theory.We shall choose the gauge condition as Φ = 1, which gives the theory in the Einstein frame form [42].This is equivalent to performing the following conformal transformation in (1): In this frame the Lagrangian is obtained without the auxiliary field Φ as where g µν is now the physical metric that couples to matter universally.This frame is equivalent for the setting where G is the Newtonian Constant.We assume in the following analyses that the gravitational vector and scalar fields do not couple directly to matter fields: Furthermore, we consider the two simplest interactions for φ and u µ , which involve u µ up to the second power and are at most linear in φ ,µ .Terms that involve higher powers of φ ,µ would modify the kinetic term of the scalar field φ, which will not be considered in this work.These assumptions on L m,int are made here mainly for simplicity.More generally, interactions of higher order as well as couplings between matter and the scalar and/or vector fields are possible.We also note that quantum corrections may induce such terms, when there is no symmetry to prevent it.

B. Field equations
Varying the Lagrangian (6) minimally coupled to the interaction term (9) with respect to g µν , λ, φ and u µ gives the field equations as where in the modified Einstein equation (10a) we have defined the energy-momentum tensors as: where the parentheses around indices denote symmetrization: which ensures the symmetric form of the energy momentum tensor.The contributions of the interactions between the scalar field and the vector field to the energymomentum tensor are given in ( 14) and (15).

III. HEALTHINESS OF THE THEORY
In this section we discuss the absence of ghosts, propagation of perturbations and stability, particularly regarding the interactions that we have introduced in (9).First we complete the Hamiltonian analysis [42] with the interactions (9).The structure of the Hamiltonian analysis remains largely unchanged, since the interactions (9) involve only one derivative of the scalar field.Hence we only present the significant changes to the Hamiltonian analysis achieved in [42].We have performed the Hamiltonian analysis of MiTeVeS with the usual method suited for constrained systems [43][44][45].
In the last subsection we discuss the relation of MiTeVeS and the Einstein-aether theory of gravity, particularly regarding the requirements of stability.

A. Hamiltonian analysis
We assume foliation of spacetime to a union of spacelike hypersurfaces and the Arnowitt-Deser-Misner parameterization of the metric, The variables of the Hamiltonian formulation are the lapse N , the shift N i , the induced metric h ij on the spatial hypersurface, the scalar fields λ, Φ and φ, and the components of the vector field u i and u n that are tangent and normal to the spatial hypersurfaces, respectively.The canonical momenta conjugate to those variables are π N , π i , π ij , p λ , p Φ , p φ , p i and p n , respectively.
The covariant derivative determined by the metric h ij is denoted with D, and R is the corresponding Ricci scalar.
We obtain the Hamiltonian as where In (17) all the fields denoted with v are Lagrange multipliers.We have denoted the De Witt metric in (18) as The primary constraints are The preservation of the primary constraints π N and π i during the time evolution of the system implies following secondary constraints The requirement of the preservation of the constraint p λ ≈ 0 implies a secondary constraint, The preservation of the constraint p n gives a secondary constraint, The constraint D is the first class constraint that generates the scaling transformation.H T , π N , H i and π i are the first class constraints which reflect of the diffeomorphism invariance of the theory.The four remaining constraints p λ , C n , p n and C λ are the second class constraints, which can be interpreted as follows.The constraints p λ and C n can be used to fix the pair of canonical variables λ and p λ , while the constraints p n and C λ can be used to fix the variables u n and p n .Overall the Hamiltonian structure remains similar to [42].
There is no sign of presence of ghosts in the Hamiltonian.

B. Propagation of perturbations
We consider the linearization of the theory.The tensor, vector and scalar fields are written as where g µ and φ (0) are the background and ḡµν , ūµ and φ are small amplitude perturbations with a length scale of variations much shorter than the length scale of variation of the background.We study the propagation of perturbations For simplicity we consider the Minkowski background: , 0, 0, 0).This background can be obtained in the limit of vanishing mimetic energy-momentum in the background: We consider the field equations (10) to lowest order in the perturbations.The Einstein field equations (10a) are linearized as Ḡµν = T tot µν , where Ḡµν is the linearized Einstein tensor and T tot µν is the total energy-momentum, where the scalar and vector fields appear like additional matter components.The propagation of the metric perturbations is similar to GR.The rest of the field equations give the linearized constraint and the linearized field equations for the scalar and vector fields where we have defined For the constant functions (52) considered in Sec.V we would have B = 0.The constraint (26) determines one of the three involved perturbations in terms of the other two.Differentiating the vector equation (27b) we obtain 2 φ = 0, which according to the scalar equation (27a) implies h 2 (φ 0 )∂ µ ūµ = 2B φ; in the case of the constant functions (52) the perturbation of the vector field is transverse.Hence the scalar and vector equations ( 27) simplify as where we have defined The functions V , f , h 1 and h 2 are chosen so that C ≥ 0. The perturbations of the tensor and scalar fields satisfy standard wave equations and propagate with the speed of light (c = 1).The perturbation ūµ of the vector field satisfies the forced wave equation (29b), where the forcing term is the gradient of the scalar perturbation.The solution of (29b) is a sum of the solution of the homogeneous wave equation 2ū µ = 0 and the solution of the inhomogeneous wave equation, ūµ = ūH µ + ūI µ .The latter can be written with a Green's function as where the Green's function G satisfies 2G(x) = δ(x) with suitable boundary conditions.Consider the plane wave solutions to the homogeneous wave equations φ = e ikµx µ and ūH µ = µ e ikµx µ .Then the solution of (29b) takes the form where we have used that the Fourier transformation Ĝ(k) of the Green's function satisfies k 2 Ĝ(k) = −1.This solution tells us three significant things.The solution (32) is a plane wave with a k-dependent amplitude and it gives the dispersion relation k 2 = 0 by substitution into (29b).The k-dependent part of the amplitude of ( 32) is parallel to k µ = (ω, k), which means that if we write ūµ = (ū 0 , ū⊥ + ū ), where ū⊥ is perpendicular to k and ū is parallel to k, the perpendicular component ū⊥ = ⊥ e ikµx µ , i.e., only the components ū0 and ū have k-dependent amplitudes.Thirdly, the scale of variation of the amplitude is given by the parameter C/µ 2 .The analysis of perturbations on a general curved background is somewhat longer.The key point to observe about the field equations (10) is that the kinetic matrix of the Lagrangian is diagonal with respect to the tensor, vector and scalar fields, and furthermore it is a direct sum of the Einstein, Klein-Gordon and Maxwell kinetic terms, respectively.Therefore, as long as the functions V , f , h 1 and h 2 are chosen reasonably (as demonstrated above), there are no instabilities or tachyons.
Finally, we remark that the interaction terms we have introduced in (9) represent attractive interactions between the vector and scalar fields.In the linearized theory, the interaction terms − 1 2 h 1 ūµ ḡµν ūν − 1 2 h 2 ūµ ḡµν φ,ν describe an interaction of the vector field with itself and an interaction between the vector and scalar fields.Both interactions are mediated by the tensor field, i.e., the graviton, which is a massless spin-2 field.It is well known that an interaction mediated by an even integer spin field is always attractive.Such an attractive interaction cannot drive the system to an illimitably excited state.So it is clear that the addition of these interactions could not have introduced instabilities.

C. Relation to Einstein-aether theory and potential issues with stability
When the function f (φ) is a constant, the vector field u µ is similar to the aether vector field of Einstein-aether theory [46].In Einstein-aether theory the vector field is constrained to be a timelike unit vector, u µ u µ = −1.Einstein-aether theory has been studied extensively as a way to consider the consequences of dynamical breaking of local Lorentz invariance such as variation of the speed of light or high frequency dispersion.In our theory, the vector field is coupled to the scalar field φ in the Lagrangian (6).In this work we only consider the two interactions in (9) in addition to the constaint (10b).
The stability of the aether vector field on Minkowski background spacetime was studied in [47].It was found that only three kinetic terms have no linear instabilities or negative-energy ghosts: the sigma model ∂ µ u ν ∂ µ u ν , the Maxwell term F µν F µν , and the scalar kinetic term (∂ µ u µ ) 2 .Only for the the sigma model the Hamiltonian of the vector field on the Minkowski background is bounded from below.
The Hamiltonian density of a field on a fixed background should not be confused with the Hamiltonian density of the full gravitational theory.The latter is a sum of constraints, as can be seen in the present case in (17), and thus it is always zero.The total energy of the system can be determined in terms of the boundary terms of the Hamiltonian (similarly as for General Relativity [48]).Positivity of the total energy for Einstein-aether theory was shown in [49] for a class of solutions where the vector field is orthogonal to the spatial hypersurfaces of constant time and has a vanishing divergence on such hypersurfaces.We do not consider the boundedness of the total energy in this work.
Instead, similar to the analysis of [47], we consider boundedness of the Hamiltonian of the vector field and the scalar field on the Minkowski background.In other words we fix the metric and ignore its dynamics.Then the Hamiltonian H = d 3 xH is obtained as where we have also gauge fixed the conformal symmetry with Φ = 1, used the second class constraints, and redefinied the potential of the scalar field as Ṽ (φ) = V (φ) + h1(φ) 2f (φ) .Completing the squares in the Hasimilarly miltonian density gives Note that u n is fixed by the constraint (23) as and the definitions of the canonical momenta are now There are two negative terms in the Hamiltonian (34), namely, The former term comes from the Maxwell kinetic term, while the latter is due to the second interaction term between φ and u µ in (9).In the previous subsection we showed that the linear perturbations of the scalar and vector fields do not grow without a limit on this background.Perturbations on other backgrounds might exhibit stability problems.For instance, if there exist a solution where h 2 2 (φ)u i u i grows to become arbitrarily large, the Hamiltonian of the vector and scalar fields is not bounded from below.Thus, solutions with instability might exist.
Stability of spherically symmetric solutions in Einstein-aether theory was studied in [50].It was found that the spherically symmetric perturbations are stable only if the coupling constant of the Maxwell kinetic term satisfies µ 2 > 4, i.e. the coupling in the Lagrangian satisfies − µ 2 4 < −1.Furthermore, in order to avoid noticeable changes in tidal effects it is also necessary to include a further kinetic term or terms in addition to the Maxwell term.The coupling constants of the kinetic terms are severely constrained by the tidal effects.In particular, we could include the scalar kinetic term (∇ µ u µ ) 2 with the value of its coupling constant set either very close to 1 or very close to µ 2 4 .Since MiTeVeS has the similar vector field (with the same kinetic structure and a similar contraint), we have confirmed that the coupling constants of the kinetic terms must be constrained similarly as in Einstein-aether gravity.
Spatially anisotropic cosmology in Einstein-aether theory, where the vector field has a nonvanishing spatial vacuum expectation value, was studied in [51].A longitudinal polarization of the linear perturbation of the vector field was found to become a ghost at the horizon crossing.Anisotropic cosmology has not been considered in MiTeVeS and it is beyond the scope of the present work.

IV. HOMOGENEOUS SOLUTION
In order to find the behavior of the theory for our universe, we consider a flat FLRW metric: and a time-dependent homogeneous ansatz for the scalar field For the vector field we also consider a time-dependent homogeneous ansatz: The energy-momentum tensor ( 13) contains the contribution of a matter fluid or fluids as: The constraint (10b) reduces to Since we assume a homogeneous solution, the kinetic term of the vector field that contains an anti-symmetric part, vanishes F µν = 0. Eq. (10d) with the cosmological background gives: The variation with respect the scalar field φ, Eq. (10c), with a cosmological background, yields: The Einstein tensor with the energy-momentum tensor terms (10a) gives: and (44b) We solve A and λ in terms of φ and φ from ( 41) and ( 42) as where we assume that A and f (φ) are positive, and obtain Ȧ by differentiating (41): .
Then we insert ( 45) and ( 46) into ( 43), (44a) and (44b), which gives: These equations summarize the homogeneous evolution of the universe considering only the scalar field φ as an external degree of freedom besides the scale parameter of the universe.One can check the consistency of the equations using the covariant conservation of the Universe by the term: For the density and the pressure from equations (47b) and (47c) we get directly (47a), when the other matter fields are conserved as well: So we end up with only one simple and consistent set of equations that depends on φ and H. Notice that we can simplify the equations using the redefinition which simplifies the set into: If h 2 (φ) = 0, the set reduces to the canonical equations for a scalar field coupled minimally to gravity: the quintessence model [52][53][54].However, the non-trivial h(φ) term that is present only in the pressure equation, compensate energy density with the other components.Therefore, choosing the function yields the exact solution for the system.

V. DARK ENERGY, DARK MATTER AND STIFF MATTER
In this section we test the impact of the second interaction in (9) with the function h 2 (φ).The interaction couples the derivative of the scalar field φ into the vector field u µ .We consider constant functions: The exact solution of Eq. (47a) for the scalar field is where C 0 > 0 is a constant of integration.Inserting the relation (53) into the Friedmann equations (47b) and (47c) yields: There are three fluid components in addition to the baryonic matter and radiation: dark matter (w = 0), dark energy (w = −1), and stiff matter with equation of state parameter w = 1, The stiff equation of state in cosmology is discussed in [55][56][57].Since the density and pressure of the stiff matter component evolves as a −6 , which is a smaller power than the radiation solution ∼ a −4 , we expect that this part would exist before radiation domination and decays rapidly in our universe.Moreover, some theories of the reheating mechanism require the kination domination era.Therefore, this picture of dark energy, dark matter and stiff matter will predict different scenarios with different potentials.
The coincidence problem rises the question why observable values of dark energy and dark matter densities in the late universe are of the same order of magnitude.The dark matter and the dark energy parts in this model are connected to each other through the parameters f and h 2 .This connection may explain the coincidence problem.
There is a correlation between the shift symmetry and an interaction between the components in the Universe.For constant functions the action has a hidden shift symmetry, that produces a Noether Symmetry: which is covariantly conserved ∇ µ j µ = 0.When we introduce a nonconstant potential V (φ) or the function h(φ), the current is no longer conserved and an interaction between the components emerges.This correspondence was introduced in [58,59] for a unified dark energy and dark matter models with modified measures.But the principle applies here and also for the Mimetic DM models.

VI. COSMOLOGICAL CONSTRAINTS
In order to constrain the additional equation of state we use standard candles (SCs), Baryon Acoustic Oscillations (BAO) and the Cosmic Microwave Background (CMB) distance priors (three data points that contain information of certain features of the power spectrum such as the position of a peak).Radiation is included in order to constraint the components also in the early universe.The normalized Hubble function reads: where H 0 is the current value of the Hubble parameter, and the other constants related to the MiTeVeS model are where Ω φ is the density parameter of the stiff matter, Ω m is the density parameter of usual matter and dark matter, and Ω Λ is the density parameter of dark energy.Ω r is the density parameter of radiation.The matter density  incorporates the dark matter density from the MiTeVeS fields as C 0 h 2 /3 √ f H 2 0 and the density of baryons Ω b .
In order to constrain our model, we use a few data sets.
First we impose the average bound on the possible variation of the Big Bang Nucleosynthesis (BBN) speed-up factor, defined as the ratio of the expansion rate pre- dicted in a given model versus that of the ΛCDM model at the BBN epoch (z BBN ∼ 10 9 ).This amounts to the limit [60] with (∆H/H ΛCDM ) 2 < 0.01, where H ΛCDM is the Hubble rate for the ΛCDM model and the ∆H is the difference between the MiTeVeS Hubble rate and the ΛCDM Hubble rate, ∆H = H − H ΛCDM .Imposing the relation (59) gives a bound on Ω φ : which gives a range of 10 −22 -10 −30 for the Planck values of Ω m and Ω Λ .Therefore we set the value of Ω φ during BBN to be between zero and Ω φ (M ax).
The Cosmic Chronometers (CC) exploit the evolution of differential ages of passive galaxies at different redshifts to directly constrain the Hubble parameter [61].We use uncorrelated 30 CC measurements of H(z) discussed in [62][63][64][65].
For Standard Candles (SC) we use measurements of the Pantheon Type Ia supernova dataset [66] from different binnes.The parameters of the model are fitted by comparing the observed value µ obs i to the theoretical value µ th i of the distance moduli, which is given by: where m and M are the apparent and absolute magnitudes and µ 0 = 5 log H −1 0 /M pc + 25 is the nuisance parameter that has been marginalized.The distance moduli is given for different redshifts µ i = µ(z i ).The luminosity distance is defined by Here, we assume that Ω k = 0 (flat space-time).We use uncorrelated data points from different Baryon Acoustic Oscillations (BAO) [67][68][69][70][71][72][73][74][75][76][77][78].Studies of the BAO feature in the transverse direction provide a measurement of D H (z)/r d = c/H(z)r d , where r d is the sound horizon at the drag epoch and it is taken as an independent parameter and with the comoving angular diameter distance being: In our database we also use the angular diameter distance D A = D M /(1+z) and D V (z)/r d , which is a combination of the BAO peak coordinates above, namely The BAO data represent the absolute distance measurements in the Universe.From the measurements of correlation function or power spectrum of large scale structure, we can use the BAO signal to estimate the distance scales at different redshifts.The BAO data are analyzed based on a fiducial cosmology and the sound horizon at the drag epoch.
For the CMB data, we use the CMB shift parameters [79].The values of the shift parameters with the corresponding covariant forms are called the CMB distant priors.The parameters read: R = Ω m H 0 r(z * )/c, l a = πr(z * )/r s (z * ), (66) which are based on the Planck data release [4].where r s (z) is the comoving sound horizon and z * is the redshift with respect to the photon-decoupling surface.The comoving sound horizon is defined as It is determined by the matter-radiation equality relation and z eq = 2.5 × 10 4 Ω m h 2 (T CMB /2.7K) −4 and with Rb = 31500w b (T CMB /2.7 K) −4 , where the CMB temperature is T CMB = 2.7255 K.
Regarding the problem of likelihood maximization, we use an affine-invariant Markov Chain Monte Carlo (MCMC) sampler, as it is implemented within the opensource package Polychord [80] with the GetDist package [81] to present the results.The prior we choose is with a uniform distribution, where For the baryonic matter we use the value reported by Planck 2018.
The posterior distribution is presented in Fig. 1, which describes the best fit of the MCMC analysis.The differences between the parameters of MiTeVeS and ΛCDM are rather small.For the MiTeVeS the Hubble parameter is 66.52 ± 0.53 km/sec/Mpc and for the ΛCDM model 66.58 ± 0.47 km/sec/Mpc.The matter density is 0.31 ± 0.021 for MiTeVeS and 0.311 ± 0.016 for ΛCDM.The dark energy density is 0.698±0.019for MiTeVeS and 0.698 ± 0.01 for the ΛCDM model.Because of the error bars of all of these parameters, one cannot distinguish the difference from the current cosmological measurements.However, according to current measurements, the stiff matter may exist in our Universe, but its density portion is constrained to be of order 10 −25 or smaller.To be more precise, the MCMC gives: The value is compatible with zero with 2σ.In order to test the preference of the models, we use the Bayesian Evidence that is calculated from the Polychord directly [82].The difference between the Bayesian Evidence yields 0.24, which is an indistinguishable preference for the ΛCDM model.Although the value of the stiff matter density Ω φ has to be very small in order to be consistent with the cosmological observations, being constrained below a maximum value during BBN (61) and having the present value (68) in our model, stiff matter may play a role in the very early universe.Also it could be interesting to mention that a cosmological fluid with the stiff equation of state appears in gravitational theories where the connection has torsion [83-85].Torsion appears in several types of gravitational theories, such as in the Einstein-Cartan and the gauge theories of gravity (for some comprehensive reviews, see, e.g.[86,87] and the references therein) as well as in supergravity (see for example [88]).[89] If the stiff fluid has a negative energy density, which is the the case for example in the Einstein-Cartan theory of gravity [85], there is no initial singularity [56,90], since the scale factor remains nonzero at the beginning.In our model, when there are no interactions between stiff matter, dark matter and dark energy due to the choice of constant functions (52), the energy density of stiff matter is positive (57), and hence the initial singularity exists [56] similarly as in the ΛCDM model.

VII. INTERACTING DARK ENERGY AND DARK MATTER
In this section we discuss the further possibilities of the MiTeVeS theory with general interaction forms.For the solution discussed in sections V and VI the mimetic dark matter component is not distinguishable from dust at the background level.However, for nonconstant functions (Eq.52) one can get interacting dynamical dark energy and dark matter, which emerges naturally from the MiTeVeS action.
In order to track the evolution of the Universe (numerically and analytically) in the general case we rewrite the field equations in terms of the redshift instead of time, which is done using the relation: The modified scalar field equation is given as: where H = H(z) and φ = φ(z) are functions of the redshift.In principle, we need two equations to describe the full evolution.Equations (51b) and (51c) describe the evolution of φ(z) and E(z) = H(z)/H 0 : where Ω m (z) = ρ m /3H 2 0 and Ω p (z) = p m /3H 2 0 .In order to determine the evolution of the variables for a general form of the functions, E(z) and φ(z) are evaluated via these two equations.For the case h(z) = 0 the system recovers to the standard solution, where φ2 ∼ (1 + z) −3 .
A nonconstant function h(φ) can give an interaction between dark energy and dark matter.As an example we take: h(φ) = βφ.The numerical evaluation of the matter density (in red) and the dark energy density (in blue) is described in Fig. 3.The smooth line shows the evaluation without any interaction and the dashed line shows the evaluation with interaction, where we assume β = 1 and φ(z = 0) = 0 for simplicity.The interaction changes the rates between the dark energy and dark matter.The possibility for such an interaction appears naturally in the MiTeVeS action.

VIII. DISCUSSION AND CONCLUSION
In this paper we have uncovered some cosmological implications of the MiTeVeS theory.The particular formulation of the theory includes both a scalar field and a vector field, which is stable and free from ghosts.The two simplest interactions between the scalar field and the vector field were considered.In the homogeneous FLRW background, the first interaction in (9) results only to a redefinition of the potential of the scalar field (50).However, the second interaction in (9) has interesting cosmological implications.For constant functions (52), a conserved Noether current emerges.From the scalar field dynamics, we obtained dark energy, dark matter and stiff matter.When the shift symmetry is broken, the current is not covariantly conserved.Consequently, possible interactions between these three components emerge.
We use the latest observations of the Big Bang Nucleosynthesis, Cosmic Chronometers, Type Ia Supernova, Baryon Acoustic Oscillations and the Cosmic Microwave Background distant prior to constrain the additional stiff matter contribution, which is not present in the standard ΛCDM model.It is clear that the dynamics of the Universe does not change significantly as long as the energy density of the stiff matter part is around 10 −25 or smaller.
There are a few directions to continue the research.It is possible to explore different potentials and functions that approach constant functions asymptotically.We show one simple example of an interaction that changes the rates between dark energy and dark matter naturally from the action, which could be investigated in the future.Such functions produce interaction between the stiff matter, the dark matter and the dark energy parts, and may resolve the cosmic tensions regarding H 0 and others.Moreover, it is important to track the perturbation solutions for this theory in order to compare with further data sets and to validate the formulation in different systems.
The third interesting option is to explore the effect of the scalar Φ in this theory.This scalar is suggested in the theory and gives the behavior of dynamical G ∼ Φ −2 .Many versions of MOND and other modified theories of gravity suggest this effect that naturally emerges in this formulation.So the full theory gives additional effect of emerging gravity that should be studied in the future.

Figure 3 .
Figure 3.The evolution of the Ωm (red) and V (φ) = ΩΛ (blue) for low and high redshifts, without any interaction (smooth line) or with interaction where β = 1 (smooth line).
The posterior distribution for the simplest case of the MiTeVeS with constant functions (red curve) and for ΛCDM model (blue curve).The ratios of the matter density Ωm, dark energy ΩΛ, radiation density Ωr and the stiff part Ω φ summarize to one.The final results for cosmological parameters for the MiTeVeS and the ΛCDM models are summarized in the table.In order to compare the models we calculate the Bayes factor.
The posterior distribution for the simplest case of the MiTeVeS with constant functions (red curve) and for ΛCDM model (blue curve), with H0 vs. Ωm.