Failures of homogeneous and isotropic cosmologies in Extended Quasi-Dilaton Massive Gravity

We analyze the Extended Quasi-Dilaton Massive Gravity model around a Friedmann-Lema\^itre-Robertson-Walker cosmological background. We present a careful stability analysis of asymptotic fixed points. We find that the traditional fixed point cannot be approached dynamically, except from a perfectly fine-tuned initial condition involving both the quasi-dilaton and the Hubble parameter. A less-well examined fixed-point solution, where the time derivative of the 0-th St\"uckelberg field vanishes $\dot\phi^0=0$, encounters no such difficulty, and the fixed point is an attractor in some finite region of initial conditions. We examine the question of the presence of a Boulware-Deser ghost in the theory. We show that the additional constraint which generically allows for the elimination of the Boulware-Deser mode is only present under special initial conditions. We find that the only possibility corresponds to the traditional fixed point and the initial conditions are the same fine-tuned conditions that allow the fixed point to be approached dynamically.


I. INTRODUCTION
In the standard cosmological model, the accelerated expansion of the universe is attributed to the cosmological constant Λ. However, to match the observed expansion, Λ must be of the order of 10 −122 in Planck units, which raises a fine-tuning problem. A possible alternative is to modify general relativity (GR) at large distances or low momenta. A massive spin-2 field theory, known as the dRGT theory [1,2], is a theoretically well-motivated modification of GR. However, the dRGT theory does not admit a flat Friedmann-Lemaître-Robertson-Walker (FLRW) solution with an expanding scale factor [3]. A modification to dRGT gravity, known as Quasi-Dilaton Massive Gravity (QDMG), was proposed in [4] and provides homogeneous and isotropic expanding solutions. It was later shown [5] that the parameters of QDMG have to be finely tuned in order to match the observed expansion history of the universe. More disastrously, the results in [6,7] indicate that the scalar perturbations in QDMG acquire a wrong sign kinetic term at short scales. A further modification, Extended QDMG (EQDMG), was proposed in [8] and has scalar perturbations that are thought to be stable at all momentum scales. The standard fixedpoint cosmological solution of EQDMG has a de Sitter metric, and thus appears to be a good candidate for latetime cosmology.
EQDMG only differs from QDMG by the addition to * Electronic address: stefano.anselmi@iap.fr † Electronic address: saurabh.kumar@case.edu ‡ Electronic address: diana.laura.lopez.nacir@cern.ch § Electronic address: glenn.starkman@case.edu the action of one operator involving the Quasi-Dilaton (QD) field and a new free parameter α σ . Naively, in the limit α σ → 0, EQDMG reduces to QDMG, but actually the limit is very subtle. Indeed, there are controversial results in the literature regarding whether or not EQDMG contains an unavoidable additional degree of freedom, which would correspond to the Boulware-Deser ghost (BD) [9][10][11].
In this paper, for the first time, we: 1. assess the stability of the standard fixed-point solutions (referred to as Case 1 in this paper) and show that this assessment requires a non-standard approach; 2. demonstrate that the Case 1 fixed points cannot be approached dynamically, due to an unavoidable singularity in the dynamical equations; 3. perform a comprehensive study of a new branch of solutions (referred to as Case 2), first proposed by [12] but largely ignored in the literature, and show that it provides stable and dynamically attainable fixed-point solutions; 4. show for a flat FLRW universe, the fact that the background equations are satisfied does not guarantee the presence of the additional constraint necessary to eliminate the BD mode (in agreement with the results of [9], but in disagreement with the computations in [11]); 5. find that the only branch of solutions for which the additional constraint exists corresponds to Case 1; 6. argue that, in order to avoid a BD ghost, the initial values of certain EQDMG dynamical variables must be extremely fine-tuned; 7. verify that the same fine-tuned initial conditions also allow the fixed point to be approached dynamically.
The paper is organized as follows. In Section II and III, we summarize the theory of EQDMG.
In Section IV, we define the dynamical variables and provide the relevant background equations.
In Section V, we find the fixed-point and de Sitter solutions of the dynamical equations, and show that they are equivalent to one another (provided the Hubble rate is positive). We identify four independent fixed-point cases, each of which is studied in greater detail in the sections that follow.
In Section VI, we introduce the standard procedure for analyzing the stability of the fixed-point solutions for the background. We discuss the inadequacy of this procedure for the Case 1 fixed points and provide an augmented framework.
In Section VII, we present the results of our numerical search for viable parameters for the EQDMG theory. We find that except for a very specific, precisely fine-tuned initial displacement away from the fixed-point values of Case 1 (explained in Section VIII), the fixed points cannot be reached in the asymptotic future.
In Section VIII, we further study the perturbative stability of the scalar sector of the theory, both in the vacuum case and with matter. We identify the conditions on the dynamical variables required to avoid the BD ghost, and show that the fact the background equations are satisfied does not guarantee the validity of the conditions. Cases 1 turns out to be the only case for which the additional constraint necessary to eliminate the BD mode can be obtained by an appropriate choice of initial conditions. However, those conditions appear to represent a difficult fine-tuning of all the degrees of freedom.
In Section IX we present our conclusions. We provide some detailed calculations and consider the special case of Minkowski solutions in the appendices.

