Tracker and scaling solutions in DHOST theories

In quadratic-order degenerate higher-order scalar-tensor (DHOST) theories compatible with gravitational-wave constraints, we derive the most general Lagrangian allowing for tracker solutions characterized by $\dot{\phi}/H^p={\rm constant}$, where $\dot{\phi}$ is the time derivative of a scalar field $\phi$, $H$ is the Hubble expansion rate, and $p$ is a constant. While the tracker is present up to the cubic-order Horndeski Lagrangian $L=c_2X-c_3X^{(p-1)/(2p)} \square \phi$, where $c_2, c_3$ are constants and $X$ is the kinetic energy of $\phi$, the DHOST interaction breaks this structure for $p \neq 1$. Even in the latter case, however, there exists an approximate tracker solution in the early cosmological epoch with the nearly constant field equation of state $w_{\phi}=-1-2p\dot{H}/(3H^2)$. The scaling solution, which corresponds to $p=1$, is the unique case in which all the terms in the field density $\rho_{\phi}$ and the pressure $P_{\phi}$ obey the scaling relation $\rho_{\phi} \propto P_{\phi} \propto H^2$. Extending the analysis to the coupled DHOST theories with the field-dependent coupling $Q(\phi)$ between the scalar field and matter, we show that the scaling solution exists for $Q(\phi)=1/(\mu_1 \phi+\mu_2)$, where $\mu_1$ and $\mu_2$ are constants. For the constant $Q$, i.e., $\mu_1=0$, we derive fixed points of the dynamical system by using the general Lagrangian with scaling solutions. This result can be applied to the model construction of late-time cosmic acceleration preceded by the scaling $\phi$-matter-dominated epoch.


