# Failures of homogeneous and isotropic cosmologies in extended quasidilaton massive gravity

Stefano Anselmi, Saurabh Kumar, Diana López Nacir, Glenn D. Starkman

► To cite this version:

Stefano Anselmi, Saurabh Kumar, Diana López Nacir, Glenn D. Starkman. Failures of homogeneous and isotropic cosmologies in extended quasidilaton massive gravity. Physical Review D, 2017, 96 (8), pp.084001. 10.1103/PhysRevD.96.084001. hal-01645804

HAL Id: hal-01645804

https://hal.science/hal-01645804v1

Submitted on 22 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Failures of homogeneous and isotropic cosmologies in extended quasidilaton massive gravity

Stefano Anselmi,$^{1,2,3,*}$ Saurabh Kumar,$^{3,†}$ Diana López Nacir,$^{4,‡}$ and Glenn D. Starkman$^{3,§}$

$^{1}$LUTH, UMR 8102 CNRS, Observatoire de Paris, PSL Research University, Université Paris Diderot, 92190 Meudon, France

$^{2}$Institut d'Astrophysique de Paris, CNRS UMR 7095 and UPMC, 98bis, bd Arago, F-75014 Paris, France

$^{3}$Department of Physics/CERCA/Institute for the Science of Origins, Case Western Reserve University, Cleveland, Ohio 44106-7079, USA

$^{4}$Theoretical Physics Department, CERN, CH-1211 Genève 23, Switzerland
(Received 25 June 2017; published 3 October 2017)

We analyze the extended quasidilaton massive gravity model around a Friedmann-Lemaître-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 quasidilaton and the Hubble parameter. A less-well examined fixed-point solution, where the time derivative of the zeroth Stückelberg 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 that 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.

DOI: 10.1103/PhysRevD.96.084001

## I. INTRODUCTION

In the standard cosmological model, the accelerated expansion of the Universe is attributed to the cosmological constant $\Lambda$. However, to match the observed expansion, $\Lambda$ 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 quasidilaton 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 fixed-point cosmological solution of EQDMG has a de Sitter metric, and thus appears to be a good candidate for late-time cosmology.

EQDMG only differs from QDMG by the addition to the action of one operator involving the quasidilaton (QD) field and a new free parameter $\alpha_{\sigma}$. Naively, in the limit $\alpha_{\sigma} \to 0$, EQDMG reduces to QDMG, but actually the limit is very subtle. Indeed, there are controversial results in the literature regarding whether EQDMG contains an unavoidable additional degree of freedom, which would correspond to the Boulware-Deser ghost (BD) [9–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 nonstandard 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 that 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

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>* stefano.anselmi@iap.fr</footnote>

<footnote>† saurabh.kumar@case.edu</footnote>

<footnote>‡ diana.laura.lopez.nacir@cern.ch</footnote>

<footnote>§ glenn.starkman@case.edu</footnote>
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 Secs. II and III, we summarize the theory of EQDMG.

In Sec. IV, we define the dynamical variables and provide the relevant background equations.

In Sec. 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 Sec. 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 Sec. 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 Sec. VIII), the fixed points cannot be reached in the asymptotic future.

In Sec. 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. Case 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 Sec. IX we present our conclusions. We provide some detailed calculations and consider the special case of Minkowski solutions in the appendixes.

## II. FORMALISM

We consider the action for the extended quasidilaton theory [8],

$$S = S_{\text{EH}} + S_m + S_{\sigma} = \frac{M_{\text{Pl}}^2}{2} \int d^4 x \sqrt{-g} \left[ R - \frac{\omega}{M_{\text{Pl}}^2} \partial_{\mu} \sigma \partial^{\mu} \sigma + 2m_g^2 (\mathcal{L}_2 + \alpha_3 \mathcal{L}_3 + \alpha_4 \mathcal{L}_4) \right] + S_m, \quad (1)$$

where $M_{\text{Pl}}$ is the Planck mass and, in addition to the Einstein-Hilbert action $S_{\text{EH}}$, we have the contribution $S_m$ of the matter sector and a quasidilaton contribution $S_{\sigma}$. Here

$$\mathcal{L}_2 \equiv \frac{1}{2} ([\mathcal{K}]^2 - [\mathcal{K}^2]), \quad (2)$$

$$\mathcal{L}_3 \equiv \frac{1}{6} ([\mathcal{K}]^3 - 3[\mathcal{K}][\mathcal{K}^2] + 2[\mathcal{K}^3]), \quad (3)$$

$$\mathcal{L}_4 \equiv \frac{1}{24} ([\mathcal{K}]^4 - 6[\mathcal{K}]^2[\mathcal{K}^2] + 3[\mathcal{K}^2]^2 + 8[\mathcal{K}][\mathcal{K}^3] - 6[\mathcal{K}^4]), \quad (4)$$

with square brackets denoting a trace, and

$$\mathcal{K}_{\nu}^{\mu} = \delta_{\nu}^{\mu} - e^{\sigma/M_{\text{Pl}}} (\sqrt{g^{-1}} f)^{\mu}_{\nu}, \quad (5)$$

$$f_{\mu\nu} \equiv \eta_{ab} \partial_{\mu} \phi^a \partial_{\nu} \phi^b - \frac{\alpha_{\sigma}}{M_{\text{Pl}}^2 m_g^2} e^{-2\sigma/M_{\text{Pl}}} \partial_{\mu} \sigma \partial_{\nu} \sigma. \quad (6)$$

$S_{\sigma}$ includes five new fields: $\sigma$ is the quasidilaton scalar field and $\phi^a$ ($a = 0, \dots, 3$) are the four Stückelberg fields. It also depends on the coupling constants $\alpha_{\sigma}$, $\alpha_2$, and $\alpha_3$, and on the graviton mass $m_g$. For $\alpha_{\sigma} = 0$ one recovers the standard quasidilaton theory.

In the space of Stückelberg fields, the theory enjoys the Poincaré symmetry

$$\phi^a \rightarrow \phi^a + c^a, \quad \phi^a \rightarrow \Lambda_b^a \phi^b, \quad (7)$$

and a global symmetry given by

$$\sigma \rightarrow \sigma + \sigma_0, \quad \phi^a \rightarrow e^{-\sigma_0/M_{\text{Pl}}} \phi^a, \quad (8)$$

with $\sigma_0$ an arbitrary constant.

## III. BACKGROUND

We consider a spatially flat FLRW ansatz, for which

$$ds^2 = -N(t)^2 dt^2 + a(t)^2 \delta_{ij} dx^i dx^j, \quad (9)$$

$$\phi^0 = \phi^0(t), \quad (10)$$

$$\phi^i = x^i, \quad (11)$$

$$\sigma = \bar{\sigma}(t). \quad (12)$$

The fiducial metric $f_{\mu\nu}$ reduces to

$$f_{00} = -n(t)^2, \quad f_{ij} = \delta_{ij}, \quad (13)$$

where
$$n(t)^2 = (\dot{\phi}^0)^2 + \frac{\alpha_\sigma}{M_{\text{Pl}}^2 m_g^2} e^{-2\bar{\sigma}/M_{\text{Pl}}\dot{\bar{\sigma}}^2}. \quad (14)$$

The minisuperspace action for the background can be written as

$$S/V = M_{\text{Pl}}^2 \int dt \left[ -3 \frac{a^3}{N} \left(\frac{\dot{a}}{a}\right)^2 + a^3 \frac{w}{M_{\text{Pl}}} \frac{\dot{\sigma}^2}{2N} + N a^3 m_g^2 (\mathcal{L}_2 + \alpha_3 \mathcal{L}_3 + \alpha_4 \mathcal{L}_4) \right], \quad (15)$$

where

$$\mathcal{L}_2 = 3(X-1)(-2+X(1+r)), \quad (16)$$

$$\mathcal{L}_3 = -(X-1)^2(-4+X(1+3r)), \quad (17)$$

$$\mathcal{L}_4 = (X-1)^3(-1+rX), \quad (18)$$

and we have defined

$$X = \frac{e^{\bar{\sigma}/M_{\text{Pl}}}}{a}, \quad (19)$$

$$r = \frac{n}{N} a. \quad (20)$$

It is worth pointing out here that in (5),

$$(\sqrt{g^{-1}f})_{\nu}^{\mu} = \begin{bmatrix} \pm n/N & 0 & 0 & 0 \\ 0 & \pm 1/a & 0 & 0 \\ 0 & 0 & \pm 1/a & 0 \\ 0 & 0 & 0 & \pm 1/a \end{bmatrix},$$

and we make the (+, +, +, +) choice following [13]. This gives us

$$r > 0. \quad (21)$$