II. FORMALISM
We consider the action for the extended quasi-dilaton theory [8]: where M Pl is the Planck mass and, in addition to the Einstein Hilbert action S EH , we have the contribution S m of the matter sector, and a quasi-dilaton contribution S σ . Here with square brackets denoting a trace, and S σ includes five new fields: σ is the quasi-dilaton scalar field and φ a (a = 0, · · · , 3) are the four Stückelberg fields. It also depends on the coupling constants α σ , α 2 and α 3 , and on the graviton mass m g . For α σ = 0 one recovers the standard quasi-dilaton theory.
In the space of Stückelberg fields, the theory enjoys the Poincare symmetry and a global symmetry given by with σ 0 an arbitrary constant.

III. BACKGROUND
We consider a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) ansatz, for which The fiducial metric f µν reduces to where The minisuperspace action for the background can be written as and we have defined It is worth pointing out here that in (5), and we make the (+, +, +, +) choice following [13]. This gives us r = 0 represents a determinant singularity in either or both f µν (if n = 0) or g µν (if a = 0) -a spacelike hypersurface where the dimensionality of the metric changes. The stability of the theory across that hypersurface is unclear [14]. Indeed, we find that when r approaches too close to zero, our numerical integrations of the dynamical equations become unstable; the instability in the numerical noise may be due to an underlying instability in the theory. We insist here that r = 0.

IV. DYNAMICAL EQUATIONS AND VARIABLES
We next set out the dynamical variables describing the background, and the dynamical equations describing their time evolution.
Varying the action with respect to φ 0 (t) leads to the constraint equation where This suggests that it will be useful to introduce as one of the dynamical variables The solution of (22) is immediately We see that, in any reasonable cosmological context, y → 0 in the asymptotic future, and that this can be achieved by one (or more) of four quantities approaching or equalling zero: J(X(t)),φ 0 (t), X(t), or X(t) − 1. These four cases will drive our analysis.
Anticipating that it will be convenient to regard X as a dynamical variable, we differentiate (19) with respect to time to getẊ H(t) ≡ȧ/a is the Hubble parameter.
Varying the action with respect to the lapse N (t) and using time reparametrization invariance to set N (t) = 1, we obtain the Friedmann equation: • From (22), we have immediately • The equations forΩ m andΩ r are similarly easily obtained from (31) and (32): • Equation (26) can be rewritten using (29), (30) and the definition ofΩ DE : The ± represents the possibility thatσ can be positive or negative 3 .
• Equation (33) can be rewritten using (29) and (30): In the above set of equations one must replace •h 2 by the Friedmann equation, which now takes the simple form • and (combining (14), (20) and (25)) r with 4 The argument of the square root on the right-hand side of (40) must be positive for r to be real. The reality of r is a condition on the dynamical variables that must be checked, in case (as we find below) it is not automatically satisfied. In particular, we see that problems may arise if (y/G 2 (X)) 2 → 1. 3 We will focus our attention below on the positive sign, because the negative sign leads to only an X = 0 fixed point. 4 For ασ = 0, r cannot be determined from (40), because (24) gets reduced to y = G 2 (X). Thus, we can no longer use (38), and the above system of evolution equations is not well-equipped to handle this case. In fact, this limit gives us the Quasi-Dilaton theory and the evolution of the dynamical variables have been previously studied by [5] and [15].