I. INTRODUCTION
There have been numerous attempts to modify or extend General Relativity (GR) at large distances [1][2][3][4][5][6]. One of such motivations is to explain the observational evidence of late-time cosmic acceleration by introducing a new ingredient beyond the scheme of standard model of particle physics. The simple candidate for such a new degree of freedom (DOF) is a scalar field φ [7][8][9][10][11][12][13], which has been widely exploited to describe the dynamics of dark energy.
The theories in which the scalar field is directly coupled to gravity (with two tensor polarized DOFs) are generally called scalar-tensor theories [14,15]. It is known that Horndeski theories [16] are the most general scalar-tensor theories with second-order equations of motion [17][18][19]. The second-order property ensures the absence of an Ostrogradsky instability [20] associated with a linear dependence of the Hamiltonian arising from extra DOFs.
Horndeski theories can be extended to more general theoretical schemes without increasing the number of propagating DOFs [21]. For example, Gleyzes-Langlois-Piazza-Vernizzi (GLPV) expressed the Horndeski Lagrangian in terms of scalar quantities arising in the 3+1 decomposition of spacetime [22] and derived a beyond-Horndeski Lagrangian without imposing two conditions Horndeski theories obey [23]. The Hamiltonian analysis in the unitary gauge showed that the GLPV theories do not increase the number of DOFs relative to those in Horndeski gravity [24][25][26].
One can further perform a healthy extension of Horndeski theories by keeping one scalar and two tensor DOFs. Even if Euler-Lagrange equations contain derivatives higher than second order in the scalar field and the metric, it is possible to maintain the same number of propagating DOFs by imposing the so-called degeneracy conditions of their Lagrangians [27][28][29][30][31]. They are dubbed degenerate higher-order scalar-tensor (DHOST) theories, which encompass GLPV theories as a special case. The absence of an extra DOF was confirmed by the Hamiltonian analysis [26,28] as well as by the field definition linking to Horndeski theories [29,30,32].
The DHOST theories contain the products of covariant derivatives of the field which are quadratic and cubic in ∇ µ ∇ ν φ, say, ( φ) 2 and ( φ) 3 , respectively. If we apply the DHOST theories to dark energy and adopt the bound of the speed c t of gravitational waves constrained from the GW170817 event [33] together with the electromagnetic counterpart [34], the Lagrangians consistent with c t = 1 (in the unit where the speed of light is equivalent to 1) are up to quadratic in ∇ µ ∇ ν φ with one of the terms vanishing (A 1 = 0) [35] among six coefficients of derivative interactions. From the degeneracy conditions there are three constraints among the other five coefficients [36][37][38], so we are left with two quadraticorder free functions. If we take into account the decay of gravitational waves to dark energy [39], we have an ad-ditional constraint on the quadratic-order Lagrangian 1 . Hence there is one free DHOST interaction containing the term B 4 (φ, X)R (where R is the Ricci scalar) besides the Horndeski Lagrangian L = G 2 (φ, X) − G 3 (φ, X) φ up to cubic order.
If we apply shift-symmetric Horndeski theories to dark energy, there are self-accelerating solutions preceded by a constant tracker equation of state w φ withφ ∝ H p (p is a constant). For example, the covariant Galileon [41,42] gives rise to the value w φ = −2 with p = −1 during the matter era [43,44], but it is disfavored from the joint data analysis of supernovae type Ia, cosmic microwave background, and baryon acoustic oscillations [45]. The extended Galileon proposed in Ref. [46] can accommodate the tracker equation of state w φ closer to −1, in which case the model can be consistent with the observational data [47]. In DHOST theories, (approximate) tracker solutions were found for particular models [38,48], but the general conditions for its existence have been unknown.
In Horndeski theories, there is a special kind of tracker called the scaling solution [49][50][51][52][53][54][55][56][57][58][59][60][61][62] along which the field density ρ φ is proportional to the background matter density ρ m . If the scalar field has a constant coupling Q with matter, the scaling solution satisfying the relatioṅ φ ∝ H exists for the cubic-order Horndeski Lagrangian [63]. In this case, it is possible to construct viable dark energy models with a scaling φ-matter dominated epoch (φMDE). In DHOST theories, the conditions for realizing the scaling solution were not derived yet.
In this paper, we will derive the Lagrangian allowing for tracking and scaling solutions in DHOST theories compatible with observational constraints of gravitational waves. We impose the conditionφ ∝ H p by assuming that the quantity h = −Ḣ/H 2 is nearly constant. For p = 1, we show the existence of approximate tracker solutions characterized by the field equation of state w φ ≃ −1 + 2ph/3 in the early cosmological epoch. The scaling solution with the power p = 1 is the special case in which the exact scaling behavior of the field density (ρ φ ∝ ρ m ∝ H 2 ) can be realized without assuming the dominance of ρ m over ρ φ . We also extend the analysis to the case in which a field-dependent coupling Q(φ) between φ and matter is present and obtain the most general Lagrangian allowing for scaling solutions.
This paper is organized as follows. In Sec. II, we derive the background equations of motion in DHOST theories in the presence of the field-dependent coupling Q(φ) with matter. In Sec. III, we constrain the forms of DHOST Lagrangians allowing for the existence of tracker and scaling solutions for Q = 0. In Sec. IV, we obtain the most general Lagrangian with scaling solutions for the fielddependent coupling Q(φ) and also derive the fixed points of the dynamical system for constant Q. Sec. V is devoted to conclusions.