$r = 0$ represents a determinant singularity in either or both $f_{\mu\nu}$ (if $n = 0$) or $g_{\mu\nu}$ (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 \neq 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 $\phi^0(t)$ leads to the constraint equation

$$\partial_t \left[ \frac{\dot{\phi}^0}{n} a^4 G_2(X) \right] = 0, \quad (22)$$

where $G_2(X) = X(1-X)J(X)$, with

$$J(X) = 3 + 3(1-X)\alpha_3 + (1-X)^2\alpha_4. \quad (23)$$

This suggests that it will be useful to introduce as one of the dynamical variables

$$y(t) = \frac{\dot{\phi}^0(t)}{n} G_2(X). \quad (24)$$

The solution of (22) is immediately

$$y = \frac{C}{a^4}. \quad (25)$$

We see that, in any reasonable cosmological context, $y \to 0$ in the asymptotic future, and that this can be achieved by one (or more) of four quantities approaching or equaling zero: $J(X(t))$, $\dot{\phi}^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

$$\dot{X} = X \left( \frac{\dot{\sigma}}{M_{\text{Pl}}} - H \right). \quad (26)$$

$H(t) \equiv \dot{a}/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,

$$3H^2 = \frac{\omega}{2} \left( \frac{\dot{\sigma}}{M_{\text{Pl}}} \right)^2 + 3m_g^2 G_1(X) + \frac{\rho_m}{M_{\text{Pl}}^2} + \frac{\rho_r}{M_{\text{Pl}}^2}, \quad (27)$$

where

$$G_1(X) = \frac{1}{3}(X-1)(J(X) + (X-1)(\alpha_3(X-1) - 3)). \quad (28)$$

The form of the Friedmann equation suggests regarding the first two terms on the right-hand side (RHS) of (27) as the dark energy density (divided by $M_{\text{Pl}}^2$). For future convenience, we represent them separately as

$$\Omega_{\Lambda} \equiv \frac{m_g^2}{H^2} G_1(X) \quad (29)$$
and

$$ \Omega_{\sigma} \equiv \frac{\omega}{6H^2} \left( \frac{\dot{\sigma}}{M_{Pl}} \right)^2, \qquad (30) $$

and define $\Omega_{DE} \equiv \Omega_{\Lambda} + \Omega_{\sigma}$.

From the conservation of the stress-energy tensor of matter and of radiation we get

$$ \dot{\rho}_m = -3H\rho_m, \qquad (31) $$

$$ \dot{\rho}_r = -4H\rho_r. \qquad (32) $$

From the conservation of the stress-energy tensor<fnref n="1" /> obtained from $S_{\sigma}$, we get

$$ (\ddot{\sigma} + 3H\dot{\sigma})\omega\dot{\sigma} + 3M_{Pl}m_g^2(\dot{\sigma} - rHM_{Pl})XG_1'(X) = 0, \qquad (33) $$

where a prime here stands for a derivative with respect to the argument X of the function.

We select y, X, $\tilde{\Omega}_i = \Omega_i\tilde{h}^2$ ($i = DE, m, r$) as our dynamical variables, where $\tilde{h} \equiv H/m_g$ is the dimensionless Hubble parameter. We can now express the background evolution equations in terms of $N \equiv \log(a)$ (giving us $dN = Hdt$ assuming $H \neq 0$)<fnref n="2" />:

(i) From (22), we have immediately

$$ \frac{dy}{dN} = -4y. \qquad (34) $$

(ii) The equations for $\tilde{\Omega}_m$ and $\tilde{\Omega}_r$ are similarly easily obtained from (31) and (32),

$$ \frac{d\tilde{\Omega}_m}{dN} = -3\tilde{\Omega}_m, \qquad (35) $$

$$ \frac{d\tilde{\Omega}_r}{dN} = -4\tilde{\Omega}_r. \qquad (36) $$

(iii) Equation (26) can be rewritten using (29), (30), and the definition of $\tilde{\Omega}_{DE}$,

$$ \frac{dX}{dN} = X \left( \pm \sqrt{\frac{6(\tilde{\Omega}_{DE} - G_1(X))}{\tilde{h}^2\omega}} - 1 \right). \qquad (37) $$

The $\pm$ represents the possibility that $\dot{\sigma}$ can be positive or negative.<fnref n="3" />

(iv) Equation (33) can be rewritten using (29) and (30),

$$ \frac{d\tilde{\Omega}_{DE}}{dN} = -6(\tilde{\Omega}_{DE} - G_1(X)) + XG_1'(X)(r-1). \qquad (38) $$

In the above set of equations one must replace

(i) $\tilde{h}^2$ by the Friedmann equation, which now takes the simple form

$$ \tilde{h}^2 = \tilde{\Omega}_{DE} + \tilde{\Omega}_m + \tilde{\Omega}_r, \qquad (39) $$

(ii) and [combining (14), (20), and (25)] r with<fnref n="4" />

$$ r = +\sqrt{\frac{6(\tilde{\Omega}_{DE} - G_1(X))\alpha_{\sigma}}{\omega X^2 (1 - \frac{y}{G_2(X)})^2}}. \qquad (40) $$

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 \to 1$.

## 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 Eqs. (34)–(36) vanish, giving us $y_{FP} = 0$, and also $\tilde{\Omega}_{m,FP} = \tilde{\Omega}_{r,FP} = 0$. From (39), we learn that $\tilde{\Omega}_{DE,FP} = \tilde{h}_{FP}^2$.

The solutions to (34)–(36) are

$$ y = y_0 e^{-4N}, \qquad \tilde{\Omega}_m = \tilde{\Omega}_{m0} e^{-3N}, \qquad \tilde{\Omega}_r = \tilde{\Omega}_{r0} e^{-4N}, \qquad (41) $$

where $y_0$, $\tilde{\Omega}_{m0}$, and $\tilde{\Omega}_{r0}$ are the corresponding initial values. Thus fixed points occur in the asymptotic future, i.e., as $N \to \infty$ and so $a \to \infty$.

The left-hand side of (37) vanishes at the fixed point, implying that $X_{FP}$ is a constant. If $X_{FP} \neq 0$, the right-hand side of (37) [and if $X_{FP} = 0$, then the right-hand side of (38)] provides us with

<footnote n="1">¹Note that, using the constraint equation (22), one can show that the equation obtained by taking the variation of $S_{\sigma}$ with respect to $\sigma$ is not an independent equation.</footnote>

<footnote n="2">²The special case $H = 0$, i.e., Minkowski space, is discussed in Appendix C.</footnote>

<footnote n="3">³We will focus our attention below on the positive sign, because the negative sign leads to only an $X = 0$ fixed point.</footnote>

<footnote n="4">⁴For $\alpha_{\sigma} = 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 quasidilaton theory, and the evolution of the dynamical variables have been previously studied by the authors of [5,15].</footnote>
$$ \tilde{\Omega}_{\text{DE,FP}} = \tilde{h}_{\text{FP}}^2 = \begin{cases} \frac{G_1(X_{\text{FP}})}{(1-\frac{\omega}{6})}, & X_{\text{FP}} \neq 0, \\ G_1(0), & X_{\text{FP}} = 0, \end{cases} \quad (42) $$

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 < \omega < 6$ and $\omega \ge 6$ are in principle suitable regions in the parameter space, since $G_1(X_{\text{FP}})$ can be either positive or negative. A special value is $\omega = 6$, in which cases (37) and (42) demand $G_1(X_{\text{FP}}) = 0$.

Observing that $\tilde{h}^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:

$$ (\tilde{\Omega}_{\text{DE}}, \tilde{\Omega}_{m}, \tilde{\Omega}_{r}) = (\tilde{h}_{\text{FP}}^2, 0, 0), \quad (43) $$

where $\tilde{h}_{\text{FP}}^2$ 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:

(i) Case 1 (the standard case):

$$ J(X_{\text{FP}}) = 0, \quad (44) $$

and hence the fixed-point solutions are

$$ X_{\text{FP}} = X_{\pm} = 1 + \frac{3}{2} \frac{\alpha_3}{\alpha_4} \pm \sqrt{\frac{9 \alpha_3^2}{4 \alpha_4^2} - \frac{3}{\alpha_4}}. \quad (45) $$

(ii) Case 2:

$$ \left[ \frac{\dot{\phi}^0}{n} \right]_{\text{FP}} = \sqrt{1 - \alpha_{\sigma} \frac{6(\tilde{\Omega}_{\text{DE,FP}} - G_1(X_{\text{FP}}))}{\omega r_{\text{FP}}^2 X_{\text{FP}}^2}} = 0. \quad (46) $$

Since the left-hand side of (38) vanishes, (40) provides us with the following equation in $X_{\text{FP}}$:

$$ \begin{aligned} & 6(\tilde{h}_{\text{FP}}^2 - G_1(X_{\text{FP}})) + X_{\text{FP}} G_1'(X_{\text{FP}}) \\ & = G_1'(X_{\text{FP}}) \sqrt{\frac{6}{\omega} \alpha_{\sigma} (\tilde{h}_{\text{FP}}^2 - G_1(X_{\text{FP}}))}. \end{aligned} \quad (47) $$

Squaring (47) gives us a polynomial equation for $X_{\text{FP}}$. $X_{\text{FP}}$ can be any of the roots of that polynomial.

(iii) Case 3:

$$ X_{\text{FP}} = 0. \quad (48) $$

(iv) Case 4:

$$ X_{\text{FP}} = 1. \quad (49) $$

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(X_{\text{FP}}) = 0$

As will become clear below, this case is very subtle. Note that these fixed points are the same as the ones analyzed in [5,15] for the QD theory (i.e., the EQD with $\alpha_{\sigma} = 0$).

Requiring that $X_{\pm}$ be real means<fnref n="5" />

$$ \alpha_4 \le \frac{3\alpha_3^2}{4}. \quad (51) $$

From (38) and (42), we get the same expression as in the QD theory for the fixed-point limit of $r$ (assuming $\omega \ne 6$),

$$ r_{\text{FP}} = 1 + \frac{\omega G_1(X_{\text{FP}})}{X_{\text{FP}} G_1'(X_{\text{FP}}) (1 - \omega/6)}. \quad (52) $$

Note that this expression is valid for Cases 1, 2, and 4 (not $X = 0$), provided $\omega \ne 6$ and $G_1'(X) \ne 0$.

Unlike QD theory, since $\alpha_{\sigma} \ne 0$ we can use (24) to obtain (40), which gives $r$ in terms of the dynamical variables. If the system is to evolve toward its expected fixed point, $r$ must approach $r_{\text{FP}}$. Therefore, at the fixed point $z^2 = (\frac{\nu}{G_2(X)})^2$ should approach

$$ \begin{aligned} z_{\text{FP}}^2 &\equiv \left( \frac{y_{\text{FP}}}{G_2(X_{\text{FP}})} \right)^2 \\ &= 1 - \alpha_{\sigma} \frac{G_1(X_{\text{FP}})(1 - \omega/6)(G_1'(X_{\text{FP}}))^2}{[X_{\text{FP}} G_1'(X_{\text{FP}})(1 - \omega/6) + \omega G_1(X_{\text{FP}})]^2}. \end{aligned} \quad (53) $$

$$ (X_{\text{FP}}, r_{\text{FP}}, \omega) = \left( 1 + \frac{3}{\alpha_3}, 1 + \frac{6\tilde{h}^2}{X_{\text{FP}}^2}, 6 \right). \quad (50) $$

This solution for $\omega = 6$ and $\tilde{h}^2$ is indeterminate because of simplifications occurring in (37). Therefore, on the $\alpha_4 = \frac{2\alpha_3^2}{3}$ hypersurface in the parameter space, the dynamical equations lose predictive power.

<footnote n="5">⁵A special case occurs when $\alpha_4 = \frac{2\alpha_3^2}{3}$, in which case $G_1(X) = 0 = J(X)$. We are left with</footnote>
For this approach to be smooth, one also needs

$$ \frac{dz}{dN} = - \left[ 4 + \frac{G_2'(X)}{G_2(X)} \frac{dX}{dN} \right] z \to 0. \quad (54) $$

One possibility is that $z \to 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.<fnref n="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 fixed-point 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 $\alpha_\sigma \neq 0$.

### B. Case 2: $[\dot{\phi}_n^0]_{FP} = 0$

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_{FP}$, the sign of $6(\tilde{h}_{FP}^2 - G_1(X_{FP})) + X_{FP}G_1'(X_{FP})$ should be the same as the sign of $G_1'(X_{FP})$.<fnref n="7" />

It is worthwhile noting that for Case 2, $\alpha_\sigma > 0$. (We omit the QD case, $\alpha_\sigma = 0$.) This can be seen by inspection of Eq. (46), recalling that

$$ \tilde{\Omega}_{\text{DE}} - G_1(X) = \frac{\omega}{6} \left( \frac{\dot{\sigma}}{m_g M_{\text{Pl}}} \right)^2 \ge 0. $$

### C. Case 3: $X_{FP} = 0$

From (42) and the definition of $\tilde{\Omega}_{\text{DE}}$, (29) and (30), we get

$$ \dot{\sigma}_{\text{FP}} = 0. \quad (55) $$

We find that $r_{FP}$ is indeterminate from both (38) and (40).

Examining (19), we see that there are two ways to get $X_{FP} = 0$. The first possibility is that $\sigma \to -\infty$, in which case we can draw no conclusion from (14) about the value of $\dot{\phi}_n^0$ at the fixed point. The second possibility is that $n^2 = (\dot{\phi}_n^0)^2$ at the fixed point.

### D. Case 4: $X_{FP} = 1$

Since $G_1(1) = 0$ and $G_1'(1) = 1$, and since $\tilde{\Omega}_{\text{m}} = \tilde{\Omega}_r = 0$ at the fixed point, from (37) and (39), we find that $X_{FP} = 1$ requires $\omega = 6$. Substituting $X = X_{FP} = 1$ in (38), we get

$$ r_{\text{FP}} = 1 + 6\tilde{h}_{\text{FP}}^2, \quad (56) $$

where $\tilde{h}_{FP}$ 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 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 $\mathbf{V} = [y, X, \tilde{\Omega}_{\text{DE}}, \tilde{\Omega}_{\text{m}}, \tilde{\Omega}_r]^T$ denote the dynamical variables and $\mathbf{f}(\mathbf{V})$ be the RHS of the first order differential equations. Thus we can express Eqs. (34)–(38) as

$$ \frac{d\mathbf{V}}{dN} = \mathbf{f}(\mathbf{V}). \quad (57) $$

Assuming small perturbations $\delta\mathbf{V}$ around any point, $\mathbf{V}_0$, a Taylor expansion of the functions $\mathbf{f}(\mathbf{V})$ gives us

$$ \frac{d}{dN} \delta \mathbf{V} = \mathbf{M} \delta \mathbf{V} + \mathbf{f}(\mathbf{V}_0), \quad (58) $$

where $\mathbf{M}$ is the stability matrix. Its elements are given by

$$ M_{ij} = \left[ \frac{\partial f_i(\mathbf{V})}{\partial V_j} \right]_{\mathbf{V}=\mathbf{V}_0}, \quad i, j = 1, \dots, 5. \quad (59) $$

In Appendix A, we provide the analytical expressions for the elements of $\mathbf{M}$ [Eqs. (A1)–(A7)].

If $\mathbf{V}_0$ is a fixed point, then the second term of RHS of (58) would vanish. Using the eigenvectors of $\mathbf{M}$, one can then find matrix $\mathbf{P}$ such that

$$ \mathcal{D} = \mathbf{P}^{-1} \mathbf{M} \mathbf{P}, \quad (60) $$

$$ = \text{Diag}(\lambda_1, \dots, \lambda_5), \quad (61) $$

where $(\lambda_1, \dots, \lambda_5)$ are the eigenvalues given by (A16). We define $\delta\tilde{\mathbf{V}}$ as

$$ \mathbf{P} \delta \tilde{\mathbf{V}} \equiv \delta \mathbf{V}. \quad (62) $$

Thus the solution to (58) in the new basis would be

$$ \delta \tilde{V}_i = e^{\lambda_i N} C_i, \quad (63) $$

<footnote n="6">⁶This constraint is Eq. (7) of [16].</footnote>

<footnote n="7">⁷There is a special solution of (47) where $G_1(X) = G_1'(X) = 0$. This is possible only when $3+2\alpha_3+3\alpha_3^2-4\alpha_4=0$ and the common root of $G_1(X)$ and $G_1'(X)$ is $X_{FP} = (3+5\alpha_3+2\alpha_4)/(2\alpha_3+2\alpha_4)$, which corresponds to $\tilde{h}_{FP} = 0$. Thus we return to the special case $h=0$ discussed in Appendix C.</footnote>
where $C_i$'s are integration constants. Multiplying the above equation by $\mathbf{P}$ and thereby returning to the original basis, we get

$$\delta V_i = \sum_{j=1}^{5} P_{ij} e^{\lambda_j N} C_j. \qquad (64)$$

Then, for the fixed points to be stable, we require $\delta\mathbf{V}$ to approach zero as $N \to \infty$. It can be seen from (64) that if the eigenvalues of $\mathbf{M}$ are either real and negative or imaginary with a 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, $\lambda_4$ and $\lambda_5$ must be real and negative or complex with negative real parts. This requires the elements of $\mathbf{M}$ to satisfy the condition

$$0 \le (3 + 2M_{22})^2 + 4M_{23}M_{32} < 9 \qquad (65)$$

for stable solutions or

$$(3 + 2M_{22})^2 + 4M_{23}M_{32} < 0 \qquad (66)$$

for stable spiral solutions.

### A. Case 1: $J(X_{FP}) = 0$

The above-described standard method will not suffice to evaluate fixed-point stability in all cases. In particular, for $J(X_{FP}) = 0$ the stability matrix $\mathbf{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 $\mathbf{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 $\delta\mathbf{V}$ around an arbitrary point $\mathbf{V}_0$; using (58) and diagonalizing $\mathbf{M}$, we solve for $\delta\mathbf{V}$. For the fixed point to be an attractor we require the following conditions to be satisfied. If $\mathbf{V}_0$ is infinitesimally close to the fixed point $\mathbf{V}_{FP}$, then, as $N \to \infty$:

(A) the perturbations $\delta\mathbf{V}$ are infinitesimally close to zero, and therefore we require

$$\lim_{\substack{N \to \infty \\ V_0 \to V_{FP}}} \delta \mathbf{V} \to 0; \qquad (67)$$

(B) the derivatives of perturbations, $d\delta\mathbf{V}/dN$, are infinitesimally close to zero, and therefore we must verify that

$$\lim_{\substack{N \to \infty \\ V_0 \to V_{FP}}} \frac{d}{dN} \delta \mathbf{V} \to 0; \qquad (68)$$

(C) $y/G_2(X)$ is infinitesimally close to its fixed-point limit (53)— in compact form this translates to

$$\lim_{\substack{N \to \infty \\ V_0 \to V_{FP}}} \left(\frac{y}{G_2(X)}\right)^2 \to 1 - \alpha_{\sigma} \frac{G_1(X_{FP})(1 - \omega/6)(G_1'(X_{FP}))^2}{[X_{FP}G_1'(X_{FP})(1 - \omega/6) + \omega G_1(X_{FP})]^2}. \qquad (69)$$

In more detail, after diagonalizing the matrix $\mathbf{M}$, from (58) we obtain

$$\frac{d}{dN} \delta \bar{\mathbf{V}} = \mathcal{D} \delta \bar{\mathbf{V}} + \mathbf{B}, \qquad (70)$$

where

$$\mathbf{B} = P^{-1} \mathbf{f}(\mathbf{V}_0), \qquad (71)$$

Solving (70), we are left with

$$\delta \bar{\mathbf{V}} = \mathrm{Diag} \left( -\frac{1}{\lambda_1}, \dots, -\frac{1}{\lambda_5} \right) \mathbf{B} + \mathrm{Diag} (e^{\lambda_1 N}, \dots, e^{\lambda_5 N}) \mathbf{C}, \qquad (72)$$

where $\mathbf{C}$ is the integration constant vector.

Upon returning to the original basis, we find

$$\delta \mathbf{V} = \mathbf{P} \mathrm{Diag} \left( -\frac{1}{\lambda_1}, \dots, -\frac{1}{\lambda_5} \right) \mathbf{P}^{-1} \mathbf{f}(\mathbf{V}_0) + \mathbf{P} \mathrm{Diag} (e^{\lambda_1 N}, \dots, e^{\lambda_5 N}) \mathbf{C}. \qquad (73)$$

In Appendix B we show that the requirement (A) could be satisfied if the eigenvalues $\lambda_i$'s are either real and negative or complex with a negative real part. It is also shown that if $X_0$ is a point infinitesimally close to its fixed-point value $X_{FP}$, it must satisfy the following relations$^8$:

(1) $\alpha_{\sigma} > 0$

$$G'_1(X_0)G_2(X_0)G'_2(X_0) > 0, \quad 0 < \omega < 6, \qquad (74)$$

$$G'_1(X_0)G_2(X_0)G'_2(X_0) < 0, \quad \omega > 6, \qquad (75)$$

(ii) $\alpha_{\sigma} < 0$

$$G'_1(X_0)G_2(X_0)G'_2(X_0) < 0, \quad 0 < \omega < 6, \qquad (76)$$

$$G'_1(X_0)G_2(X_0)G'_2(X_0) > 0, \quad \omega > 6. \qquad (77)$$

<footnote>$^8$In special cases where $G_1(X_{FP}) = 0$, such as $\omega = 6$, or $2\alpha_3^2 = \alpha_4$, the eigenvalues involve ratios of zeros that we are unable to resolve, so we cannot determine the stability.</footnote>
In the numerical investigation performed in Sec. 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_{FP}$ and never both.<fnref n="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) should be performed. In practice the analysis of $(\frac{y}{G_2(X)})^2$ is enough. Indeed, we find that, even though the dynamical variables approach their fixed-point values in the asymptotic future, $(\frac{y}{G_2(X)})^2$ oscillates with maxima that grow in time. Eventually, it reaches the critical value of 1, making $r$ in Eq. (40) singular. We discuss this phenomenon in greater detail in Sec. VII.

### B. Case 2: $[\frac{\partial^0}{n}]_{FP} = 0$

In this case, all elements of the stability matrix are well defined and finite in the fixed-point limit. Hence we use the standard approach of fixed-point analysis. The analytical form of the elements of $M$ around the fixed-point solutions are provided in Appendix A [Eqs. (A9)–(A15)]. For the fixed points to be attractors, the elements of $M$ must satisfy the conditions (65) for stable solutions and (66) for stable spiral solutions.

### C. Case 3: $X_{FP} = 0$

Recalling (55), we find that the presence of terms $\frac{X}{\Omega_\sigma}$ and $\frac{y}{G_2(X)}$ in the stability matrix make the elements (A1), (A2), (A5), (A6), and (A7) and consequently the eigenvalues (A16) indeterminate. Hence the stability of $X_{FP} = 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 $(\alpha_3, \alpha_4, \omega, \alpha_\sigma)$, we check numerically that the five 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 Sec. VI to guarantee linear stability. We look further for values of the parameters for which

(i) $X_{FP} \ge 0$;

(ii) $\tilde{h}_{FP}^2$ given by (42) is positive;

(iii) the fixed-point linear stability conditions [given by (74)–(77) for Case 1, and (65) or (66) for Case 2] are satisfied.

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(X_{FP}) = 0$

Our goal in this section is nominally to identify values of the parameters $(\alpha_3, \alpha_4, \omega, \alpha_\sigma)$ 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 begin the discussion with an example by setting $(\alpha_3, \alpha_4, \omega, \alpha_\sigma) = (3.1, 3.1, 5.5, 10.0)$. Solving (45), we find that $X_+$ satisfies the conditions (i) and (ii). This set of parameter values also satisfies the tensor, vector, and scalar stability at the $X_+$ fixed point (see Appendix E).

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 < \omega < 6$ and $\alpha_\sigma > 0$, the initial value of $X$ must be of the form $X_0 = X_+ + \epsilon$, where $0 < \epsilon \ll 1$ so that (74) is satisfied. To study the behavior of the dynamical variables close to the fixed point, we set $\epsilon = 10^{-6}$ and the initial conditions to be $\tilde{\Omega}_{DE,0} = \tilde{\Omega}_{DE,FP} + \epsilon$, $\tilde{\Omega}_{m,0} = \epsilon$, $\tilde{\Omega}_{r,0} = \epsilon$. 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 $\epsilon$ and initialize $y$ using (B49). For a detailed discussion, we refer the reader to Appendix B. Notice that $\tilde{\Omega}_{DE,0}$, $r_0$, $y_0$

<footnote n="9">⁹We exclude particular cases for which in a neighborhood of $X_{FP}$ not only $G_2(X)$ changes sign but also either $G'_2(X)$ or $G'_1(X)$.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. X(N) (top panel) and $\tilde{\Omega}_{DE}(N)$ (bottom panel) starting from initial values of dynamical variables $\epsilon = 10^{-6}$ away from their fixed-point values. The parameters are set to $(\alpha_3, \alpha_4, \omega, \alpha_\sigma) = (3.1, 3.1, 5.5, 10.0)$. X and $\tilde{\Omega}_{DE}$ approach their fixed-point values until they hit singularities.</figcaption></figure>

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,

$$y = y_0 e^{-4N}, \quad \tilde{\Omega}_m = \tilde{\Omega}_{m,0} e^{-3N}, \quad \tilde{\Omega}_r = \tilde{\Omega}_{r,0} e^{-4N}. \qquad (78)$$

The evolution equations for X and $\tilde{\Omega}_{DE}$ have no analytical solutions. Therefore, we use (78) in (39) and (40), and we

solve Eqs. (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 $\tilde{\Omega}_{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) we can identify the possible sources of the singular behavior. The only possibilities are that, as the fixed point is approached, the square-root term $(\tilde{\Omega}_{DE} - G_1(X))$ in (37) or $(1 - (\frac{y}{G_2(X)})^2)$ in (38) approaches zero. From the top panel of Fig. 3 we see that $(\tilde{\Omega}_{DE} - G_1(X))$ does not vanish as the fixed point is approached. On the other hand in Fig. 2 we show that $z^2 = (\frac{y}{G_2(X)})^2$ exhibits an oscillating behavior where the peaks grow in time monotonically. Thus $z$ does not approach its expected fixed-point value. Moreover, as shown in the bottom panel of Fig. 3, it eventually approaches unity, causing a breakdown in the coupled dynamical equations for X and $\tilde{\Omega}_{DE}$.

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 \to 1$ behavior was robust to increase in the numerical precision demanded from the integrator, increase in the order of the integrator, change in the integration scheme, and change in the initial conditions. Therefore, we have conclusive evidence that the evolution of the dynamical variables drives $z^2$ toward 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.

<figure><img src="image_3.png" /><figcaption>FIG. 2. $z^2(N)$ evaluated from solving the full equations (non-perturbative) and linear perturbation theory. Parameter values and initial conditions are the same as for Fig. 1.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. $(\tilde{\Omega}_{DE} - G_1(X))$ (top panel) and $z^2$ (bottom panel) in function of $N$. $(\tilde{\Omega}_{DE} - G_1(X))$ is well behaved. $z^2$ shows oscillations with peaks that grow in time and eventually reach unity making the equations singular.</figcaption></figure>

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 $(\alpha_3, \alpha_4, \omega, \alpha_\sigma)$ 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 nonperturbative 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 large— outside 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 Eqs. (B25), (B26), (B27), (B28), and (B29) where the $C_i$ coefficients are determined by requiring that $\delta V(N_{in} \equiv 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} < |\alpha_3|, |\alpha_4|, \alpha_\sigma < 10^5, \\ 0 < \omega < 6 \quad \text{or} \quad 6 < \omega < 10^5, \qquad (79)$$

and evaluate $X_\pm$. 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 $\epsilon = 10^{-8}$ to $\epsilon = 10^{-12}$ away from their fixed-point values.<fnref n="10" /> We repeated the same procedure using the numerical solutions to the full (nonperturbative) equations with a few thousand randomly selected points. We find that the peaks of $z^2$ always grow monotonically both in linear perturbation theory and by using the exact equations.

Notice that for $\alpha_\sigma < 0$, Eqs. (14) and (24) imply $z^2 > 1$. However, even for negative values of $\alpha_\sigma$, $r$ from (40) still becomes singular because the troughs of $z^2$ monotonically decrease and approach unity.

<footnote n="10">¹⁰For simplicity we only consider $\tilde{\Omega}_{DE,0}$, $r_0$, $y_0$ larger than their fixed-point values. We numerically inspect a few examples with initial conditions smaller than the fixed-point values, and we find the same singular behavior.</footnote>
Hence we come to the following conclusion: the standard (Case 1) fixed-point solutions of extended quasidilaton 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.

### B. Case 2: $[\dot{\phi}/n]_{FP} = 0$

Recalling that $\alpha_\sigma$ 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 seventh 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<fnref n="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 ($\epsilon = 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 (Sec. VI), we find that the dynamical variables reach their fixed-point values in the asymptotic future.

### C. Case 3: $X_{FP} = 0$

For the Case 3 fixed points we do not have indications from the linear stability analysis. We perform a random search selecting $64 \times 10^4$ points in the parameter space. Starting from small perturbations ($\epsilon = 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 $(\tilde{\Omega}_{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.

Because of disagreements in the literature mentioned above (see [9–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 $\dot{\sigma} = 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 $\dot{\sigma} = 0$ corresponds to a Minkowski background metric (provided $X \neq 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 \forall N > 0$. This possibility can also be ruled out because $X = 0$ only makes sense as a fixed-point limit. Hence from now on, we will not consider Case 3 any further.

Following the standard treatment in this section only, we take the action of the matter sector to be<fnref n="12" />

$$S_m = \int d^4 x \sqrt{-g} P(Y, \chi), \qquad (80)$$

which corresponds to the addition of a scalar field $\chi$ with a noncanonical kinetic term, given by the function $P(Y, \chi)$, where

$$Y \equiv -\frac{1}{2} \partial_{\mu} \chi \partial^{\mu} \chi. \qquad (81)$$

The fluid variables (pressure $p$, energy density $\rho$, and sound speed $c_s$) associated with $\chi$ can be written as [18]

$$p = P(Y, \chi), \quad \rho \equiv 2P_{,Y}Y - P, \quad c_s^2 \equiv \frac{P_{,Y}}{\rho_{,Y}}. \qquad (82)$$

To study the perturbations in $\chi$, we replace $\chi$ by $\chi + M_{Pl}\delta\chi$.

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

$$\begin{align}
\delta g_{00} &= -2\Phi, \nonumber \\
\delta g_{0i} &= a(\partial_i B + B_i), \nonumber \\
\delta g_{ij} &= a^2 \left[ 2\delta_{ij}\psi + \left( \partial_i \partial_j - \frac{\delta_{ij}}{3} \partial^k \partial_k \right) E + \partial_{(i} E_{j)} + h_{ij} \right], \tag{83}
\end{align}$$

where $\delta^{ij}h_{ij} = \partial^i h_{ij} = \partial^i E_i = \partial^i B_i = 0$ and Latin indices are raised with $\delta^{ij}$. For the quasidilaton field we write the background plus perturbation by replacing $\sigma$ by $\sigma + M_{Pl}\delta\sigma$.

<footnote n="11">¹¹We keep this analysis agnostic to the stability conditions for scalar perturbations since, as we show below, the fact that $\dot{\phi}^0 = 0$ does not guarantee the absence of the BD ghost.</footnote>

<footnote n="12">¹²Clearly, the action (80) does not describe dustlike matter. However, we consider it only for simplicity because it will be enough for the purposes of this section.</footnote>
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 we relegate the remaining details to Appendix D.

We focus first on the case $\dot{\phi}^0 \neq 0$ and, following literature, we adopt the unitary gauge $\phi^a = \phi^0(t)\delta_0^a + \delta_i^a x^i$. The case $\dot{\phi}^0(t) = 0$, for which we cannot work in this gauge, is analyzed later.

$$A. \dot{\phi}^0 \neq 0$$

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 $\Phi$ and $B$, we obtain constraint equations that allow us to eliminate these two variables. The solution for $\Phi$ and $B$ can be found for a generic background and without making use of the background equations.<fnref n="13" /> Introducing the solutions for $\Phi$ and $B$ back into the action, we can write the kinetic part as

$$S_{\text{scalar}}^{(2)} \supset \frac{M_{\text{Pl}}^2}{2} \int d^3k dta^3 \dot{Z}^\dagger Q \dot{Z}, \quad (84)$$

with $Z = \{\psi, \delta\sigma, E, \delta\chi\}$. One combination corresponds to the Boulware-Deser mode. For this mode to be non-dynamical, the determinant of the matrix $Q$ must vanish.

We first consider the vacuum case, and then we see how the inclusion of matter changes the results:

(i) Vacuum case: In the absence of the additional field $\chi$, $Q$ is a $3 \times 3$ matrix and $Z = \{\psi, \delta\sigma, E\}$. The determinant $\text{Det}(Q)$ can be computed analytically in Fourier space. After expanding in powers of comoving wave number $k$, keeping only the leading order terms (the infrared part), and using the background equations to express $H$ and $\dot{H}$ in terms of the other dynamical variables,

$$\text{Det}(Q) = \frac{4\omega\alpha_{\sigma}M_{\text{Pl}}^4 a^2 J(X) \dot{\phi}^{02} \dot{\sigma}^2 k^4}{r^3 X m_g^2 (J(X) + \alpha_3 (X-1)^2 - 3X + 3)} + \mathcal{O}(k^6). \quad (85)$$

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 $\text{Det}(Q)$ vanishes when $J(X) = 0$ or $\dot{\sigma} = 0$. Moreover, in that case it can be shown that the determinant vanishes at all orders in $k$. As we

show next, these results are robust to the addition of matter.

(ii) With matter: Proceeding as above, at leading order in the wave number $k$, after using the background equations to replace $\dot{H}$ and $H$ in terms of the other dynamical variables, the determinant of the now $4 \times 4$ matrix $Q$ can be written as

$$\begin{aligned} \text{Det}(Q) = & \frac{32\omega}{r^3 X} \alpha_{\sigma} M_{\text{Pl}}^6 a^2 J(X) \dot{\phi}^{02} (X-1) P_{,Y} (Y, \chi) \dot{\sigma}^2 k^4 \\ & \times \{p + \rho(1 - 2c_s^2) - 2c_s^2(X-1)m_g^2 M_{\text{Pl}}^2 \\ & \times (J(X) + \alpha_3 (X-1)^2 - 3X + 3)\}^{-1} \\ & + \mathcal{O}(k^6). \end{aligned} \quad (86)$$

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 that the determinant vanishes at all orders in $k$.

Of course, the determinant also vanishes in the case where the matter is just a cosmological constant, $P_{,Y} = 0$.

Therefore, $\text{Det}(Q) = 0$ at the Case 1 fixed points. If we perturb $X$ away from that fixed-point value, then $J(X) \neq 0$, and $\text{Det}(Q) \neq 0$ either. However, from (22), we can conclude that if $J \equiv J(X) = 0$ at some initial time $t_0$, then also $\dot{J} = 0$ at that time, and also $\dot{J}$. Thus $J(X)$ remains zero at all times once $X = X_{\text{FP}}$. Consequently $\text{Det}(Q)$ also remains equal to zero. Therefore, to eliminate the BD ghost in both the vacuum and the matter contexts, we must impose special initial conditions on the dynamical variables, namely, that $X = X_{\pm}$ is exactly satisfied. Setting $X$ to this fixed-point value, and thus setting $J(X) = 0, y = 0$, and $\dot{X} = 0$, we can solve (35) through (40) to obtain

$$\begin{aligned} \tilde{\Omega}_{\text{DE}} &= \frac{\frac{\omega}{6} (\tilde{\Omega}_{m0} e^{-3N} + \tilde{\Omega}_{r0} e^{-4N}) + G_1(X_{\pm})}{1 - \frac{\omega}{6}}, \\ \tilde{h}^2 &= \frac{6}{\omega} (\tilde{\Omega}_{\text{DE}} - G_1(X_{\pm})). \end{aligned} \quad (87)$$

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 quasidilaton 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 \to -\infty$, the matter and radiation terms will dominate over $G_1(X_{\pm})$. This restricts the values of $\omega$ to $0 < \omega < 6$ or it means that the theory is not valid arbitrarily far into the past.

<footnote n="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.</footnote>
$$B. \dot{\phi}^0 = 0$$

As mentioned earlier, in the special case $\dot{\phi}^0 = 0$ (Case 2) we cannot use the unitary gauge. Assuming $H \neq 0$ we choose the gauge with $\psi = 0$ instead, while we keep $\phi^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 $\dot{\phi}^0(t_0) = 0$ at some initial time $t_0$, yielding $\ddot{\phi}^0(t_0) = 0$ and $\dot{\phi}^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 $\delta\phi^0$ (since $\phi^0 = 0$ at the background level) only enter as a contribution to $f_{\mu\nu}$ that is quadratic in $\delta\phi^0$. Therefore, the quadratic action for this scalar decouples from the other parts, and the kinetic part can be immediately computed,

$$S_{\text{scalar}}^{(2)} \supset -\frac{M_{\text{Pl}}^2}{2} \int d^3k dt a^3 \frac{m_g^2}{r} G_2(X) a^2 |\delta\dot{\phi}^0|^2. \quad (88)$$

Now, to compute the determinant of the kinetic matrix corresponding to the other degrees of freedom we proceed as above. We integrate out $\Phi$ 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:

(i) Vacuum case: In this case $Z = \{\delta\sigma, E\}$, and

$$\begin{aligned} \text{Det}(Q) = & G'_1(X)k^4m_g^2X\omega a^2(6\alpha_\sigma G_1(X) + r^2X^2\omega) \\ & \times \{4(9\alpha_\sigma G_1(X)G'_1(X)m_g^2Xa^2 \\ & + k^2(r+1)(6\alpha_\sigma G_1(X) \\ & + r^2X^2\omega))\}^{-1}; \end{aligned} \quad (89)$$

(ii) With matter: $Z = \{\delta\sigma, E, \delta\chi\}$, and

$$\begin{aligned} \text{Det}(Q) = & -H^2k^4X\omega m_g^2 G'_1(X)M_{\text{Pl}}^2 \\ & \times [Xm_g^2(r^2X\omega - (r-1)\alpha_\sigma G'_1(X)) + 2\dot{H}\alpha_\sigma] \\ & \times \{\dot{\chi}^2[a^2Xm_g^2G'_1(X)(2\alpha_\sigma(3H^2c_s^2 + \dot{H}) \\ & - r^2X^2\omega(c_s^2 - 1)m_g^2) \\ & - a^2(r-1)X^2\alpha_\sigma m_g^4(G'_1(X))^2 \\ & + 4H^2k^2(r+1)\alpha_\sigma c_s^2]\}^{-1}. \end{aligned} \quad (90)$$

Hence, we conclude that the condition $\dot{\phi}^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 toward their fixed-point 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 quasidilaton 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–22]).

As mentioned below (87), the condition $\omega < 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 $\omega < 6$ is in principle unnecessary. It would be worth exploring the parameter region $\omega > 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.

## ACKNOWLEDGMENTS

We thank especially C. de Rham for helpful insights and M. Fasiello, K. Hinterbichler, A. Tolley, and S.-Y. Zhou for discussions. S. A., S. K., and G. D. S. are partially supported by CWRU Grant No. DOE-SC0009946. The research leading to these results has received funding from the European Research Council under the European Community Seventh Framework Programme
(FP7/2007-2013 Grant Agreement No. 279954) ERC-StG “EDECS.” S. A. and G. D. S. thank the CERN theory group for their hospitality. We acknowledge the use of the xAct-xPand package for Mathematica [23].

## APPENDIX A: STABILITY MATRIX

This appendix complements the results presented in Sec. VI. In particular, we provide, in linear perturbation theory, the full expressions for the stability matrix and its eigenvalues. We then present the Case 2 fixed-point limit of the matrix elements and the conditions of stability.

### 1. General expression

The derivative/stability matrix **M** defined in (59) takes the form

$$ M = \begin{bmatrix} -4 & 0 & 0 & 0 & 0 \\ 0 & M_{22} & M_{23} & M_{24} & M_{25} \\ M_{31} & M_{32} & M_{33} & 0 & 0 \\ 0 & 0 & 0 & -3 & 0 \\ 0 & 0 & 0 & 0 & -4 \end{bmatrix}, $$

where

$$ M_{22} = -1 \pm \sqrt{\frac{6\tilde{\Omega}_{\sigma}}{\tilde{h}^2\omega}} \left(1 - \frac{XG'_{1}(X)}{2\tilde{\Omega}_{\sigma}}\right), \quad (A1) $$

$$ M_{23} = \pm \frac{\sqrt{\frac{3}{2}X(\tilde{h}^2 - \tilde{\Omega}_\sigma)}}{\sqrt{\tilde{\Omega}_\sigma \tilde{h}^6 \omega}}, \quad (A2) $$

$$ M_{24} = \mp X \sqrt{\frac{3\tilde{\Omega}_{\sigma}}{2\tilde{h}^6 \omega}}, \quad (A3) $$

$$ M_{25} = M_{24}, \quad (A4) $$

$$ M_{31} = \frac{\sqrt{\frac{6\alpha_{\sigma}\tilde{\Omega}_{\sigma}}{\omega(1-\left(\frac{y}{G_2(X)}\right)^2)}} G'_{1}(X)}{y\left(1-\left(\frac{y}{G_2(X)}\right)^2\right)} \left(\frac{y}{G_2(X)}\right)^2, \quad (A5) $$

$$ M_{32} = -\sqrt{\frac{3}{2}} \sqrt{\frac{\alpha_{\sigma}}{\omega(1 - \left(\frac{y}{G_2(X)}\right)^2)\tilde{\Omega}_{\sigma}}} (G'_{1}(X))^2 \\ + \left(5G'_{1}(X) - M_{31}G'_{2}(X)\left(\frac{y}{G_2(X)}\right)\right) \\ + XG''_{1}(X) \left(-1 + \sqrt{\frac{6\alpha_{\sigma}\tilde{\Omega}_{\sigma}}{X^2\omega(1 - \left(\frac{y}{G_2(X)}\right)^2)}}\right), \quad (A6) $$

$$ M_{33} = -6 + \sqrt{\frac{3}{2}} \sqrt{\frac{\alpha_{\sigma}}{\omega(1 - \left(\frac{y}{G_2(X)}\right)^2)\tilde{\Omega}_{\sigma}}} G'_{1}(X). \quad (A7) $$

In the above equations, one must replace $\tilde{h}^2$ by (39) and $\tilde{\Omega}_\sigma$ by $\tilde{\Omega}_{\text{DE}} - G_1(X)$.

### 2. Case 2

We evaluate **M** for the Case 2 fixed point. Notice that $\tilde{\Omega}_\sigma$, when $X \neq 0$ [using (26)], is given by

$$ \tilde{\Omega}_{\sigma} = \frac{\omega \tilde{h}^2}{6}. \quad (A8) $$

As explained in Sec. 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

$$ M_{22} = -\frac{3XG'_{1}(X)}{\tilde{h}^2\omega}, \quad (A9) $$

$$ M_{23} = \frac{3X}{\tilde{h}^2\omega} (1 - \omega/6), \quad (A10) $$

$$ M_{24} = -\frac{X}{2\tilde{h}^2}, \quad (A11) $$

$$ M_{25} = M_{24}, \quad (A12) $$

$$ M_{31} = 0, \quad (A13) $$

$$ M_{32} = -3 \sqrt{\frac{\alpha_{\sigma}}{\omega^2 \tilde{h}^2}} (G'_{1}(X))^2 + 5G'_{1}(X) \\ + XG''_{1}(X) \left(-1 + \sqrt{\frac{\alpha_{\sigma} \tilde{h}^2}{X^2}}\right), \quad (A14) $$

$$ M_{33} = -6 + 3 \sqrt{\frac{\alpha_{\sigma}}{\omega^2 \tilde{h}^2}} G'_{1}(X), \quad (A15) $$

at the fixed point. In the above equations, one must$^{14}$ use (42) to replace $\tilde{h}^2$ by $\frac{G_1(X)}{(1-\omega/6)}$.

The eigenvalues ($\lambda_i$'s) of the matrix are
$$ \lambda_i = -4, -4, -3, $$

$$ \times \frac{1}{2} \left( (M_{22} + M_{33}) \mp \sqrt{(M_{22} - M_{33})^2 + 4M_{23}M_{32}} \right), \\ i = 1, \dots, 5, \quad (A16) $$

$$ (6\tilde{h}^2 + XG'_{1}(X))^2 = (G'_{1}(X))^2\alpha_{\sigma}\tilde{h}^2. $$

<footnote>$^{14}$Except in the particular case where $\omega = 6$, in which case one must instead use a root of</footnote>
respectively. Using (47) we find that

$$M_{22} + M_{33} = -3. \qquad (A17)$$

We can thus impose the following conditions on elements (A9)–(A15) so that the fixed points become stable or stable spirals:

(i) Stable solutions

$$0 \le (3 + 2M_{22})^2 + 4M_{23}M_{32} < 9; \qquad (A18)$$

(ii) Stable spiral solutions

$$(3 + 2M_{22})^2 + 4M_{23}M_{32} < 0. \qquad (A19)$$

These are Eqs. (65) and (66) in Sec. VI.

## APPENDIX B: NONSTANDARD STABILITY ANALYSIS FOR CASE 1

In this appendix we provide the detailed computations supporting the linear stability analysis presented in Sec. VI A and the numerical analysis explained in Sec. VII A for the Case 1 fixed points.

To compute $\delta V$ from (73) and the limits (68) and (69), we split the analysis into 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 \to +\infty$ and $V_0 \to V_{FP}$. 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 Sec. VII.

### 1. Matrix P and its inverse in the fixed-point limit

To analyze the fixed-point limit of $\delta 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 = \begin{bmatrix} P_{11} & P_{12} & 0 & 0 & 0 \\ P_{21} & P_{22} & P_{23} & P_{24} & P_{25} \\ 0 & 1 & P_{33} & 1 & 1 \\ 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \end{bmatrix},$$

where

$$P_{11} = \frac{M_{25}M_{32}}{(M_{22} + 4)M_{31}}, \qquad (B1)$$

$$P_{12} = \frac{M_{23}M_{32}}{(M_{22} + 4)M_{31}} - \frac{4 + M_{33}}{M_{31}}, \qquad (B2)$$

$$P_{21} = -\frac{M_{25}}{(M_{22} + 4)}, \qquad (B3)$$

$$P_{22} = -\frac{M_{23}}{(M_{22} + 4)}, \qquad (B4)$$

$$P_{23} = -M_{24}(M_{33} + 3) \times (M_{22}M_{33} + 3M_{22} - M_{23}M_{32} + 3M_{33} + 9)^{-1}, \qquad (B5)$$

$$P_{24} = -2M_{23} \times \left( \sqrt{(M_{22} - M_{33})^2 + 4M_{23}M_{32}} + M_{22} - M_{33} \right)^{-1}, \qquad (B6)$$

$$P_{25} = 2M_{23} \times \left( \sqrt{(M_{22} - M_{33})^2 + 4M_{23}M_{32}} - M_{22} + M_{33} \right)^{-1}, \qquad (B7)$$

$$P_{33} = M_{24} \times ((M_{22}M_{33} + 3M_{22} + 3M_{33} + 9)M_{32}^{-1} - M_{23})^{-1}. \qquad (B8)$$

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,<fnref n="15" /> since both of them scale as $(\sim 1/y)$. All the other elements of $M$ are finite.

The nonzero elements of $P$ that have the terms $M_{31}$ and $M_{32}$ either depend on $M_{31}^{-1}$, $M_{32}^{-1}$ (both approaching 0) or contain the ratio $M_{32}/M_{31}$ that tends to

$$\frac{M_{32}}{M_{31}} \to -G_2'(X) \left( \frac{y}{G_2(X)} \right) \qquad (B9)$$

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

<footnote n="15">¹⁵We use the phrase “fixed-point limit” loosely in this appendix, since, as we discuss in Sec. VII, the Case 1 fixed point cannot be approached dynamically.</footnote>
$$ \mathbf{P}^{-1} = \begin{bmatrix} 0 & 0 & 0 & 0 & 1 \\ R_{21} & 0 & 0 & 0 & R_{25} \\ 0 & 0 & 0 & 1 & 0 \\ R_{41} & R_{42} & R_{43} & R_{44} & R_{45} \\ R_{51} & -R_{42} & R_{53} & R_{54} & R_{55} \end{bmatrix}, $$

where

$$ R_{21} = \frac{-(M_{22} + 4)}{D_1}, \qquad (B10) $$

$$ R_{25} = \frac{M_{25}M_{32}}{D_1}, \qquad (B11) $$

$$ R_{41} = \frac{M_{31}(N_{-})(8 + N_{+})}{4D_{1}D_{2}}, \qquad (B12) $$

$$ R_{42} = \frac{-M_{32}}{D_2}, \qquad (B13) $$

$$ R_{43} = \frac{N_{-}}{2D_{2}}, \qquad (B14) $$

$$ R_{44} = \frac{-M_{24}M_{32}(6 + N_{+})}{2D_{2}D_{3}}, \qquad (B15) $$

$$ R_{51} = \frac{-M_{31}(-8 + N_{-} - 2M_{22})}{4D_{2}D_{1}} \times (N_{+} - 2M_{22}), \qquad (B16) $$

$$ R_{53} = \frac{N_{+} - 2M_{22}}{2D_{2}}, \qquad (B17) $$

$$ R_{54} = \frac{M_{24}M_{32}(6 - N_{-} + 2M_{22})}{2D_{2}D_{3}}, \qquad (B18) $$

$$ R_{55} = \frac{M_{25}M_{32}(8 - N_{-} + 2M_{22})}{2D_{1}D_{2}}, \qquad (B19) $$

where

$$ D_1 = M_{31}M_{22}(M_{33} + 4) - M_{23}M_{32} + 4(M_{33} + 4), \qquad (B20) $$

$$ D_2 = \sqrt{(M_{22} - M_{33})^2 + 4M_{23}M_{32}}, \qquad (B21) $$

$$ D_3 = M_{22}(M_{33} + 3) - M_{23}M_{32} + 3M_{33} + 9, \qquad (B22) $$

$$ N_{-} = \sqrt{(M_{22} - M_{33})^2 + 4M_{23}M_{32} + M_{22} - M_{33}}, \qquad (B23) $$

$$ N_{+} = \sqrt{(M_{22} - M_{33})^2 + 4M_{23}M_{32} + M_{22} + M_{33}}, \qquad (B24) $$

It is easy to see that $D_1, D_3 \sim 1/y$ and $D_2, N_-, N_+ \sim 1/\sqrt{y}$, which makes $R_{41}, R_{42}, R_{51}$, and $R_{52}$ divergent at the fixed point.

## 2. Analysis of perturbations in the fixed-point limit

The perturbations $\delta V$ given by (73) are expanded to

$$ \delta y = \frac{f_1(\mathbf{V}_0)}{4} + \left( \frac{e^{-4N}}{4 + M_{22}} \right) \times \left[ M_{25} \frac{M_{32}}{M_{31}} C_1 - \left( \frac{16}{M_{31}} - M_{23} \frac{M_{32}}{M_{31}} + 4 \frac{M_{33}}{M_{31}} + \frac{M_{22}}{M_{31}} (4 + M_{33}) \right) C_2 \right], \qquad (B25) $$

$$ \begin{aligned} \delta X = & \frac{1}{12(M_{23} - M_{22} \frac{M_{33}}{M_{32}})} \\ & \times \left[ -3M_{23} \left( \frac{M_{31}}{M_{32}} f_1(\mathbf{V}_0) + \frac{4}{M_{32}} f_3(\mathbf{V}_0) \right) \right. \\ & \left. + \frac{M_{33}}{M_{32}} (12f_2(\mathbf{V}_0) + 4M_{24}f_4(\mathbf{V}_0) + 3M_{25}f_5(\mathbf{V}_0)) \right] \\ & - \left( \frac{e^{-4N}}{4 + M_{22}} \right) (M_{25}C_1 + M_{23}C_2) \\ & - M_{24} \left( \frac{3 + M_{33}}{D_3} \right) C_3 e^{-3N} \\ & - \left( \frac{2M_{23}C_4 e^{\lambda_4 N}}{N_{-}} - \frac{2M_{23}C_5 e^{\lambda_5 N}}{N_{+} - 2M_{22}} \right), \end{aligned} \qquad (B26) $$

$$ \begin{aligned} \delta \tilde{\Omega}_{\text{DE}} = & \frac{1}{12(M_{22} \frac{M_{33}}{M_{32}} - M_{23})} \\ & \times \left[ -3M_{22} \left( \frac{M_{31}}{M_{32}} f_1(\mathbf{V}_0) + \frac{4}{M_{32}} f_3(\mathbf{V}_0) \right) \right. \\ & \left. + (12f_2(\mathbf{V}_0) + 4M_{24}f_4(\mathbf{V}_0) + 3M_{25}f_5(\mathbf{V}_0)) \right] \\ & + C_2 e^{-4N} + \frac{M_{24}C_3 e^{-3N}}{D_3 M_{32}^{-1}} \\ & + (C_4 e^{\lambda_4 N} + C_5 e^{\lambda_5 N}), \end{aligned} \qquad (B27) $$

$$ \delta \tilde{\Omega}_m = \frac{f_4(\mathbf{V}_0)}{3} + C_3 e^{-3N}, \qquad (B28) $$

$$ \delta \tilde{\Omega}_r = \frac{f_5(\mathbf{V}_0)}{4} + C_1 e^{-4N}. \qquad (B29) $$

Since $\delta X_0$ and $\delta \tilde{\Omega}_{DE_0}$ 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 $\delta X$ real automatically.

We now study the limit of $\delta V$ for $N \to +\infty$ and $V_0 \to V_{FP}$. We identify two different possibilities.

(a) If and only if the eigenvalues $\lambda_4$ and $\lambda_5$ are either real and negative or complex with a 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 $f_i(\mathbf{V}_0)$ vanish for $V_0 \to \mathbf{V}_{FP}$, we have

$$ \lim_{\substack{N \to +\infty \\ V_0 \to V_{FP}}} \delta \mathbf{V} \to 0. \qquad (B30) $$

We conclude that, under the above specified conditions for the eigenvalues, the perturbations approach zero near the fixed point.

(b) If the eigenvalues $\lambda_4$ and $\lambda_5$ are real and positive (or complex with a positive real part), the last term of (B27) diverges (or is indeterminate) for $N \to \infty$ while the other terms remain finite. We conclude that $\delta\tilde{\Omega}_{DE}$ is divergent (or indeterminate).

From a close inspection to $\lambda_4$, $\lambda_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

$$ M_{23} = \frac{3X(1 - \omega/6)}{\omega \tilde{h}^2}, \qquad (B31) $$

at the fixed point that is positive if $0 < \omega < 6$ and negative when $\omega > 6$. This means one must satisfy the following condition for stable attractors:

$$ M_{32} \to \begin{cases} -\infty, & 0 < \omega < 6 \\ +\infty, & \omega > 6 \end{cases}. \qquad (B32) $$

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} \begin{cases} < 0, & 0 < \omega < 6 \\ > 0, & \omega > 6 \end{cases} \qquad (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_{32}^{\text{Div}} = -M_{31}G_2'(X_0)(\frac{y}{G_2(X_0)})$. As an illustration, consider the region $0 < \omega < 6$ and $\alpha_\sigma > 0$. From (24) and (14), we see that

$$ 1 - \left(\frac{y}{G_2(X)}\right)^2 > 0. \qquad (B34) $$

We can therefore write

$$ \text{Sign}(M_{31}) = \text{Sign}\left(\frac{G_1'(X_0)}{y}\right). \qquad (B35) $$

Plugging this expression in the divergent term of $M_{32}$ gives

$$ \text{Sign}(M_{32}^{\text{Div}}) = \text{Sign}\left(-\frac{G_1'(X_0)G_2'(X_0)}{G_2(X_0)}\right), \qquad (B36) $$

which must be negative. We arrive at the relation (74) in Sec. VI A

$$ G_1'(X_0)G_2(X_0)G_2'(X_0) > 0. \qquad (B37) $$

Similarly, we find that

(i) $\alpha_\sigma > 0$

$$ G_1'(X_0)G_2(X_0)G_2'(X_0) < 0, \quad \omega > 6; \qquad (B38) $$

(ii) and $\alpha_\sigma < 0$

$$ G_1'(X_0)G_2(X_0)G_2'(X_0) < 0, \quad 0 < \omega < 6, \qquad (B39) $$

$$ G_1'(X_0)G_2(X_0)G_2'(X_0) > 0, \quad \omega > 6, \qquad (B40) $$

which are (75), (76), and (77) in Sec. VI A. Therefore, since $G_2(X)$ changes sign at $X_{FP}$, in a neighborhood of a fixed point $V_{FP}$, there are regions that satisfy (a) and regions that satisfy (b).

As emphasized in Sec. VI A, for $V_{FP}$ to be an attractor we require that

$$ \lim_{\substack{N \to +\infty \\ V_0 \to V_{FP}}} \frac{d}{dN} \delta \mathbf{V} \to 0. \qquad (B41) $$

When condition (a) is satisfied, the last term of the $d\delta\tilde{\Omega}_{DE}/dN$ [see (B27)] is indeterminate in the limit.

The last requirement for $V_{FP}$ to be an attractor is that, in linear perturbation theory, the following limit

$$ \lim_{\substack{N \to +\infty \\ V_0 \to V_{FP}}} \left( \frac{y}{G_2(X)} \right)^2 \qquad (B42) $$

approaches dynamically its fixed-point limit given by (53). When condition (a) holds true, we can expand (B42) to

$$ \lim_{\substack{N \to +\infty \\ V_0 \to V_{FP}}} \left( \frac{y_0 + \delta y}{G_2(X_0) + G_2'(X_0)\delta X} \right)^2, \qquad (B43) $$

and substituting (B25) and (B26) we find that the limit is indeterminate.

### 3. Initial conditions for the numerical analysis

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, $\mathbf{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, \tilde{\Omega}_{DE_i}, \tilde{\Omega}_{m_i}, \tilde{\Omega}_{m_i})$ as follows:

$$X_i = X_{\pm} + \Delta X, \qquad (B44)$$

$$\tilde{\Omega}_{DE_i} = \frac{G_1(X_{\pm})}{(1 - \omega/6)} + \Delta\tilde{\Omega}_{DE}, \qquad (B45)$$

$$\tilde{\Omega}_{m_i} = \Delta\tilde{\Omega}_{m_i}, \qquad (B46)$$

$$\tilde{\Omega}_{r_i} = \Delta\tilde{\Omega}_{r_i}. \qquad (B47)$$

We observe that, if we set $y_i = \Delta y$, then the term $\frac{y}{G_2(X)} \equiv z(y, X)$ in (40) would not be close to its fixed-point value given by (53). This happens because, on expanding $z(y, X)$ around the fixed point, we obtain

$$z(\Delta y, X_{\pm} + \Delta X) = z(0, X_{\pm}) + \frac{\Delta y}{G_2(X_{\pm})} - z(0, X_{\pm}) \frac{G_2'(X_{\pm})}{G_2(X_{\pm})} \Delta X. \qquad (B48)$$

The second and third terms in the above expansion have a divergent factor of $\frac{1}{G_2(X_{\pm})}$. 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 + \frac{\omega G_1(X)}{X G_1'(X)(1-\omega/6)}) + \Delta r$ and use this value in (24) to set the initial value of $y$ as

$$y_i = \pm \sqrt{1 - \alpha_{\sigma} \frac{6(\tilde{\Omega}_{DE_i} - G_1(X_i))}{\omega r_i^2 X_i^2}} G_2(X_i). \qquad (B49)$$

## 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 $\mathcal{T} = tm_g$ and the following dimensionless quantities: $\tau = (\frac{\sigma'}{M_{Pl}})^2$, $\bar{\rho}_m = \frac{\rho_m}{M_{Pl}^2 m_g^2}$, $\bar{\rho}_r = \frac{\rho_r}{M_{Pl}^2 m_g^2}$ and $\tilde{h} = \frac{\dot{a}}{a}$, where the primes ($'$) represent derivative with respect to $\mathcal{T}$. The resulting background evolution equations are

$$y' = -4y\tilde{h}, \qquad (C1)$$

$$X' = X(\pm\sqrt{\tau} - \tilde{h}), \qquad (C2)$$

$$\tau' = -6 \left[ \tilde{h}\tau + \frac{XG_1'(X)}{\omega} (\pm\sqrt{\tau} - r\tilde{h}) \right], \qquad (C3)$$

$$\bar{\rho}'_m = -3\tilde{h}\bar{\rho}_m, \qquad (C4)$$

$$\bar{\rho}'_r = -4\tilde{h}\bar{\rho}_r, \qquad (C5)$$

where one must substitute

$$r = \sqrt{\frac{\tau\alpha_{\sigma}}{X^2(1 - (\frac{y}{G_2(X)})^2)}} \qquad (C6)$$

and

$$\tilde{h}^2 = \frac{\omega\tau}{6} + G_1(X) + \frac{\bar{\rho}_m}{3} + \frac{\bar{\rho}_r}{3}. \qquad (C7)$$

Notice that using $\mathcal{T}$ as the independent variable rather than $N$ provides us with a system of equations with no divergences as $\tilde{h}$ approaches zero, making it suitable to study the special fixed-point case $\tilde{h} = 0$. However, these equations are inconvenient for the general stability analysis, since the stability matrix has fewer diagonal terms.

### 1. Fixed-point solutions

We now focus on the $\tilde{h} = 0$ fixed-point analysis. One can easily verify from (C1), (C4), and (C5) that the derivatives of $y$, $\bar{\rho}_m$, and $\bar{\rho}_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 $\tau$ must vanish. There are three cases:

(i) If $\tau = 0$ and $X \neq 0$, then (C3) is automatically satisfied and from the Friedmann equation (C7)

$$G_1(X) + \frac{\bar{\rho}_m}{3} + \frac{\bar{\rho}_r}{3} = 0. \qquad (C8)$$

From (14), we get $n(t)^2 = (\dot{\phi}_0)^2$. From (24), we get $(\frac{y}{G_2(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 $\rho_m = \rho_r = 0$. Since the initial values of matter and radiation density are not zero, we find that $a \to \infty$.

(ii) If $X = 0$ but $\tau \neq 0$, then (C3) is automatically satisfied (meaning $\tau$ is indeterminate) and (C7) gives us

$$\frac{\omega\tau}{6} + G_1(0) + \frac{\bar{\rho}_m}{3} + \frac{\bar{\rho}_r}{3} = 0. \qquad (C9)$$

From (19), either $\sigma \to -\infty$ [in which case, from (14) and (20), $n$ and $r$ become divergent, respectively], or $a \to \infty$ [in which case, from (20), $r$ becomes divergent].
(iii) If $\tau = X = 0$, then (C7) gives us

$$G_1(0) + \frac{\bar{\rho}_m}{3} + \frac{\bar{\rho}_r}{3} = 0. \quad (C10)$$

From (19), either $\sigma \to -\infty$ [in which case, from (14), $\dot{\phi}_0/n$ is indeterminate], or $a \to \infty$ [in which case $n^2 = (\dot{\phi}^0)^2$].

We then conclude that $\Omega_{DE}$ balances the sum of matter and radiation energy densities.

## 2. Fixed-point linear stability

To study the linear stability of the fixed points we follow (C1)–(C5) defining $\mathbf{V} = [y, X, \tau, \bar{\rho}_m, \bar{\rho}_r]^T$ as the dynamical variables and $\mathbf{f}(\mathbf{V})$ as the RHS of these evolution equations. The stability matrix $\mathbf{M}$ at any time $T$ takes the form

$$M = \begin{bmatrix} M_{11} & M_{12} & M_{13} & M_{14} & M_{15} \\ 0 & M_{22} & M_{23} & M_{24} & M_{25} \\ M_{31} & M_{32} & M_{33} & M_{34} & M_{35} \\ 0 & M_{42} & M_{43} & M_{44} & M_{45} \\ 0 & M_{52} & M_{53} & M_{54} & M_{55} \end{bmatrix},$$

where

$$M_{11} = -4\tilde{h}, \quad (C11)$$

$$M_{12} = \frac{-2yG_1'(X)}{\tilde{h}}, \quad (C12)$$

$$M_{13} = \frac{-y\omega}{3\tilde{h}}, \quad (C13)$$

$$M_{14} = \frac{-2y}{3\tilde{h}}, \quad (C14)$$

$$M_{15} = M_{14}, \quad (C15)$$

$$M_{22} = \pm\sqrt{\tau} - \tilde{h} - \frac{XG_1'(X)}{2\tilde{h}}, \quad (C16)$$

$$M_{23} = \frac{X}{2} \left( \pm \frac{1}{\sqrt{\tau}} - \frac{\omega}{6\tilde{h}} \right), \quad (C17)$$

$$M_{24} = \frac{-X}{6\tilde{h}}, \quad (C18)$$

$$M_{25} = M_{24}, \quad (C19)$$

$$M_{31} = \frac{6y\tilde{h} \sqrt{\frac{\alpha_{\sigma}\tau}{(1-(y/G_2(X))^2)^3}} G_1'(X)}{\omega(G_2(X))^2}, \quad (C20)$$

$$M_{32} = 6 \frac{\sqrt{\tau}}{\omega} \left[ \frac{(G_1'(X))^2}{2\tilde{h}} \sqrt{\frac{\alpha_{\sigma}}{(1-(y/G_2(X))^2)}} - G_1'(X) \left( \pm 1 + \frac{\sqrt{\tau\omega}}{2\tilde{h}} + \frac{y^2\tilde{h}G_2'(X)\left(\frac{\alpha_{\sigma}}{(1-(y/G_2(X))^2)}\right)^{3/2}}{\alpha_{\sigma}(G_2(X))^3} \right) - \left( \pm X - \tilde{h} \sqrt{\frac{\alpha_{\sigma}}{(1-(y/G_2(X))^2)}} \right) G_1''(X) \right], \quad (C21)$$

$$M_{33} = -6 \left( \tilde{h} + \frac{(\pm X - \tilde{h} \sqrt{\frac{\alpha_{\sigma}}{(1-(y/G_2(X))^2)}}) G_1'(X)}{2\sqrt{\tau\omega}} \right) - \frac{\tau\omega}{2\tilde{h}} + \sqrt{\frac{\alpha_{\sigma}\tau}{(1-(y/G_2(X))^2)}} \frac{G_1'(X)}{2\tilde{h}}, \quad (C22)$$

$$M_{34} = \frac{-\tau}{\tilde{h}} + \frac{G_1'(X)}{\omega\tilde{h}} \sqrt{\frac{\tau\alpha_{\sigma}}{(1-(y/G_2(X))^2)}}, \quad (C23)$$

$$M_{35} = M_{34}, \quad (C24)$$

$$M_{42} = \frac{-3\bar{\rho}_m G_1'(X)}{2\tilde{h}}, \quad (C25)$$

$$M_{43} = \frac{-\omega\bar{\rho}_m}{4\tilde{h}}, \quad (C26)$$

$$M_{44} = \frac{-3}{\tilde{h}} \left( \tilde{h}^2 + \frac{\bar{\rho}_m}{6} \right), \quad (C27)$$

$$M_{45} = \frac{-\bar{\rho}_m}{2\tilde{h}}, \quad (C28)$$

$$M_{52} = \frac{-2\bar{\rho}_r G_1'(X)}{\tilde{h}}, \quad (C29)$$

$$M_{53} = \frac{-\omega\bar{\rho}_r}{3\tilde{h}}, \quad (C30)$$

$$M_{54} = \frac{-2\bar{\rho}_r}{3\tilde{h}}, \quad (C31)$$

$$M_{55} = \frac{-4}{\tilde{h}} \left( \tilde{h}^2 + \frac{\bar{\rho}_r}{6} \right). \quad (C32)$$

We cannot find the eigenvalues of $\mathbf{M}$ analytically, because the characteristic polynomial is fifth order. Therefore we cannot apply the method described in Sec. VI A. We instead study the characteristic polynomial in the fixed-point limit. We find that in both fixed-point cases, $\tau = 0$ and $X = 0$, the polynomial has indeterminate coefficients [due to the presence of ratios of powers of $\tilde{h}$, $\tau$, or $1 - (\frac{y}{G_2(X)})^2$] or divergent coefficients [due to inverses of powers of $\tilde{h}$, $\tau$, or $1 - (\frac{y}{G_2(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 an analysis as out of the scope of our present work.

## APPENDIX D: SCALAR PERTURBATIONS

Since there are certain disagreements between our equations for the perturbations and those in the literature, for the sake of completeness, here we present the pertinent equations.

Regarding scalar sector, the discrepancy we find can be noted by comparing the constraint equations obtained from the variation of the action quadratic in the perturbations, with respect to $\Phi$ and $B$, and we obtain (in Fourier space)

$$ 
\begin{aligned}
0 ={}& 6H\dot{\psi} + \frac{2k^2}{a^2}\left(\psi + \frac{k^2 E}{6} + aBH\right) - \frac{\omega\dot{\sigma}\delta\sigma}{M_{\text{Pl}}} \\
& - \frac{\delta\chi\rho_{,\chi}}{M_{\text{Pl}}} - \frac{\dot{\delta\chi}(p+\rho)}{\dot{\chi}c_s^2 M_{\text{Pl}}} + \frac{3\psi}{M_{\text{Pl}}^2}\left(3H^2M_{\text{Pl}}^2 - \rho - \frac{\omega}{2}\dot{\sigma}^2\right) \\
& + 3\psi m_g^2(J + \alpha_3(X-1)(2X-1) + (X-6)X + 3) \\
& + \frac{3\Phi}{M_{\text{Pl}}^2}\left(\frac{p+\rho(1+c_s^2)}{3c_s^2} - 3H^2M_{\text{Pl}}^2 + \frac{\omega}{2}\dot{\sigma}^2\right) \\
& + \Phi(X-1)m_g^2(J + \alpha_3(X-1)^2 - 3X + 3) \\
& - 3Xm_g^2(J(X) + \alpha_3(X-1)X - 2X)\delta\sigma,
\end{aligned}
\quad (\text{D1})
$$

$$ 
\begin{aligned}
0 ={}& \frac{3\delta\chi(r^2 - 1)(p + \rho)}{a\dot{\chi}\omega M_{\text{Pl}}} - \frac{6(r^2 - 1)}{a\omega} \times \left(H\Phi - \dot{\psi} - \frac{k^2}{6}\dot{E}\right) \\
& + \frac{3}{\omega}B(r-1)m_g^2 \times \left[\alpha_3(X-1)(r(X-1)^2 - 2X + 1) - 3r(X-1)^2 - X^2 + 6X - 3\right] \\
& + \frac{3B(r^2 - 1)}{\omega M_{\text{Pl}}^2} \\
& \times \left(-3H^2M_{\text{Pl}}^2 + \rho + \frac{\omega}{2}\dot{\sigma}^2\right) \\
& + \frac{3}{\omega}BJ(r-1) \\
& \times m_g^2(r(X-1) - 1) + \frac{3\delta\sigma\dot{\sigma}}{aX\omega M_{\text{Pl}}} \left[(r^2 - 1)X\omega - (r-1)\alpha_\sigma \times (J + \alpha_3(X-1)X - 2X)\right].
\end{aligned}
\quad (\text{D2})
$$

In deriving these equations we have not made use of the background equations. One can show that using the assumed background evolution given in (28) and (29) of [12], these equations are equivalent to (48) of [12]. Even after using the background equations, one can show that our Eq. (D2) differs from Eq. (4.22) of [11] (we have an additional contribution proportional to $\alpha_\sigma\delta\sigma$), and (48) of [12] also differs from those in [11] in their regime of applicability. For this reason we believe our equations are correct.

## APPENDIX E: STABILITY OF PERTURBATIONS

In this appendix we provide the conditions obtained from requiring the stability of the perturbations, which we used to restrict the values of the parameters of the theory in the numerical study of Sec. VII.

The mass of tensor perturbations is given by

$$ m_{\text{GW}}^2 = m_g^2 X(3 + 3\alpha_3 + \alpha_4 + rX^2(\alpha_3 + \alpha_4) - (1 + r)X(1 + 2\alpha_3 + \alpha_4)). \quad (\text{E1}) $$

To avoid tachyonic instabilities, we require $m_{\text{GW}}^2 > -H^2$. For tensor perturbations, we do not reproduce Eq. (4.7) of [11], nor the underlying equation of motion for the $\sigma$ field, Eq. (3.13). Equation (E1) agrees with [8,11] at the Case 1 fixed point, but not elsewhere.

For vector perturbations, we confirm the expression [Eq. (4.17)] in [11], for the coefficient of the vector kinetic term, and rewrite it as

$$ \kappa_V^2 = \frac{2m_g^2 a^2 X G_1'(X)}{(k^2(1+r) + 2m_g^2 a^2 X G_1'(X))}. \quad (\text{E2}) $$

For the perturbations to be stable, we must have $k_V^2 > 0$ and a vector mass squared that limits any tachyonic instability $m_V^2 = \frac{m_{\text{GW}}^2}{\kappa_V^2} > -H^2$. We see from (E2) that if we impose

$$ G_1'(X) > 0, \quad (\text{E3}) $$

the kinetic coefficient obeys $0 < \kappa_V^2 < 1$. If $m_{\text{GW}}^2 > -H^2$, as imposed above, then this also ensures tachyonic stability for vector perturbations.

For Case 1, after integrating out the BD mode, we obtain the conditions for stability of the scalar perturbations. Our results for this case reproduce the ones in [12]. At this fixed point, (E1) reduces to

$$ m_{\text{GW}}^2 = \frac{m_g^2(r-1)^2 X^3 + \omega H^2(r(X+1)-2)}{(r-1)(X-1)}. \quad (\text{E4}) $$

Following [12], we arrive at the stability conditions for scalar perturbations,

$$ X^2 < \alpha_\sigma \tilde{h}^2 < r^2 X^2. \quad (\text{E5}) $$

It follows that

$$ r > 1 \quad (\text{E6}) $$

and

$$ \alpha_\sigma > 0. \quad (\text{E7}) $$

For simplicity, we will not evaluate the stability conditions for scalar perturbations in the region $\omega \ge 6$ for either $\alpha_\sigma > 0$ or $\alpha_\sigma < 0$.
[1] C. de Rham, G. Gabadadze, and A. J. Tolley, Resummation of Massive Gravity, Phys. Rev. Lett. **106**, 231101 (2011).

[2] S. F. Hassan and R. A. Rosen, Resolving the Ghost Problem in Nonlinear Massive Gravity, Phys. Rev. Lett. **108**, 041101 (2012).

[3] G. D'Amico, C. de Rham, S. Dubovsky, G. Gabadadze, D. Pirtskhalava, and A. J. Tolley, Massive cosmologies, Phys. Rev. D **84**, 124046 (2011).

[4] G. D'Amico, G. Gabadadze, L. Hui, and D. Pirtskhalava, Quasidilaton: Theory and cosmology, Phys. Rev. D **87**, 064037 (2013).

[5] S. Anselmi, D. López Nacir, and G. D. Starkman, Extreme parameter sensitivity in quasidilaton massive gravity, Phys. Rev. D **92**, 084033 (2015).

[6] G. D'Amico, G. Gabadadze, L. Hui, and D. Pirtskhalava, On cosmological perturbations of quasidilaton, Classical Quantum Gravity **30**, 184005 (2013).

[7] A. E. Güm>rükçüoğlu, K. Hinterbichler, C. Lin, S. Mukohyama, and M. Trodden, Cosmological perturbations in extended massive gravity, Phys. Rev. D **88**, 024023 (2013).

[8] A. De Felice and S. Mukohyama, Towards consistent extension of quasidilaton massive gravity, Phys. Lett. B **728**, 622 (2014).

[9] J. Klusoň, Consistent extension of quasidilaton massive gravity, J. Grav. **2014**, 413835 (2014).

[10] S. Mukohyama, Extended quasidilaton massive gravity is ghost free, arXiv:1309.2146 v1

[11] L. Heisenberg, Revisiting perturbations in extended quasidilaton massive gravity, J. Cosmol. Astropart. Phys. **04** (2015) 010.

[12] A. E. Gumrukcuoglu, K. Koyama, and S. Mukohyama, Role of matter in extended quasidilaton massive gravity, Phys. Rev. D **94**, 123510 (2016).

[13] D. Comelli, M. Crisostomi, K. Koyama, L. Pilo, and G. Tasinato, New branches of massive gravity, Phys. Rev. D **91**, 121502 (2015).

[14] P. Motloch, W. Hu, A. Joyce, and H. Motohashi, Self-accelerating massive gravity: Superluminality, Cauchy surfaces, and strong coupling, Phys. Rev. D **92**, 044024 (2015).

[15] R. Gannouji, M. W. Hossain, M. Sami, and E. N. Saridakis, Quasidilaton nonlinear massive gravity: Investigations of background cosmological dynamics, Phys. Rev. D **87**, 123536 (2013).

[16] T. Kahniashvili, A. Kar, G. Lavrelashvili, N. Agarwal, L. Heisenberg, and A. Kosowsky, Cosmic expansion in extended quasidilaton massive gravity, Phys. Rev. D **91**, 041301 (2015).

[17] E. J. Copeland, M. Sami, and S. Tsujikawa, Dynamics of dark energy, Int. J. Mod. Phys. D **15**, 1753 (2006).

[18] C. Armendariz-Picon, V. F. Mukhanov, and P. J. Steinhardt, Essentials of *k* essence, Phys. Rev. D **63**, 103510 (2001).

[19] S. Mukohyama, A new quasidilaton theory of massive gravity, J. Cosmol. Astropart. Phys. **12** (2014) 011.

[20] A. De Felice, A. E. Güm>rükçüoğlu, and S. Mukohyama, Generalized quasidilaton theory, Phys. Rev. D **88**, 124006 (2013).

[21] Q.-G. Huang, Y.-S. Piao, and S.-Y. Zhou, Mass-varying massive gravity, Phys. Rev. D **86**, 124014 (2012).

[22] C. de Rham, M. Fasiello, and A. J. Tolley, Stable FLRW solutions in generalized massive gravity, Int. J. Mod. Phys. D **23**, 1443006 (2014).

[23] C. Pitrou, X. Roy, and O. Umeh, xPand: An algorithm for perturbing homogeneous cosmologies, Classical Quantum Gravity **30**, 165002 (2013).