V. FIXED-POINT ANALYSIS
In this section, we evaluate the dynamical variables when their N −derivatives vanish in (34)-(38). We term the values of the dynamical variables in this limit as fixed points.
In the fixed-point limit, the left-hand sides of equations (34)-(36) vanish, giving us y F P = 0, and alsõ Ω m,F P =Ω r,F P = 0. From (39), we learn that Ω DE,F P =h 2 F P . The solutions to (34)-(36) are where y 0 ,Ω m0 andΩ r0 are the corresponding initial values. Thus fixed points occur in the asymptotic future, i.e. as N → ∞ and so a → ∞.
The left-hand side of (37) vanishes at the fixed point, implying that X F P is a constant. If X F P = 0, the righthand side of (37) (and if X F P = 0, then the right-hand side of (38)) provides us with , X F P = 0 , at the fixed point. In arriving at (42), we take the '+' sign in (37), since the '−' sign leads to X = 0 as the only fixed-point solution.
Notice that, as for getting background fixed-point solutions, both 0 < ω < 6 and ω ≥ 6 are in principle suitable regions in the parameter space, since G 1 (X F P ) can be either positive or negative. A special value is ω = 6, in which case (37) and (42) demand G 1 (X F P ) = 0.
Observing thath 2 is also a constant at the fixed point, we conclude that the fixed points of the evolution equations are de Sitter.
We find that the converse is also true: the de Sitter solutions of the evolution equations are fixed points as we approach the asymptotic future. To prove this, we require that the dynamical variables attain the following de Sitter values in the future, whereh 2 F P is a constant different from zero. In this situation, the left-hand sides of (38), (35) and (36) become zero, meaning they are fixed points.
From (25), we learn that y = 0 in the asymptotic future, which means that the left-hand side of (34) is also zero. The only point left to establish is that X approaches a constant in the future.
From the definitions (14), (20) and (24), one can split the fixed-point solution into four cases: • Case 1 (the standard case): and hence the fixed-point solutions are • Case 2: Since the left-hand side of (38) vanishes, (40) provides us with the following equation in X F P Squaring (47) gives us a polynomial equation for X F P . X F P can be any of the roots of that polynomial.
• Case 3: • Case 4: For all cases, X approaches a constant and thus the left-hand side of (37) vanishes in the asymptotic future proving it is a fixed point.
We analyze the fixed-point solutions in more detail below.
A. Case 1: J(XF P ) = 0 As will become clear below, this case is very subtle. Note that these fixed points are the same as ones analyzed in [5], [15] for the QD theory (i.e., the EQD with α σ = 0).
Requiring that X ± be real, means 5 (51) 5 A special case occurs when α 4 = 2α 2 3 3 , in which case G 1 (X) = 0 = J(X). We are left with (X F P , r, F P , ω) = 1 + 3 This solution for ω = 6 andh 2 is indeterminate because of simplifications occurring in (37). Therefore, on the α 4 = From (38) and (42), we get the same expression as in the QD theory for the fixed-point limit of r (assuming ω = 6): Note that this expression is valid for Cases 1, 2 and 4 (not X = 0), provided ω = 6 and G 1 (X) = 0. Unlike QD theory, since α σ = 0 we can use (24) to obtain (40), which gives r in terms of the dynamical variables. If the system is to evolve towards its expected fixed point, r must approach r F P . Therefore, at the fixed point In order for this approach to be smooth, one also needs One possibility is that z → 0, which requires a fine tuned relation among the parameters of the model. Since this is subsumed in Case 2 anyway, we will not analyze this particular case any further. A second possibility is that the quantity in square brackets approaches zero in the fixed-point limit. This implies a constraint equation for the dynamical variables near the fixed point. We note that in previous literature this constraint equation has been assumed to hold for the full dynamics with no justification 6 . In such a case, the evolution of the dynamical system near the fixed point could be described in terms of 4 instead of 5 dynamical variables; i.e., near the fixedpoint limit the evolution would be driven by the same dynamical equations as in the QD theory. However, a priori there is no reason to expect this condition to be valid for α σ = 0.
B. Case 2: For Case 2, one has to solve (46) to get the fixed-point value of X. We must take care that, after solving for X F P , the sign of 6(h 2 F P − G 1 (X F P )) + X F P G 1 (X F P ) should be the same as the sign of G 1 (X F P ) 7 .
It is worthwhile noting that for Case 2, α σ > 0. (We omit the QD case, α σ = 0.) This can be seen by inspection of equation (46), recalling that From (42) and the definition ofΩ DE , (29) and (30), we getσ We find that r F P is indeterminate from both (38) and (40). Examining (19), we see that there are two ways to get X F P = 0. The first possibility is that σ → −∞, in which case we can draw no conclusion from (14)  Since G 1 (1) = 0 and G 1 (1) = 1, and since Ω m = Ω r = 0 at the fixed point, from (37) and (39), we find that X F P = 1 requires ω = 6. Substituting X = X F P = 1 in (38), we get whereh F P is indeterminate. Therefore the theory loses its predictive power. For this reason we will not consider this case anymore in the following analysis.

VI. FIXED-POINT LINEAR STABILITY
We wish to check whether or not the fixed-point solutions are attractors in the asymptotic future. This would be the case if any small perturbation around the fixed point decays to zero asymptotically.
We start with the prescription given by [17] to evaluate the fixed-point stability.
Let V = y, X,Ω DE ,Ω m ,Ω r T denote the dynamical variables and f (V) be the RHS of the first order differential equations. Thus we can express equations (34)-(38) as Assuming small perturbations δV around any point, V 0 , a Taylor expansion of the functions f (V) gives us where M is the stability matrix. Its elements are given by In Appendix A, we provide the analytical expressions for the elements of M (Equations (A1)-(A7)).
If V 0 is a fixed point, then the second term of RHS of (58) would vanish. Using the eigenvectors of M, one can then find matrix P such that where (λ 1 , . . . , λ 5 ) are the eigenvalues given by (A16). We define δV as Thus the solution to (58) in the new basis would be where C i 's are integration constants. Multiplying the above equation by P and thereby returning to the original basis, we get Then, for the fixed points to be stable, we require δV to approach zero as N → ∞. It can be seen from (64) that if the eigenvalues of M are either real and negative or imaginary with negative real part, the fixed points will be stable or form a stable spiral respectively. We find in Appendix A, that to obtain attractor solutions, λ 4 and λ 5 must be real and negative or complex with negative real parts. This requires the elements of M to satisfy the condition for stable solutions or for stable spiral solutions.
A. Case 1: J(XF P ) = 0 The above-described standard method will not suffice to evaluate fixed-point stability in all cases. In particular, for J(X F P ) = 0 the stability matrix M has divergent terms at the fixed point (see Appendix A for the exact expressions of its elements). This makes the evaluation of the matrix P and its inverse indeterminate. One can thus no longer diagonalize the system of equations as in (60) and come up with solutions for (58) given by (64).
We devise the following scheme to assess stability: we introduce small perturbations δV around an arbitrary point V 0 ; using (58) and diagonalizing M, we solve for δV. In order for the fixed point to be an attractor we require the following conditions to be satisfied. If V 0 is infinitesimally close to the fixed point V F P , then, as N → ∞: (A) the perturbations δV are infinitesimally close to zero, therefore we require (B) the derivatives of perturbations, d δV/dN , are infinitesimally close to zero, therefore we must verify that In more detail, after diagonalizing the matrix M, from (58) we obtain where Solving (70), we are left with where C is the integration constant vector. Upon returning to the original basis, we find + P Diag e λ1N , · · · , e λ5N C.
In Appendix B we show that the requirement (A) could be satisfied if the eigenvalues λ i 's are either real and negative or complex with negative real part. It is also shown that, if X 0 is a point infinitesimally close to its fixed-point value X F P , it must satisfy the following relations 8 In the numerical investigation performed in Section VII, we will take X 0 as the initial condition for the dynamical variable X. Therefore, relations (74)-(77) provide the viable initial conditions needed to have linear stable solutions for the dynamical variables. We note here that (74)-(77) can only be satisfied for X 0 either greater or less than X F P and never both 9 . For V FP to be an attractor, we are required to verify under which conditions (B) and (C) are satisfied. However, as explained in more details in Appendix B, in linear perturbation theory, (B) and (C) cannot be determined. Therefore a numerical investigation on both (68) and (69)  Recalling (55), we find that the presence of terms X Ωσ and y G2(X) in the stability matrix make the ele- 8 In special cases where G 1 (X F P ) = 0, such as ω = 6, or 2α 2 3 = α 4 , the eigenvalues involve ratios of zeros that we are unable to resolve, so we cannot determine the stability. 9 We exclude particular cases for which in a neighborhood of X F P not only G 2 (X) changes sign but also either G 2 (X) or G 1 (X). ments (A1), (A2), (A5), (A6) and (A7) and consequently the eigenvalues (A16) indeterminate. Hence the stability of X F P = 0 is unclear.