II. BACKGROUND EQUATIONS IN DHOST THEORIES
Let us consider the quadratic-order DHOST theories given by the action [27][28][29][30][31]: where g is the determinant of metric tensor g µν , R is the Ricci scalar, and 2) Here, we use the unit where the reduced Planck mass M pl is equivalent to 1. The functions G 2 , G 3 , B 4 , A 4 depend on φ and X = −∇ µ φ∇ µ φ/2, with the covariant derivative operator ∇ µ and the d'Alembert operator = ∇ µ ∇ µ . The quantity Z is defined by The function A 4 is related to B 4 according to with the notation B 4,X ≡ ∂B 4 /∂X. The full DHOST theories contain the other four Lagrangians Requiring that the speed c t of gravitational waves is equivalent to 1, it follows that A 1 = 0 [35]. The degeneracy conditions constrain the coupling A 2 to be 0 and the functions A 3 and A 5 are related to each other according to To avoid the decay of gravitational waves to dark energy perturbations [39], these functions are constrained to be A 3 = 0 = A 5 . On using the other degeneracy condition, we end up with the Lagrangian (2.2) with the particular relation (2.4).
The theory (2.2) can be also obtained from the cubicorder Horndeski Lagrangian L = P (φ, X)+Q(φ, X) φ+ f (φ)R after performing an invertible conformal transformation g µν → C(φ, X)g µν [39]. We also note that the GLPV theories [23] correspond to A 3 = −A 4 = −B 4,X /X = 0 and A 5 = 0, so they do not belong to the Lagrangian (2.2). In other words, the GLPV theories do not satisfy the bound arising from the decay of gravitational waves to dark energy.
For the matter action S m , we consider a barotropic perfect fluid which can be coupled to the scalar field φ. In scalar-tensor theories without vector propagating DOFs, the matter sector can be described by the Schutz-Sorkin action [63][64][65][66][67][68]: where the matter density ρ m depends on the fluid number density n as well as on φ. The four vector J µ is related to n, as n = J µ J ν g µν /g, and ℓ is a scalar quantity. We define the coupling between φ and matter, as where ρ m,φ ≡ ∂ρ m /∂φ. To study the background cosmological dynamics, we consider a flat Friedmann-Lemaître-Robertson-Walker spacetime given by the line element where N (t) is a lapse, and a(t) is a scale factor. The Lagrangian in the action (2.1) is given by with √ −g = N a 3 and H ≡ȧ/a, where a dot represents the derivative with respect to t ≡ N dt.
For the matter sector, the temporal component of J µ is related to the background number density n 0 , as J 0 = n 0 a 3 , so the matter action is expressed as The variations of S m with respect to n 0 and ℓ lead tȯ ℓ = −ρ m,n where ρ m,n ≡ ∂ρ m /∂n, anḋ respectively. Varying the total action (2.1) with respect to N and a, we obtain the modified Friedmann equations: Here, ρ φ and P φ correspond to the field density and pressure defined, respectively, by where The expressions of ρ φ and P φ derived above are valid even for DHOST theories with non-vanishing functions A 3 and A 5 , in which case the function B 1 is given by [38]. Since we are now considering the theories with A 3 = 0, B 1 is directly related to B 4 according to Eq. (2.15).
On using the propertyρ m = ρ m,nṅ0 +Q(φ)ρ mφ and the matter pressure (2.16), the conservation (2.10) of total fluid number translates tȯ where w m = P m /ρ m . Varying the total action S with respect to φ, it follows thaṫ where w φ = P φ /ρ φ . One can also derive Eq. (2.18) by taking the time derivative of Eq. (2.13) and using Eqs. (2.14) and (2.17). From Eq. (2.11), the density parameters Ω φ = ρ φ /(3H 2 ) and Ω m = ρ m /(3H 2 ) obey where w eff is the effective equation of state defined by We note that there are time derivatives ... φ in Eqs. (2.11)-(2.12) as well as .... φ andḦ in Eq. (2.18). As we will discuss in Sec. IV C, however, the background equations reduce to the dynamical system containing the time derivatives of φ and a up to second order thanks to the degeneracy conditions.