VII. NUMERICAL INVESTIGATION OF FIXED-POINT STABILITY
As noted above, fixed-point linear stability conditions are necessary but not sufficient to guarantee that a fixed point can be reached by evolving from an initial configuration that is displaced from that fixed point. Moreover, as described above, for Case 1 and Case 3 some or all the relevant quantities that appear in the background equations cannot be analytically assessed in linear perturbation theory. Therefore we need to perform numerical tests of the fixed-points stability.
For suitable values of the parameters (α 3 , α 4 , ω, α σ ), we check numerically that the 5 dynamical variables approach their fixed-point values if initially perturbed from them. Crucially, this includes verifying that z = y/G 2 (X) approaches the value given by (53).
In our numerical investigations, we use the results of Section (VI) to guarantee linear stability. We look further for values of the parameters for which: (42)  Note that we could further constrain the parameter space by selecting regions where the scalar, vector and tensor perturbations are stable [8,11,12]. Although it is not necessary for the present analysis, to simplify our search, we use some of the necessary restrictions imposed by the stability of the perturbations, which we summarize in Appendix E.
As discussed below, we find that the traditional (Case 1) fixed points cannot be reached from any neighboring configuration because z does not approach the corresponding fixed-point value! The Case 2 fixed points behaved as expected from the linear stability analysis. The Case 3 fixed points always encounter a singularity before reaching the fixed-point values.
A. Case 1: J(XF P ) = 0 Our goal in this section is nominally to identify values of the parameters (α 3 , α 4 , ω, α σ ) such that (i)-(iii) hold true and verify that the dynamical variables approach their fixed-point values in the asymptotic future. One of the central results of this paper is that we fail to do so. We show that, for all the choices of parameters and initial conditions we consider, the dynamical variables never approach their fixed-point values if they do not start there to begin with.
We remind the reader that since G 2 (X) vanishes at the fixed point, the stability conditions (74)-(77) only hold true in the vicinity of the fixed point. In this example, since 0 < ω < 6 and α σ > 0, the initial value of X must be of the form X 0 = X + + , where 0 < 1 so that (74) is satisfied. To study the behavior of the dynamical variables close to the fixed point, we set = 10 −6 and the initial conditions to beΩ DE,0 =Ω DE,F P + ,Ω m,0 = , Ω r,0 = . We recall that for Case 1, we cannot freely perturb y, because the initial value of r from (40) would not then be close to its fixed-point value (52). Instead, we perturb r by and initialize y using (B49). For a detailed discussion, we refer the reader to Appendix B. Notice thatΩ DE,0 , r 0 , y 0 could be either greater or less than the fixed-point values. In this example we chose the former.
After setting the initial conditions as described above, we study the behavior of the dynamical variables with time. Equations (34), (35) and (36) have simple solutions, The evolution equations for X andΩ DE have no analytical solutions. Therefore, we use (78) in (39) and (40), and solve equations (37) and (38) numerically. We find that (37) and (38) evolve until they reach a singularity. As shown in Fig. 1, the evolution of X andΩ DE before the singularity turns out to be exactly what we expect from perturbation theory: they oscillate around their fixed-point values with decaying amplitude and increasing frequency.
By inspecting the RHSs of (37) and (38)  To ensure that the phenomenon described above is not a numerical artifact, we performed numerous numerical and analytic tests. Numerically, we confirmed that the z 2 → 1 behavior was robust to: increasing the numerical precision demanded from the integrator, increasing the order of the integrator, changing the integration scheme, and changing the initial conditions. Therefore, we have conclusive evidence that the evolution of the dynamical variables drives z 2 towards unity. We therefore conclude that for this value of the EQDMG parameters, in the neighborhood of a Case 1 fixed point the evolution equations (37) and (38) reach a singularity in the asymptotic future.
In order to confirm that this conclusion holds true for all regions of the parameter space we would need  to evolve the differential equations for a high number of different values of (α 3 , α 4 , ω, α σ ) and initial conditions close to the fixed-point solutions. This is numerically too costly, therefore we solve the exact equations for a limited number of points in the parameter space (few thousands). For a more extensive scan (millions of points) we rely on a sensible approximation based on perturbation theory and on the following argument. If the fixed point is an attractor, perturbation theory should be increasingly accurate the closer the initial conditions are to the fixed-point value; if the fixed point is not an attractor, then perturbation theory may or may not work. If in perturbation theory we were to find that z approached its fixed-point value, then we would want to verify that conclusion with the full non-perturbative solution. If z does not approach its fixed-point value in perturbation theory, then that may be because the initial perturbation away from the fixed point was too largeoutside the basin of attraction -so one should decrease the magnitude of the perturbation as much as possible. The smaller the perturbation, the more one would expect to trust perturbation theory. If this still fails, it is always possible that the basin of attraction is extremely small and difficult to find numerically; but in any case, while we will not have arrived at a mathematical proof that this fixed point is not an attractor, we will certainly have shown that it is not a suitable candidate for a cosmological model.
We start by testing the prediction of z 2 from perturbation theory in the previously considered example. We rely on equations (B25), (B26), (B27), (B28), (B29) where the C i coefficients are determined by requiring that δV(N in ≡ 0) = 0. In Fig. 2 we plot the behavior of z 2 for the first ∼ 20 periods comparing the full and perturbation theory solutions. We verify a consistent growth in the peaks of z 2 . Hence, in the subsequent analysis, we employ the following supporting argument: if in perturbation theory the first 20 peaks of z 2 grow monotonically with time, we conclude that z 2 does not approach its fixed-point value. We randomly selected ∼ 10 7 values of the parameters in the following ranges 10 −5 < |α 3 | , |α 4 | , α σ < 10 5 , 0 < ω < 6 or 6 < ω < 10 5 and evaluate X ± . We evaluate the first ∼ 20 peaks of z 2 for all points that satisfied the conditions (i)-(iii) together with the pertinent ones in Appendix E using linear perturbation theory. We set the initial conditions of the dynamical variables to be from = 10 −8 to = 10 −12 away from their fixed-point values 10 . We repeated the same procedure using the numerical solutions to the full (nonperturbative) equations with few thousand randomly selected points. We find that the peaks of z 2 always grow monotonically in both linear perturbation theory and using the exact equations.
Hence we come to the following conclusion: the standard (Case 1) fixed-point solutions of Extended Quasi-Dilaton Massive Gravity are dynamically unattainable due to an unavoidable singularity while approaching them. Hence, we rule out the suitability of this fixed point as a cosmological model. Recalling that α σ is positive, we randomly select points in the ranges (79).
To obtain the fixed-point values of X, we use (47) and (42) to arrive at a 7th order polynomial in X with no analytical solutions. After numerically evaluating the roots of the polynomial, we demand that the fixed-point values satisfy (i)-(iii) and that the vector and tensor perturbations are stable 11 (see Appendix E for details). We find parameter values for which both stable and spiral solutions are allowed. We selected a few points in the allowed parameter space to verify the attractor behavior of the fixed points. Starting from small perturbations ( = 10 −6 ) around the fixed-point values of the dynamical variables, we study the evolution of the differential equations with time by solving (34), (35) and (36) analytically and (37) and (38) numerically. As predicted from linear stability analysis (Section VI), we find that the dynamical variables reach their fixed-point values in the asymptotic future. C. Case 3: XF P = 0 For the Case 3 fixed points we do not have indications from the linear stability analysis. We perform a random search selecting 64 × 10 4 points in the parameter space. Starting from small perturbations ( = 10 −6 ) around the fixed-point values we find that the dynamical variables are approaching the fixed-point limit. However, we find that the quantity Ω DE − G 1 (X) always approaches zero at finite time. We find numerically that this makes r vanish which is not allowed in the theory.