III. TRACKER AND SCALING SOLUTIONS FOR Q = 0
We derive the Lagrangian L allowing for the tracking solution satisfyingφ where p and α are constants. We focus on the case and impose the condition so that P φ scales in the same way as ρ φ . The tracker solution found for covariant Galileons [43,44] corresponds to p = −1, with the field equation of state w φ = −2 during the matter era. The scaling solution found for cubic-order Horndeski theories [63] corresponds to p = 1, with w φ = w m . Now, we are extending the analysis to a more general power p.
We take into account the canonical kinetic term X in G 2 and search for the theories in which each term in ρ φ and P φ evolves in the same way as X, i.e., The terms associated with the couplings G 2 and B 4 in the Lagrangian (2.8) appear in the expressions of ρ φ or P φ . Moreover, the G 3 -dependent contributions to Eq. (2.8) reduce to the term −φ 2 (G 3,φ +φ G 3,X ) in P φ after the integration by parts. Then, the Lagrangian should follow the same time dependence as ρ φ and P φ , i.e., In the following, we obtain the form of the Lagrangian allowing for the property (3.5). Since there are terms in ρ φ and P φ which are absent in L, we need to confirm whether each term in ρ φ and P φ obeys the property (3.4) after deriving the Lagrangian satisfying the condition (3.5). The relation (3.5) translates tȯ On using the relation (3.1), the quantities associated with the time derivatives of φ, X, φ = −φ − 3Hφ, R = 6(2H 2 +Ḣ), and Z = −φ 2φ2 can be expressed, respectively, asφ Substituting the Lagrangian (2.2) into Eq. (3.7) and treating A 4 as an independent function from B 4 , it follows that the couplings G = G 2 , G 3 , B 4 , A 4 need to separately obey the partial differential equations: where s is a constant given by In the following, we discuss the cases p = 1 and p = 1, separately.
A. Tracker solutions: p = 1 For p = 1 (i.e., n = 0), the general solution to Eq. (3.14) is given by where g is an arbitrary function of Since we are now considering the case in which h = −Ḣ/H 2 is approximately constant, the integration of this relation gives where t 0 is a constant. Then, we can integrate Eq. (3.1) to give where φ 0 is an integration constant. Since X n = 2 −n α 2n h 1−p (t − t 0 ) 1−p , it follows that Y = nλφ 0 = constant. Then, the function g(Y) does not vary in time along the tracker solution. Let us study whether each term in ρ φ and P φ following from the Lagrangian (3.16) obeys the property (3.4). First of all, the quadratic Lagrangian is given by G 2 = Xg 2 (Y), where g 2 (Y) is an arbitrary function of Y. Since g 2 (Y) does not vary in time along the tracker solution, the term G 2 in ρ φ and P φ evolves as G 2 ∝ X ∝ H 2p . The contributionφ 2 G 2,X to ρ φ has the dependencė φ 2 G 2,X =φ 2 (g 2 + nX n g 2,Y ), so it satisfies the property (3.4) for g 2,Y = 0, i.e., g 2 (Y) = c 2 = constant. Hence the quadratic Lagrangian obeying the relation (3.4) is constrained to be (3.20) The integrated solution to (3.14) for the cubic Lagrangian is given by G 3 = X n g 3 (Y). In order to have the relatioṅ φ 2 G 3,φ ∝φ 2 , we require that G 3,φ = nλX n g 3,Y does not change in time and hence g 3 (Y) = c 3 = constant. This restricts the Lagrangian to the form In this case, both the terms 3Hφ 3 G 3,X and −φ 2φ G 3,X are proportional toφ 2 , so all the cubic-order contributions to ρ φ and P φ satisfy the relation (3.4). The cubic Galileon G 3 = c 3 X corresponds to p = −1, in which case the tracker solution characterized byφH = constant is present during the radiation-and matter-dominated epochs [43,44]. The coupling B 4 following from the solution (3.16) is given by B 4 = X 2n b 4 (Y). The term HφB 4,φ in ρ φ and P φ is in proportion to H 3p−1 b 4,Y . Since we are considering the case p = 1, this term is consistent with the dependence (3.4) only for b 4 (Y) = c 4 = constant. Then we have (3.23) The field density ρ φ and the pressure P φ contain the terms like H 2 B 4 proportional to H 2p . On the other hand, there exists the term which does not behave as ∝ H 2p for p = 1. Under the condition |2c 4 X 2n | ≪ 1, there are also contributions to ρ φ and P φ proportional to c 3 4 H 6p−4 . If we demand the exact tracking behavior along which all the terms in ρ φ and P φ have the dependence ∝ H 2p , we have c 4 = 0 and hence For c 4 = 0, this is at odds with the integrated solution A 4 = X −(p+1)/p a 4 (Y). Even in the case c 4 = 0, there exists an approximate tracker solution in the early cosmological epoch. Provided that |2c 4 X 2n | ≪ 1, the leading-order terms to ρ φ and P φ for p < 1 correspond to those proportional to H 2p , which arise from the couplings G 2 = c 2 X, G 3 = c 3 X (p−1)/(2p) as well as B 4 = c 4 X (p−1)/p . The existence of the coupling B 4 = c 4 X (p−1)/p gives rise to terms with the different power-law dependence of H. The nextto-leading contributions to ρ φ and P φ are in proportion to H 4p−2 . Then, it follows that (3.27) where α 1,2 and β 1,2 are constants, and the abbreviation means the terms which are next order to H 4p−2 . Substituting Eqs. (3.27)-(3.28) and the time derivativeρ φ into the continuity equationρ φ + 3H(ρ φ + P φ ) = 0, we can solve it for β 1 . On using this relation, the field equation of state w φ = P φ /ρ φ yields where we picked up the terms up to the order H 4p−2 in Eqs. (3.27) and (3.28). For c 4 = 0, we have w φ = −1 + 2ph/3. Indeed, the cubic Galileon corresponds to p = −1, in which case w φ = −7/3 during the radiation dominance (h = 2) and w φ = −2 during the matter dominance (h = 3/2) [43,44]. For general values of p, we have w φ = −1 + p during the matter era. In this case, for p closer to 0, the model can be compatible with the observational data associated with the background expansion history.
The non-vanishing coupling B 4 = c 4 X (p−1)/p gives rise to the variation of w φ . For p < 1, the terms in the parentheses of Eqs. (3.27) and (3.28) evolve faster than H 2p , so they are suppressed relative to the former in the asymptotic past. In this limit, we recover the tracker equation of state w φ = −1 + 2ph/3. As long as the terms in the parentheses of Eqs. (3.27) and (3.28) catch up with their first terms, w φ starts to deviate from the tracker value −1 + 2ph/3. Thus, in the presence of the coupling B 4 = c 4 X (p−1)/p , the tracking behavior can be approximately realized in the early cosmological epoch during which the terms proportional to H 4p−2 and H 6p−4 are subdominant to the H 2p contributions to ρ φ and P φ . From Eq. (3.29), we observe that, in the limit p → 1, the field equation of state reduces to the tracker value w φ → −1 + 2h/3 = w eff even in the presence of the DHOST term B 4 = c 4 X (p−1)/p . This limit corresponds to the scaling solution along which w φ is equivalent to w m with constant Ω φ . For p = 1, the solution to Eq. (3.14) is different from Eq. (3.16), so we discuss this case separately in the following.

B. Scaling solutions: p = 1
If p = 1, then the solution to Eq. (3.14) is given by where s is given by Eq. (3.15), and g is an arbitrary function of and λ is a constant. Each coefficient in the Lagrangian can be written in the form:

34)
From the degeneracy condition (2.4), the function a 4 (Y ) is determined from b 4 (Y ), as Unlike the case p = 1, the Lagrangian A 4 derived above is consistent with the other scaling Lagrangian B 4 . From Eq. (2.15), the quantity B 1 is given by which depends on Y alone. For the solution satisfying the condition (3.1), the scalar field evolves as where φ 0 is an integration constant. Since X ∝ H 2 and e λφ ∝ a 2h ∝ H −2 , the quantity Y remains constant.

IV. SCALING LAGRANGIAN FOR GENERAL FIELD-DEPENDENT COUPLING Q(φ)
In Sec. III, we showed that the power p = 1 is the special case in which all the terms of the background equations of motion scale in the same manner (ρ φ ∝ P φ ∝ H 2 ) for Q = 0. Now, we extend the analysis to the field-dependent coupling Q(φ) and derive the Lagrangian whose equations of motion obey the scaling relations ρ φ ∝ P φ ∝ H 2 . Since the scaling solution satisfies the relation ρ φ /ρ m = constant, both Ω φ and Ω m are constant. Then, the effective equation of state w eff and the quantity h = −Ḣ/H 2 do not vary in time in the scaling regime.
In the scaling regime, the quantityλ is constant. From Eq. (4.1), the field derivative has the dependenceφ ∝ H/Q(φ).
As we mentioned in Sec. III, the Lagrangian L contains terms which are present in ρ φ and P φ . We first derive the form of L consistent with the condition L ∝ H 2 and study whether some additional conditions are required to satisfy the scaling relations of each term in ρ φ and P φ . Then, the Lagrangian obeyṡ which corresponds to p = 1 in Eq. (3.7). On using Eq. (4.1), it follows thaṫ where , (4.8) . (4.9) Substituting the Lagrangian (2.2) into Eq. (3.7), it follows that the couplings G = G 2 , G 3 , B 4 , A 4 need to separately obey the partial differential equations: (4.11) The integrated solution to Eq. (4.10) is generally given by whereg is an arbitrary function of 13) and ψ is defined by (4.14) From Eq. (4.12), each coupling is restricted to be , (4.18) whereg 2 ,g 3 ,b 4 ,ã 4 are arbitrary functions ofỸ , and The Lagrangians G 2 and G 3 agree with those derived in Refs. [58] and [63], respectively. The couplings A 4 and B 4 are related to each other according to the degeneracy condition (2.4). On using Eqs. (4.17) and (4.18), it follows that Q ,φ /Q 2 = constant. This is integrated to give where µ 1 and µ 2 are constants. Thus, the degeneracy condition restricts the coupling to be of the form (4.21). In this case, both q 1 (φ) and q 2 (φ) are proportional to Q 2 (φ). Absorbing the proportionality constant into the definitions ofg 3 (Ỹ ) andã 4 (Ỹ ), the Lagrangian corresponding to the functions (4.15)-(4.18) yields where, from the degeneracy condition (2.4), the functioñ a 4 (Ỹ ) is constrained to bẽ a 4 (Ỹ ) = 3Ỹb 2 4,Ỹ (Ỹ ) 2 1 + 2b 4 (Ỹ ) . (4.23) From Eq. (2.15), we have which depends onỸ alone. On using Eq. (4.1), we find that the quantity ψ defined by Eq. (4.14) has the dependence ψ = (2h/λ) ln a + ψ 0 , where ψ 0 is a constant. Then, we have eλ ψ ∝ a 2h ∝ H −2 and hencẽ Y ∝ Q 2φ2 H −2 = constant. This means that the couplings B 4 and B 1 do not change in time in the scaling regime. Exploiting Eq. (4.1) together with the property Q ,φ /Q 2 = −µ 1 = constant, it follows thaẗ Moreover, the φ and X derivatives of B 4 =b 4 (Ỹ ) have the dependence where the function B 1 also satisfies similar relations. Then, all the terms appearing in ρ φ and P φ are in proportion to H 2 . We have thus shown that the Lagrangian (4.22) with the coupling (4.21) ensures the existence of scaling solutions. For the Lagrangian G 2 = Xg 2 (Ỹ )Q 2 (φ), there exist scaling solutions for any arbitrary coupling Q(φ). In theories containing the functions G 3 and B 4 , however, the coupling is constrained to be of the form (4.21). In cubic Horndeski theories the degeneracy conditions are absent, so the coupling Q(φ) does not seem to be constrained. Substituting Eq. (4.16) into the G 3 -dependent terms of ρ φ and P φ , however, they are proportional to H 2 only for Q ,φ /Q 2 = constant. Hence we obtain the same coupling as Eq. (4.21) [63].