VIII. PERTURBATIONS: AVOIDING THE BOULWARE-DESER GHOST
One of the outstanding concerns in any theory that adds a scalar degree of freedom to Einstein-Hilbert gravity is the possibility that the theory includes a Boulware-Deser ghost -a dynamical degree of freedom with a wrong-sign kinetic term.
Due to disagreements in the literature mentioned above (see [9][10][11]), in this section we reconsider the analysis of the scalar perturbations, and determine under what conditions there is or not a necessary additional constraint equation that generically allows one to eliminate the Boulware-Deser mode.
We will show that such a constraint exists when J(X) = 0 (Case 1) orσ = 0, giving the possibility of being ghost-free at the corresponding fixed points. In both cases, there is also a well-defined set of initial conditions of the dynamical variables for which this virtue extends beyond the fixed point. The first case is the only one potentially relevant for cosmology, and we explore the consequences of the initial conditions below. The second caseσ = 0 corresponds to a Minkowski background metric (provided X = 0) which we consider only in Appendix C.
Given that in the previous section, for Case 3 we have shown X = 0 cannot be approached dynamically from any neighborhood, one might consider the possibility when X(N ) = 0 ∀ N > 0. This possibility can also be ruled out because X = 0 only makes sense as a fixedpoint limit. Hence from now on, we will not consider Case 3 any further.
Following the standard treatment and in this section only, we take the action of the matter sector to be 12 which corresponds to the addition of a scalar field χ with a non-canonical kinetic term, given by the function The fluid variables (pressure p, energy density ρ and sound speed c s ) associated with χ can be written as [18]: To study the perturbations in χ, we replace χ by χ + M Pl δχ.
To facilitate comparisons, we adopt notations as close as possible to those most used in the literature. We decompose the metric into tensor, vector, and scalar as For the sake of clarity, in this section we present only those results that are either more generic or different from those in previous literature, and relegate the remaining details to Appendix D.
We focus first in the caseφ 0 = 0 and, following literature, we adopt the unitary gauge φ a = φ 0 (t)δ a 0 + δ a i x i . The caseφ 0 (t) = 0, for which we cannot work in this gauge, is analyzed later.
The part of the action (2) that is quadratic in the perturbations can be split into tensor, vector and scalar contributions. We focus here on the scalar part.
From the variation of the quadratic action for the scalar sector with respect to Φ and B, we obtain constraint equations that allow us to eliminate these two variables. The solution for Φ and B can be found for a generic background and without making use of the background equations 13 . Introducing the solutions for Φ and B back into the action, we can write the kinetic part as 13 There are some disagreements in the literature regarding these results. In the appropriate limit, our results reduce to those of [12] rather than those of [11]. For more details see Appendix D.
with Z = {ψ, δσ, E, δχ}. One combination corresponds to the Boulware-Deser mode. For this mode to be nondynamical, the determinant of the matrix Q must vanish. We first consider the vacuum case, and then see how the inclusion of matter changes the results.
• Vacuum case: In the absence of the additional field χ, Q is a 3 × 3 matrix and Z = {ψ, δσ, E}. The determinant Det(Q) can be computed analytically in Fourier space. After expanding in powers of comoving wavenumber k, keeping only the leading order terms (the infrared part), and using the background equations to express H andḢ in terms of the other dynamical variables, The fact that the background equations are satisfied does not guarantee that the determinant vanishes, as previous computations suggested [11]. However, it is clear that Det(Q) vanishes when J(X) = 0 orσ = 0. Moreover, in that cases it can be shown that the determinant vanishes at all order in k. As we show next, these results are robust to the addition of matter.
• With Matter: Proceeding as above, at leading order in the wavenumber k, after using the background equations to replaceḢ and H in terms of the other dynamical variables, the determinant of the now 4 × 4 matrix Q can be written as Therefore, ignoring the particular cases for which X = 1, we see that the determinant vanishes in the infrared under the same conditions as in the vacuum case. Under those conditions, we also check the determinant vanishes at all order in k.
Of course the determinant also vanishes in the case where the matter is just a cosmological constant, P ,Y = 0.
Therefore, Det(Q) = 0 at the Case 1 fixed points. If we perturb X away from that fixed-point value, then J(X) = 0, and Det(Q) = 0 either. However, from (22), we can conclude that if J ≡ J(X) = 0 at some initial time t 0 , then alsoJ = 0 at that time, and alsoJ. Thus J(X) remains zero at all times once X = X F P . Consequently Det(Q) also remains equal to zero. Therefore, in order to eliminate the BD ghost in both the vacuum and matter contexts, we must impose special initial conditions on the dynamical variables; namely, that X = X ± is exactly satisfied. Setting X to this fixed-point value, and thus setting J(X) = 0, y = 0 andẊ = 0, we can solve (35) through (40) to obtain: This appears to be just a fine-tuning of the dynamical variable X to some parameter-dependent value; however X is a function of both the quasi-dilaton and the scale factor given by (19). This thus appears to be an awkward fine-tuning, relating the initial values of many of the dynamical variables to one another. It also allows the fixed points to become stable attractors that can be approached dynamically. Furthermore, notice that in the asymptotic past, N → −∞, the matter and radiation terms will dominate over G 1 (X±). This restricts the values of ω to 0 < ω < 6 or it means that the theory is not valid arbitrarily far into the past.
B.φ 0 = 0 As mentioned earlier, in the special caseφ 0 = 0 (Case 2) we cannot use the unitary gauge. Assuming H = 0 we choose the gauge with ψ = 0 instead, while we keep φ i = x i . In principle, as first noticed by [12], this case could also be interesting beyond the fixed-point limit. This corresponds to settingφ 0 (t 0 ) = 0 at some initial time t 0 , yieldingφ 0 (t 0 ) = 0 andφ 0 (t) = 0 for all time. This enforces y = 0, but the other dynamical variables remain free to evolve.
The perturbations of the Stückelberg scalar degree of freedom δφ 0 (sinceφ 0 = 0 at the background level) only enter as a contribution to f µν that is quadratic in δφ 0 . Therefore, the quadratic action for this scalar decouples from the other parts, and the kinetic part can be immediately computed Now, to compute the determinant of the kinetic matrix corresponding to the other degrees of freedom we proceed as above. We integrate out Φ and B, we write the relevant part of the action as in (84), and we consider the vacuum case and the case with matter separately: • Vacuum case: In this case Z = {δσ, E}, and Det(Q) = G 1 (X)k 4 m 2 g Xωa 2 6α σ G 1 (X) + r 2 X 2 ω × 4 9α σ G 1 (X)G 1 (X)m 2 g Xa 2 + k 2 (r + 1) (6α σ G 1 (X) • With matter: Z = {δσ, E, δχ}, and Hence, we conclude that the conditionφ 0 = 0 is not sufficient to obtain the required additional constraint.

IX. CONCLUSIONS
In this paper we studied the fixed-point solutions of EQDMG in great detail after splitting them into four cases.
We performed a linear stability analysis of the background (homogeneous) fixed-point solutions. This stability analysis for the standard case (Case 1) required an unconventional approach. We derived necessary stability conditions for the dynamical variables. However, we verify numerically that the dynamical variables inevitably encounter a singularity while approaching their fixed-point values. This is because a function of two of the dynamical variables fails to approach its fixed-point limit, and instead oscillates, finally reaching a critical value at which the dynamical equations are no longer valid.
On the other hand, in Case 2, the dynamical variables smoothly asymptote towards their fixed-points values. However, in this case, the presence of the additional constraint that would allow one to eliminate the BD mode is not guaranteed. A numerical search of the Case 2 parameter space revealed many values of the parameters that have stable evolutions toward the fixed point.
We analyzed the conditions under which the constraint equation that generically allows for the elimination of the Boulware-Deser ghost can be obtained. We found that the constraint equation exists in Case 1 type of solutions, but not around a generic background. Moreover, in these solutions the fixed points are attractors. However, such solutions require the time-derivative of the quasi-dilaton field must be exactly tuned against the Hubble parameter.
The conclusive result of our study indicates that the EQDMG theory shows pathological behaviors when a generic FLRW solution is assumed. Only an "awkward" fine-tuned solution is healthy. From our point of view this finding makes the EQDMG less appealing as a viable model to explain the evolution of our Universe.
The extensive analytical and numerical analysis presented for the EQDMG theory must be carried out for all the proposed massive gravity theories that provide flat Friedmann-Lemaître-Robertson-Walker solutions (a non-exhaustive list includes [9,[19][20][21][22]).
As mentioned below (87), the condition ω < 6 arises in EQDMG when the J(X) = 0 solution is taken to describe the past cosmological evolution, when matter and radiation dominated. In EQDMG this solution is imposed to avoid the presence of the BD ghost. However, since in QDMG the BD mode can be eliminated without setting J(X) = 0, the restriction ω < 6 is in principle unnecessary. It would be worth exploring the parameter region ω > 6 to see whether the simpler QDMG theory allows a proper description of the expansion history of the Universe. We propose to contribute to the understanding of this issue in the future.

General expression
The derivative/stability matrix M defined in (59) takes the form In the above equations, one must replaceh 2 by (39) and Ω σ byΩ DE − G 1 (X).

Case 2
We evaluate M for the Case 2 fixed point. Notice thatΩ σ , when X = 0 (using (26)), is given bỹ As explained in Section V, in Case 2 we must take the upper sign in (37), this implies that in (A1)-(A4) the sign is also positive. The elements of the M matrix are given by at the fixed point. In the above equations, one must 14 use (42) to replaceh 2 by G1(X) (1− ω /6) . The eigenvalues (λ i 's) of the matrix are: (A17) We can thus impose the following conditions on elements (A9)-(A15) so that the fixed-points become stable or stable spirals: • Stable spiral solutions These are equations (65) and (66) in Section VI.
14 Except in the particular case where ω = 6, in which case one must instead use a root of In this appendix we provide the detailed computations supporting the linear stability analysis presented in Section VI A and the numerical analysis explained in Section VII A for the Case 1 fixed points.
In order to compute δV from (73) and the limits (68) and (69), we split the analysis in two steps: 1. we provide the analytical expression for the matrix P and its inverse and analyze the fixed-point limit of their elements; 2. we expand (73) and show that, under certain conditions, perturbations are infinitesimally close to zero in the limit N → +∞ and V 0 → V F P . We then show that in linear perturbation theory the limits (68) and (69) cannot be assessed.
The third subsection of this appendix deals with the initial conditions needed in numerical analysis presented in Section VII.

Matrix P and its inverse in the fixed-point limit
In order to analyze the fixed-point limit of δV we previously need to study the matrix P and its inverse.
The column vectors in the matrix P are composed of eigenvectors of M. The matrix P reads P =      P 11 P 12 0 0 0 P 21 P 22 P 23 P 24 P 25 0 1 P 33 1 1 0 We first need to address the fixed-point limit of the matrix M elements. We recall that in Case 1, although both y and G 2 (X) vanish, the ratio y/G 2 (X) is finite and is given by (53). Therefore, one can see that the elements of the matrix, M 31 and M 32 ((A5) and (A6) respectively) are divergent in the fixed-point limit 15 , since both of them scale as (∼ 1/y). All the other elements of M are finite.
The non-zero elements of P which have the terms M 31 and M 32 , either depend on M −1 31 , M −1 32 (both approaching 0) or contain the ratio M 32 /M 31 which tends to at the fixed point. Therefore the elements of P are convergent at the fixed point. We now show that some of the elements of P −1 are divergent. The matrix P −1 reads where It is easy to see that D 1 , D 3 ∼ 1 /y and D 2 , N − ,N + ∼ 1 / √ y which makes R 41 , R 42 , R 51 and R 52 divergent at the fixed point.

Analysis of perturbations in the fixed-point limit
The perturbations δV given by (73) are expanded to Since δX 0 and δΩ DE0 are real quantities, we must set the imaginary terms in (B26) and (B27) occurring in the last two terms to zero. From (B27), this amounts to C 4 = C * 5 . Additionally, it can be shown that this requirement also makes δX real automatically.
We now study the limit of δV for N → +∞ and V 0 → V F P . We identify two different possibilities.
(a) If and only if the eigenvalues λ 4 and λ 5 are either real and negative or complex with negative real part, none of the two limits give rise to divergent terms. Thus the two limits commute and the final limit is well defined. Recalling that We conclude that, under the above specified conditions for the eigenvalues, the perturbations approach zero near the fixed point.
(b) If the eigenvalues λ 4 and λ 5 are either real and positive (or complex with positive real part) the last term of (B27) diverges (or is indeterminate) for N → ∞ while the other terms remain finite. We conclude that δΩ DE is divergent (or indeterminate).
From a close inspection to λ 4 , λ 5 given by (A16) we find that, using (53), M 22 + M 33 = −3 (as in (A17)). Thus, being M 32 divergent, condition (a) is realized only if M 23 M 32 < 0. We notice that at the fixed point which is positive if 0 < ω < 6 and negative when ω > 6. This means one must satisfy the following condition for stable attractors: To satisfy the above conditions at the fixed point, we first evaluate M 32 in (A6) at a point X 0 that is infinitesimally close to the fixed-point value of X. Demanding that M 32 < 0, 0 < ω < 6 > 0, ω > 6 (B33) at the fixed-point limit, (B32) will automatically be satisfied. Evaluating (A6) at X 0 , we find that the divergent, hence dominating term, is M Div 32 ≡ −M 31 G 2 (X 0 ) y G2(X0) . As an illustration, consider the region 0 < ω < 6 and α σ > 0. From (24) and (14), we see that We can therefore write Plugging this expression in the divergent term of M 32 gives which must be negative. We arrive at the relation (74) in Section VI A Similarly, we find that • and α σ < 0 which are (75), (76) and (77) in Section VI A. Therefore, since G 2 (X) changes sign at X F P , in a neighborhood of a fixed point V F P , there are regions that satisfy (a) and regions that satisfy (b). As underlined in Section VI A, for V F P to be an attractor we require that When condition (a) is satisfied, the last term of the dδΩ DE /dN (see (B27)) is indeterminate in the limit. The last requirement for V F P to be an attractor is that, in linear perturbation theory, the following limit approaches dynamically its fixed-point limit given by (53). When condition (a) holds true we can expand (B42) to and substituting (B25) and (B26) we find that the limit is indeterminate.

Initial conditions for the numerical analysis
In order to numerically investigate whether a solution of the dynamical equations is approaching the fixed point we need to set the initial conditions for the dynamical variables in a neighborhood close to the fixed-point value. In this regard we perturb the dynamical variables around their fixed-point values instead of around their values close to the fixed point, V 0 . This is possible because any point in the neighborhood of a solution infinitesimally close to the fixed point can be treated as a perturbation around the fixed-point value. We set the initial values X i ,Ω DEi ,Ω mi ,Ω mi as follows: Ω ri = ∆Ω r (B47) We observe that, if we set y i = ∆y, then the term y G2(X) ≡ z(y, X) in (40) would not be close to its fixedpoint value given by (53). This happens because, on expanding z(y, X) around the fixed point, we obtain z(∆y, X ± + ∆X) = z(0, X ± ) (B48) The second and third terms in the above expansion have a divergent factor of 1 G2(X±) . Hence the initial value of r obtained from (40) is not close to its fixed-point value (52). We instead add a small perturbation to the fixed-point value of r, r i = 1 + ωG1(X) XG 1 (X)(1− ω /6) + ∆r and use this value in (24) to set the initial value of y as Appendix C: Minkowski limit During evolution of the dynamical variables, if we encounter H = 0, then the infinitesimal dN is zero, making the standard evolution equations indeterminate. In such a case, we revert to equations written in cosmic time. For notational convenience, we define T = tm g and the following dimensionless quantities: τ = σ M P l 2 , ρ m = ρm M 2 P l m 2 g ,ρ r = ρr M 2 P l m 2 g andh = a a , where the primes ( ) represent derivative with respect to T . The resulting background evolution equations are: where one must substitute r = τ α σ Notice that using T as the independent variable rather than N provides us with a system of equations with no divergences ash approaches zero, making it suitable to study the special fixed-point caseh = 0. However, these equations are inconvenient for the general stability analysis, since the stability matrix has fewer diagonal terms.

Fixed point solutions
We now focus on theh = 0 fixed-points analysis. One can easily verify from (C1), (C4) and (C5) that the derivatives of y,ρ m andρ r approach zero. These variables therefore approach constant values (and not necessarily 0) at the fixed point.
From (C2), at least one of X and τ must vanish. There are three cases: • If τ = 0 and X = 0, then (C3) is automatically satisfied and from the Friedmann equation (C7) From (14), we get n(t) 2 = (φ 0 ) 2 . From (24), we get y G2(X) 2 = 1, which makes r from (C6) indeterminate at the fixed point. r is also indeterminate from (C3). A special case occurs when G 1 (X) = 0. From (C8), we get ρ m = ρ r = 0. Since the initial values of matter and radiation density are not zero, we find that a → ∞.
We cannot find the eigenvalues of M analytically, because the characteristic polynomial is fifth order. Therefore we cannot apply the method described in Section VI A. We instead study the characteristic polynomial in the fixed-point limit. We find that in both fixed-point cases, τ = 0 and X = 0, the polynomial has indeterminate coefficients (due to the presence of ratios of powers ofh, τ or 1 − y G2(X) 2 ), or divergent coefficients (due to inverses of powers ofh, τ or 1 − y G2(X)

2
). This makes the eigenvalues and thus the stability of this special case indeterminate. Therefore, the stability of this fixed point can only be assessed by a numerical analysis. However, we consider such analysis as out of the scope of our present work. the kinetic coefficient obeys 0 < κ 2 V < 1. If m 2 GW > −H 2 , as imposed above, then this also ensures tachyonic stability for vector perturbations.
For simplicity, we will not evaluate the stability conditions for scalar perturbations in the region ω ≥ 6 for neither α σ > 0 nor α σ < 0.