B. Constant Q
If µ 1 = 0, then the matter coupling (4.21) is constant (Q = 1/µ 2 ). Since ψ = Qφ, the function (4.17) yields We derive the fixed points for the dynamical system given by the Lagrangian (4.28) with the application to dark energy in mind. The background Eqs. (2.11) and (2.12) contain the time derivativesḢ and ... φ , but they can be eliminated to give where f 1 , f 2 , f 3 are functions of their arguments. The branch with an expanding Universe corresponds to The quantity f 2 2 − 4f 1 f 3 does not possess the second derivativeφ, whereas f 2 /(2f 1 ) contains the term proportional toφ. Then, the Hubble parameter can be expressed in the form where P m is related to ρ m according to w m = P m /ρ m = constant. Taking the time derivative of Eq. (4.30) and using the continuity Eq. (2.17),Ḣ and ... φ appear again. However, they can be eliminated by using Eq. (2.11). The resulting equation can be combined with Eq. (4.30) to solve for the second-order field derivativeφ in the form The above discussion shows that the dynamical system is kept up to second order in time derivatives for both φ and a.
To derive the fixed points of DHOST theories given by the Lagrangian (4.28), it is convenient to introduce the following dimensionless variables: where the quantity Y can be expressed as Y = x 2 /y 2 .
Since both x and Y are constant along the scaling solution, y does not vary in time either. The variables x and y obey where ǫ φ =φ/(Hφ) and h = −Ḣ/H 2 , and a prime represents a derivative with respect to N = ln a. The quantities ǫ φ and h are known from Eqs. (4.33) and (4.34).
The fixed points of the dynamical system (4.36)-(4.37) can be derived by setting x ′ = 0 and y ′ = 0. The scaling solution obtained for constant Q corresponds to where the subscript "c" represents the value on the critical point. Since the relationsφ = − √ 6H 2 hx c and ... φ = 2 √ 6H 3 h 2 x c hold on the fixed points, we substitute them into Eqs. (2.11)-(2.12) and solve them for g 2,Y and b 4 , respectively. On using Eq. (4.33), the fixed points satisfying the condition (4.38) obey There are the following two fixed points.
• (a) Scaling solution: This corresponds to the case in which Ω φ and Ω m are non-vanishing constants. Along this solution, w φ and w eff are given, respectively, by For the vanishing coupling (Q = 0), it follows that w φ = w eff = w m . In the presence of the coupling Q, the scaling solution can lead to the cosmic acceleration for w eff < −1/3, but we need |Q| to be larger than the order |λ| for achieving this purpose.
There is another fixed point satisfying Ω φ = 1. In this case, we have where x c is known for given functions g 2 (Y ), g 3 (Y ), and b 4 (Y ). If λx c < √ 6/3, then the point (b) can be used for the late-time cosmic acceleration. For the dynamical system (4.36)-(4.37), there exist other kinetic-type fixed points satisfying The functions g 2 (Y ), g 3 (Y ), b 4 (Y ) consistent with the background equations of motion are given by where c n , d n , e n are constants and n is an integer. In g 3 (Y ), we do not include a constant d 0 since it is just a total derivative. The termd 1 ln Y can be taken into account in g 3 (Y ) as in Refs. [61,63]. Here, we do not do so since we are interested in the effect of the function b 4 (Y ) on the fixed points. A constant e 0 is not included in b 4 (Y ) by reflecting the fact that this is merely a shift of the reduced Planck mass. We substitute Eqs. (4.44)-(4.46) and their Y derivatives into Eqs. (2.11)-(2.12) and solve them for Ω m and h. Plugging these relations into Eq. (4.33), we find that there are the following two fixed points.
This is characterized by with (4.50) The constant e n (n ≥ 1) in b 4 (Y ) does not modify the values of w φ , w eff , and Ω φ of the standard φMDE [69]. For w m = 0, we have w eff = Ω φ = 2Q 2 /(3c 0 ). Provided that |Q| ≪ 1, the φMDE can replace the standard matter era.
There exists another kinetic point satisfying This can be used for neither radiation/matter eras nor the cosmic acceleration.
In cubic-order Horndeski theories, it was shown in Ref. [63] that there exist viable dark energy models with the φMDE followed by the fixed point (b). In the presence of the coupling b 4 (Y ) of the form (4.46), it is of interest to study in detail how the cosmological dynamics and the evolution of perturbations are subject to change compared to cubic-order Horndeski theories.

V. CONCLUSIONS
In this paper, we considered quadratic-order DHOST theories satisfying degeneracy conditions to avoid the Ostrogradsky instability, the constraint on the speed of gravitational waves, and the bound on the decay of gravitational waves to dark energy perturbations. The Lagrangian of this class is given by Eq. (2.2), where A 4 is related to B 4 according to Eq. (2.4). We derived the most general Lagrangians that are able to reproduce separately tracking and scaling behaviors under the condition that h = −Ḣ/H 2 is approximately constant. In the absence of coupling Q between the scalar field and matter, we obtained the Lagrangian of tracking solutions satisfying the conditions L ∝ H 2p andφ ∝ H p . In particular, the scaling behavior corresponds to the choice p = 1.
In Sec. III A, we showed that, for Q = 0, the exact tracker solution exists up to the cubic-order Horndeski Lagrangian with the functions G 2 = c 2 X and G 3 = c 3 X (p−1)/(2p) . We verified that these contributions to the background equations obey the relations ρ φ ∝ P φ ∝ H 2p . In the presence of the DHOST Lagrangian, we found that the function B 4 = c 4 X (p−1)/p leads to the approximate tracker solution at early times when the terms proportional to H 2p (with p < 1) are the dominant contributions to ρ φ and P φ . At late times, the other terms in ρ φ and P φ , which grow faster than H 2p , give rise to a variation in the field equation of state w φ . For p = 1, we found that the exact scaling solution can be realized by the DHOST Lagrangian given by Eqs. (3.32)-(3.35) with (3.36).
In Sec. IV, we extended the analysis of scaling solutions to the case of a field-dependent coupling Q(φ). The most general Lagrangian with scaling solutions is of the form (4.22), withã 4 (Ỹ ) related tob 4 (Ỹ ) according to Eq. (4.23). We showed that the degeneracy condition (2.4) fixes the form of the coupling to be Q(φ) = 1/(µ 1 φ + µ 2 ), including the constant Q as a special case. Indeed, we verified that all the terms in ρ φ and P φ are in proportion to H 2 . The coupling Q(φ) can be arbitrary for the quadratic Lagrangian L = Xg 2 (Ỹ )Q 2 (φ) alone, but the existence of cubic and quartic Lagrangians restricts the coupling to be of the above form to satisfy the scaling property of each term in ρ φ and P φ (as shown in Ref. [63] for the cubic Lagrangian).
For a non-vanishing constant Q, the Lagrangian with scaling solutions reduces to the form (4.28) with (4.29), which matches with the result found for Q = 0 in Sec. III B. In Sec. IV C, we derived the fixed points of the dynamical system described by this Lagrangian. In particular, we obtained four fixed points: (a) a scaling critical point, (b) a scalar-field dominated point, (c) a φMDE point, and (d) a purely kinetic critical point. The points (c) and (d) arise for the models given by the functions (4.44)-(4.46). The point (a) is unlikely to be responsible for the late-time cosmic acceleration with w eff close to −1 because one would need a large value for the coupling |Q|, while the observations of temperature anisotropies in cosmic microwave background place the upper bound |Q| < O(0.1) [70]. On the other hand, the other scaling point (c) can replace the standard matter era. Moreover, the point (b) can be used for driving the cosmic acceleration.
It would be of interest to apply the Lagrangians with tracking and scaling solutions to the construction of concrete dark energy models. In particular, one can investigate whether there exists a viable cosmology allowing for the φMDE point (c) followed by the accelerated point (b) without ghost and Laplacian instabilities. In such a case, one can explore the differences with the dark energy model in cubic-order Horndeski theories where a viable cosmological sequence exists [63]. The analysis of cosmological perturbations is also important to compare those models with the observations associated with the cosmic growth history.