Limits on Clustering and Smooth Quintessence from the EFTofLSS

We apply the Effective Field Theory of Large-Scale Structure (EFTofLSS) to analyze cosmological models with clustering quintessence, which allows us to consistently describe the parameter region in which the quintessence equation of state $w<- 1$. First, we extend the description of biased tracers in redshift space to the presence of clustering quintessence, and compute the one-loop power spectrum. We solve the EFTofLSS equations using the exact time dependence, which is relevant to obtain unbiased constraints. Then, fitting the full shape of BOSS pre-reconstructed power spectrum measurements, the BOSS post-reconstruction BAO measurements, BAO measurements from 6DF/MGS and eBOSS, the Supernovae from Pantheon, and a prior from BBN, we bound the clustering quintessence equation of state parameter $w=-1.011_{-0.048}^{+0.053}$ at $68\%$ C.L.. Further combining with Planck, we obtain $w=-1.028_{-0.030}^{+0.037}$ at $68\%$ C.L.. We also obtain constraints on smooth quintessence, in the physical regime $w \geq -1$: combining all datasets, we get $-1\leq w<- 0.979$ at $68\%$ C.L.. These results strongly support a cosmological constant.


Introduction and Summary
Introduction The analysis of the Full Shape (FS) of the BOSS galaxy power spectrum with the Effective Field Theory of Large-scale Structure (EFTofLSS) at one loop has provided us with a measurement of all parameters in ΛCDM with just a Big Bang Nucleosynthesis (BBN) prior [1,2,3] (see also [4] for other prior choices and [1] for a joint analysis with the BOSS bispectrum using the tree-level prediction). The combination with BOSS reconstructed measurements and baryon acoustic oscillations (BAO) from eBOSS, as well as with supernovae redshift-distance or cosmic microwave background (CMB) measurements, has further allowed us to bound the total neutrino mass, and put limits on the effective number of relativistic species, on smooth dark energy, or on curvature [1,3,5,6,7,8]. In particular, the FS analysis can help constrain models invented to address the Hubble tension as it provides measurements independent on the CMB or local distance ladders [9,10,11,12]. All these results were made possible thanks to the development of the EFTofLSS, which is a powerful tool to extract cosmological information from Large-Scale Structure surveys. A long line of study was necessary to bring the framework to the level where it can be applied to the data. We therefore find fair to add the following footnote where we acknowledge a fraction of its important developments, though not all intermediate results are used in the present analysis 1 . 1 The initial formulation of the EFTofLSS was performed in Eulerian space in [13,14], and subsequently extended to Lagrangian space in [15]. The dark matter power spectrum has been computed at one-, twoand three-loop orders in [14,16,17,18,19,20,21,22,23,24,25]. These calculations were accompanied by some theoretical developments of the EFTofLSS, such as a careful understanding of renormalization [14,26,27] (including rather-subtle aspects such as lattice-running [14] and a better understanding of the velocity field [16,28]), of several ways for extracting the value of the counterterms from simulations [14,29], and of the non-locality in time of the EFTofLSS [16,18,30]. These theoretical explorations also include enlightening In this paper, we analyze the BOSS FS power spectrum using the EFTofLSS at one loop in the context of clustering quintessence [64,65,66] and smooth quintessence. In clustering quintessence, dark energy is made of a scalar field (the quintessence field) whose fluctuations have effectively zero speed of sound, c s , and therefore 'cluster', as they can fall into gravitational potentials. It is a particularly appealing model since the dark energy equation of state parameter w can cross the so-called phantom divide, w = −1 and consistently describe the regime w < −1. This is allowed thanks to the presence of higher-derivative operators in the Lagrangian that stabilize gradient instabilities, but this can only happen if c 2 s → 0 such that they remain not parametrically suppressed. Clustering quintessence has been considered within the context of structure formation in [67,68] and in the EFTofLSS in [53] (see also [24,54,55] for embeddings of other dark energy theories in the EFTofLSS). In this work, we extend the description to biased tracers in redshift space with exact-time dependence in order to apply it to data from galaxy surveys. We remark that we find it quantitatively important to solve the EFTofLSS equations with the exact time dependence, rather than with the approximate, so-called 'EdS', approximation. As for smooth quintessence, which has already been analyzed in light of the BOSS FS and LSS data in [7], we will perform here the analysis by imposing a physical flat prior −1 ≤ w on the smooth quintessence equation of state parameter. By wCDM, we refer to a Universe that includes a smooth dark energy component, i.e. a scalar quintessence field with c 2 s → 1, whose perturbations can be neglected since the sound horizon is of the size of the cosmological horizon. In this picture, w < −1 is an unphysical region where the vacuum is unstable, therefore we should analyze wCDM excluding this region (see discussions in e.g. [69,64]). This paper is organized as follows. We compute the power spectrum at one loop in redshift space for biased tracers with exact time dependence for the clustering quintessence model in Section 2. Further details concerning this derivation are given in the appendices. In Section 3, we apply our framework to LSS data.

Data sets
We analyze the FS of BOSS DR12 pre-reconstructed power spectrum measurements [70], baryon acoustic oscillations (BAO) of BOSS DR12 post-reconstructed power spectrum measurements [71], 6DF [72] and SDSS DR7 MGS [73], as well as high redshift Lyman-α forest auto-correlation and cross-correlation with quasars from eBOSS DR14 measurements [74,75]. We also consider combinations with Supernovae (SN) measurements from the Pantheon sample [76] and with Planck2018 TT,TE,EE+lowE+lensing [77]. studies in 1+1 dimensions [29,31]. An IR-resummation of the long displacement fields had to be performed in order to reproduce the Baryon Acoustic Oscillation (BAO) peak, giving rise to the so-called IR-Resummed EFTofLSS [32,33,34,35,36]. An account of baryonic effects was presented in [37,38]. The dark-matter bispectrum has been computed at one-loop in [39,40], the one-loop trispectrum in [41], and the displacement field in [42]. The lensing power spectrum has been computed at two loops in [43]. Biased tracers, such as halos and galaxies, have been studied in the context of the EFTofLSS in [30,44,45,46,47,48] (see also [49]), the halo and matter power spectra and bispectra (including all cross correlations) in [30,45]. Redshift space distortions have been developed in [32,50,47]. Neutrinos have been included in the EFTofLSS in [51,52], clustering dark energy in [53,24,54,55], and primordial non-Gaussianities in [45,56,57,58,50,59]. The exact-time dependence in the loop has been clarified in [60,61]. Faster evaluation schemes for the calculation of some of the loop integrals have been developed in [62]. Comparison with high-fidelity N -body simulations to show that the EFTofLSS can accurately recover the cosmological parameters have been performed in [1,3,63]. Methodology We analyze the BOSS FS using the galaxy power spectrum in redshift space at one loop in the EFTofLSS [47] following the methodology described in [1,3]. The description of the likelihood, including the covariances and priors used, can be found in [1]. The theory of biased tracers in redshift space with exact time dependence in clustering quintessence cosmology at one loop is derived in Section 2 (see also [61] which has already derived the same expressions, but just in real space, with a different approach), and the scale cut up to which the FS is analyzed is discussed in Sec. 3.1. The power spectrum is IR-resummed [32,34,35,7], and includes corrections to observational systematics: the Alcock-Paczynski effect [78], window functions [79], and fiber collisions [80]. We sample over the following cosmological parameters: the abundance of baryons ω b , the abundance of cold dark matter ω cdm , the Hubble constant H 0 , the amplitude of the primordial fluctuations ln(10 10 A s ), the tilt of the primordial power spectrum n s , and the quintessence equation of state parameter w. We impose no prior on the cosmological parameters but a BBN prior on ω b : a Gaussian prior centered on 0.02235 with σ BBN = 0.0005, obtained by adding up the theory and statistical errors of [81]. We use the Planck prescription of one single massive neutrino with mass 0.06 eV as done in [77]. Allowing the EFT parameters to vary only within physical ranges, we impose priors on them as in [7].
The BAO measurements from the post-reconstructed BOSS power spectrum are correlated with BOSS pre-reconstructed (FS) measurements. The joint analysis is described in [7] (see also [6]). When adding BAO from 6DF/MGS or eBOSS, SN from Pantheon, or Planck data, we simply add the loglikelihoods as these measurements are uncorrelated among each other. We neglect the small crosscorrelation between LSS data with Planck weak lensing and the integrated Sachs-Wolfe (ISW) effect.  84 , and also obtain ln(10 10 A s ) = 3.046 +0.014 −0.014 and n s = 0.9665 +0.0042 −0.0036 at 68% C.L.. All analyses performed here show that our Universe is consistent with ΛCDM. First, clustering quintessence in the limit w = −1 reduces to ΛCDM, and we find that w is consistent with −1 at 68% C.L. Second, the values obtained for the other cosmological parameters in clustering quintessence are consistent within 68% C.L. with the ΛCDM ones obtained by fitting BOSS FS with the EFTofLSS [1,2,3], in combination with other probes [5,6,7], or fitting Planck alone [77]

Main Results
A similar observation applies when fitting wCDM with a flat prior on the dark energy equation of state parameter of w ≥ −1. Fitting BOSS data with a BBN prior, we find in this case Ω m = 0.337 +0.017 −0.022 and H 0 = 68.6 ± 1.8, and we bound −1 ≤ w < −0.91 at 68% C.L. (−1 ≤ w < −0.81 at 95% C.L.). We also get ln(10 10 A s ) = 2.77 ± 0.19 and n s = 0.885 +0.069 −0.058 at 68% C.L.. Adding BAO measurements, Pantheon SN and Planck data we obtain the very stringent constraint −1 ≤ w < −0.979 at 68% C.L. (−1 ≤ w < −0.956 at 95% C.L.). Thus, allowing wCDM only within the physical region gives tight posteriors that are also consistent with the ones obtained on ΛCDM fitting BOSS or Planck. This is illustrated in Fig. 2.
We end this summary of the main results with a note of warning. It should be emphasized that in performing this analysis, as well as the preceding ones using the EFTofLSS by our group [1,3,7,9], we have assumed that the observational data are not affected by any unknown systematic error, such as, for example, line of sight selection effects or undetected foregrounds. In other words, we have simply analyzed the publicly available data for what they were declared to be: the power spectrum of the galaxy density in redshift space. Given the additional cosmological information that the theoretical modeling of the EFTofLSS allows us to exploit in BOSS data, it might be worthwhile to investigate if potential undetected systematic errors might affect our results. We leave an investigation of these issues to future work.  Table 2 in [3]. This plot illustrates the consistency of the datasets as well as the consistency of the present analyses with a cosmological constant.

Biased tracers with exact time dependence in clustering quintessence
In this section, we extend the study of biased tracers in redshift space with exact time dependence, first studied in [60,61], to clustering quintessence.

Review of the EFTofLSS with clustering quintessence
We start by reviewing the underlying equations of motion for dark matter and the dark energy component. For a more detailed discussion, we refer the reader to [53]. In the EFT of dark energy, previously studied in [64,86,65,87], the dark energy degree of freedom is assumed to be the Goldstone boson arising from the spontaneous breaking of time diffeomorphisms. To write the most general theory, we work in unitary gauge where the scalar degree of freedom appears in the metric. The gravitational action will contain operators that break time diffeomorphisms, while remaining invariant under timedependent spatial diffeomorphisms. Up to second order in perturbations, and at leading order in the derivatives, the action reads where we use the 'u' subscript, to emphasize that the metric in the action above is in unitary gauge.
Here δK ij is the perturbation of the extrinsic curvature tensor, and δK is its trace. For simplicity, in the following we work withm 1 = 0, but it can be checked [64,65] that this operator describes a clustering quintessence at cosmological scales. The operators proportional tom i are negligible on large scales as they scale as ∼ k 4 , but are necessary to guarantee the stability of perturbations, as discussed below. To S G , we add the action for matter S M , which we take to be fully diffeomorphism invariant. This guarantees that, once we explicitly reintroduce the Goldstone mode π, there will be no direct couplings of π to matter. The background equations we obtain from S G + S M are the familiar Friedmann equations: where we set the cold dark matter pressure p m = 0, and define the background dark energy density and pressure by From the Friedmann equations we obtain the background solutions for the dark matter and dark energy densities: where the sub index 0 stands for the present day value, and we use the equation of state parameter for dark energy w = p D /ρ D . In the following, we will often use the present day fractional densities Ω x,0 = ρ x,0 ρ D,0 +ρ m,0 , with x ∈ {m, D}. Starting from the action in unitary gauge, it is useful to explicitly reintroduce the Goldstone mode doing the Stueckelberg trick. We perform the time diffeomorphism x 0 → x 0 + ξ 0 ( x, t) and x i → x i , and then substitute ξ 0 (x) → −π(x). The replacement rules for the coefficients and the metric are (for details see for example [53]) g 00 u → g 00 + 2g 0µ ∂ µ π + g µν ∂ µ π∂ ν π .
Gravitational perturbations will be described by the spatially flat perturbed FLRW metric in Newtonian gauge: where Φ and Ψ are the gravitational potentials, and we ignore tensor fluctuations. We then obtain the action for the Goldstone boson π up to second order: At short distances, one can focus on the action of the Goldstone boson. We can see that the kinetic part is given by and thus the speed of sound is The theory must be free of ghosts, which implies that the denominator has to be positive. Therefore the speed of sound needs to have the same sign as 1 + w. In particular, w < −1 implies c 2 s < 0, which would produce gradient instabilities. One can circumvent this instability by including the higher derivative terms proportional tom 1,2 , which scale as k 4 and give a stable dispersion relation at small scales [64,65]. In order for the higher derivative terms not to be highly suppressed (which would make them irrelevant on cosmological scales), we need the speed of sound to be bound by |c 2 s | < 10 −30 , which means it is practically zero. These considerations hold also when a careful analysis including the mixing with gravity is performed. Similar considerations are obtained by including the higher derivative operator proportional tom 1 [64,65]. In conclusion, it is possible to have viable theories with w < −1, but they need to have c 2 s → 0, which are called clustering dark energy or clustering quintessence. We notice furthermore that in order to have a stable theory, we need to have w −2 if we use the operators inm 1,2 , or ≥ −1.17 if we use the operator inm 1 [64,65] The name stems from the fact that the dark energy can cluster with the dark matter, and they jointly contribute to the gravitational potential. Hence the adiabatic mode (i.e. the perturbations of the total energy density, which source the gravitational potential) depends on both the dark matter and dark energy perturbations. As a result, dark energy perturbations leave an imprint on biased tracers such as galaxies, which are the main interest in this work. Therefore, next we wish to give a quick overview of how we derive the equations of motion for the adiabatic mode in the presence of clustering quintessence.
Before analyzing the equations for π, it is useful to write down the EFT equations for dark matter, which couples to dark energy through gravity [53]: where δ m and v m are the dark matter overdensity and velocity,˙= d/dt and τ ij is the effective stress tensor.
Let us start analyzing the linear equations, and we will study the non-linear equations subsequently. The linear equation for π [53,65,66], which we get from (10), reads: This shows that, in the limit c s → 0, the RHS can be neglected. We can, therefore, writeπ − Φ ∝ (a 3 M 4 2 ) −1 , which is a decaying mode, assuming the speed of sound to be approximately constant. In particular, we have ∂ iπ − ∂ i Φ = 0, and, using the linear-level Euler equation (14), we get that d dt av i m + ∂ i π = 0. This means that on the growing adiabatic mode we have which implies that the two species are comoving. This will eventually allow us to write a closed set of differential equations for the adiabatic mode, defined by Using the definition of the adiabatic mode, we find We can now take the derivative of the above equation and plug in the equation of motion for π, Eq. (15), the solution for ρ D , Eq. (6), and substitute the dark matter velocity for the spatial derivatives of π, Eq. (16). We then getδ where we have introduced the dark matter velocity divergence θ m = ∂ i v i m and we have defined We now move on to the full non-linear equations of motion for the adiabatic mode, which is somewhat more technical. We will just mention the main results and refer to [53] (see also [88]) for more details. First, we can easily see that the two species remain comoving at non-linear level. Using the equations of motion, one can show that δg 00 u ∝ c 2 s also at non-linear level. Taking a spatial derivative, ∂ i δg 00 u = 0 in the limit c 2 s → 0, yields This is satisfied by simply using Eq. (16), thus the two species are comoving also at non-linear level. The full non-relativistic equation of motion for the dark energy field π is given by where we used that ∂ i δg 00 u = 0. The full Poisson equation introduces non-linearities in the definition of the adiabatic mode, which reads Now we can take a time derivative and obtain a non-linear continuity equation for the adiabatic mode. The only difference is that we have to include the non-linear terms forδ m and we have an additional term in the equations of motion for π on the right-hand side of Eq. (23). We then geṫ where in the second line we use Eq. (24), and in the last line we use Since the two species are comoving, θ A = θ m and the Euler equation for the adiabatic mode is simply obtained by using the definition of the adiabatic mode in terms of the gravitational potential in Eq. (14). We finally get the governing equations for the clustering quintessence -dark matter system (without counterterms):δ where ρ m /(2M 2 Pl ) = 3Ω m,0 H 2 0 a 0 /(2a 3 ). As explained in [53], since clustering quintessence is comoving with dark matter, there is no isocurvature mode, and the counterterms are the same as for standard dark matter. To solve the equations above perturbatively we transform into Fourier space, where they read (still neglecting the counterterms): dependence of a perturbation by powers of the growth factor, for instance δ (a i ), for some intital time a i . Instead, we will use the exact time dependence solution discussed below. As we will see later, the EdS approximation significantly biases the determination of the cosmological parameters in the presence of clustering quintessence.
Eqs. (28)- (29) are slightly different from the dark matter equations in the presence of smooth dark energy with c 2 s = 1, i.e. wCDM. In fact, in the limit (1 + w) → 0, with Ω D,0 =const, we recover, at large distances where we can neglect the higher derivative terms, the equations of motion for the matter overdensity in ΛCDM. This difference in the equations of motion between the two models results in a modified definition of the time functions that appear in the exact time solutions for δ and θ. Exact solutions for the adiabatic mode δ in the presence of clustering quintessence have been previously studied in [53,67,88]. The time-dependent integral kernel solutions in Fourier space are given by [53] K (1) where repeated σ ∈ {1, 2} are summed over and λ ∈ {δ, θ}. The explicit time functions are defined in Appendix A, and the momentum functions in Appendix B. The kernels in Eqs. (30)- (32), and in the following sections are defined by where X may for instance stand for δ or θ. In the next section, we will see how the solution with exact time dependence for clustering quintessence leaves an imprint in the bias expansion of biased tracers such as galaxies.

Perturbative expansions of δ h and θ h
To find the bias expansion for the galaxy overdensity δ h following the exact time dependence solution of the adiabatic mode, we can follow a procedure similar to [60]. Ref. [61] has also recently derived the same results, using a different approach. Here equations will change with respect to [60], as a consequence of the modified equations of motion for δ A , relative to the equations for the dark matter solutions in wCDM. As has been previously studied in [30], the bias expansion for δ h is given by where we include all possible operators allowed by the equivalence principle, including stochastic contributions and higher derivative terms. Their definitions are found in Appendix B. As for the dark matter equations, since clustering quintessence is comoving with dark matter, there is no isocurvature mode, and the bias expansion depends on the same fields as for the dark-matter-only universe [53].
The time-kernels, such as c δ (a, a ), that account for the time non-locality, can be formally integrated over a after the perturbative solutions are substituted in. All operators (which are explicitly given in Appendix B) are evaluated along the fluid line-element: This results in Taylor expansions in the fields around x given by It turns out that even in the presence of clustering quintessence, once we perturbatively expand the overdensity and velocity, the time integrals in Eq. (36) can be done analytically and the solutions are given in terms of the time functions and kernels that appear in Eqs. (30)- (32). This is explicitly derived in Appendix C. Then, as mentioned before, after perturbatively expanding the fields, the time integrals in Eq. (34) are formally done, and result in the definition of coefficients such as For a complete list see Appendix B. After this procedure, the resulting halo overdensity can then be written as a sum of functions of time multiplied by functions of momentum. As was shown in [60], some of the momentum functions are degenerate and can all be expressed in terms of the basis {I, α, β, α 1 , α 2 , β 1 , β 2 , γ 1 , γ 2 }, which are the kernels that appear in Eqs. (30)- (32). This is true in wCDM as well as the clustering quintessence case, because the momentum functions are the same in both cases, and only the time functions change. We can therefore write where in the last expression a sum is implied over σ ∈ {1, 2}. The main reason that the time coefficients c i change, relative to wCDM, is because the integrals from the flow terms that stem from the Taylor expansion of Eq. (36) now have an additional dependence on C(a) (for a comparison see Appendix C). The coefficients in Eq. (38) are explicitly defined in Appendix B. For more details on the derivation of the halo overdensity kernels, see [60]. From here we can proceed in a very similar fashion to [60]. We reduce the number of coefficients by looking for degeneracies in the time coefficients. Luckily all the identities from [60] still hold in a slightly more general form. The main difference here is that we define the calculable function where in the limit G wCDM = 1 we recover the identities from [60]. Y (a) is defined by However, it is useful to defineỸ so that, taking limits, we haveỸ (a) wCDM = Y (a) EdS = 0. We can then write the final halo overdensity (see also [61]): β ( k, a) where we can see that no new C i operators have to be included compared to the exact wCDM case or EdS approximated case. The C i are defined in the same way as in [60] and are explicitly given in Appendix B. Similarly to what happens when we use the exact time dependence for smooth dark energy and ΛCDM, we see that there are additional calculable time dependencies in the final bias expansion for the galaxy overdensity. However, there are no new bias coefficients. We can take two interesting limits to see how the above expansion generalizes previous models. In the G → 1 limit, we obtain the galaxy overdensity in wCDM with exact time dependence. Furthermore, in the limit where we use the EdS approximation, the time functions in Eqs. (30)-(32) become independent of a and with a value so that G → 1 andỸ → 0. Eq. (42) can then simply be linearly transformed into the BoD basis from [45], therefore the space spanned by the kernels in Eq. (42) is the same as the one spanned by the BoD basis from [45] (for a transformation see [60]). For illustration, we plot in Fig. 3 the values ofỸ and G as functions of the redshift z = 1/(1 + a) and w. show ΛCDM and wCDM cases for comparison. Notice that, as we argued earlier, for w < −1 we need c 2 s → 0 and thus for c 2 s = 1, i.e. wCDM, w < −1 is not allowed in the EFT of dark energy. We, nevertheless, plot it here for illustration.
In a last step, we write the expansion for θ h , which appears in the redshift space expansion. For the velocity divergence, there is no bias [47], up to higher derivative terms. We can thus model the velocity divergence as a species of biased tracer. Specifically, we obtain the velocity divergence by plugging in the following choice of functions into Eq. (42): The counterterms will take the exact same form as for wCDM [53,60]. We will now transform into redshift space and compute the power spectrum.

Galaxy Power spectrum in redshift space
As the next step, we wish to compute the full galaxy power spectrum in redshift space, which we will later use to fit the data. As shown in [60], the EdS approximation has no influence on the transformation into redshift space 7 . This means we can proceed in the same way as described in [47]. The galaxy overdensity kernels in redshift space in terms of the real space quantities δ h and θ h are given by (without counterterms) δ h,r ( q 1 , q 2 , q 3 , µ, a) = K δ h ( q 1 , q 2 , q 3 , a) + f + µ 2 123 K θ h ( q 1 , q 2 , q 3 , a) δ h,r ( q 1 , q 2 , µ 123 , a) where δ h,r is the halo overdensity in redshift space. Usingẑ as the line of sight unit vector, we have defined µ = q ·ẑ/q, with q = q 1 + · · · + q n , and µ i 1 ...in = q i 1 ...in ·ẑ/q i 1 ...in , q i 1 ...im = q i 1 + · · · + q im . As we mentioned previously, the counterterms and stochastic terms that come from real and redshift space (see [47,60] for a discussion) do not change in the presence of clustering quintessence. Therefore, the final expression for the galaxy power spectrum in redshift space, including the counterterms, reads P g (k, µ, a) = K (1) δ h,r (µ, a) 2 P 11 (k, a) δ h,r ( q, k − q, µ, a) 2 P 11 (|k − q|, a)P 11 (q, a) (45) δ h,r ( q, − q, k, µ, a)P 11 (q, a)

+ 2K
(1) δ h,r (µ, a)P 11 (k, a) c ct where P 11 (k, a) is the time-dependent linear power spectrum for the adiabatic mode, k m k NL is the comoving wavenumber which controls the bias derivative expansion, and n g is the background galaxy number density. In the first line, we have the linear power spectrum in redshift space. In the second and third line, we have the P 13 and P 22 contributions of the loop and in the fourth and fifth line we have the counterterms and stochastic terms, respectively. Finally, the power spectrum is IR-resummed following [32,34,35,7]: as quintessence is comoving with dark matter, the same equations hold. We then apply corrections to take into account the Alcock-Pacszynski effect [78], window functions [79], and fiber collisions [80]. Notice that, as we argued earlier, for w < −1 we need c 2 s → 0 and thus for c 2 s = 1, i.e. wCDM, w < −1 is not allowed in the EFT of dark energy. We, nevertheless, plot it here for illustration.
In Fig. 4, we show the difference between the one-loop galaxy power spectrum multipoles = 0, 2 evaluated in different cosmologies: ΛCDM, wCDM and clustering quintessence, for w = −0.95 and w = −1.05. We also show the difference between the evaluation with and without the EdS approximation for clustering quintessence. It is apparent that the difference between wCDM and clustering quintessence is important with respect to the BOSS error bars. The difference between the evaluation with and without the EdS approximation for clustering quintessence is clearly important, especially in the monopole. Given how large the differences in the power spectrum are, we expect to see differences at the level of the posteriors of the cosmological and EFT parameters.

LSS data analysis
In this section, after calibrating the scale cut of the theory against simulations, we present the results from fitting clustering and smooth quintessence to the BOSS FS, and its combinations with BAO, SN and CMB measurements.

Tests against simulations
To assess the theory-systematic error of the FS analysis, we fit the power spectrum multipoles measured from large-volume N-body simulations on clustering quintessence with a BBN prior. We consider two independent realizations of the BOSS 'lettered' challenge simulations, which are boxes of side length 2.5 Gpc/h, described in e.g. [1]. The first realization is made of four boxes, labelled A, B, F, and G, populated by four different Halo Occupation Distribution (HOD) models. The second realization, labelled D, is populated by another HOD model. Using one box, we can measure for each cosmological parameter the theory-systematic error as the distance in the 1D posterior of the 1σ region to the truth of the simulation. Therefore, the theory-systematic error is zero if the truth lies within the 1σ region. For A, B, F, and G, which are correlated, we average the posteriors for the cosmological parameters, that we label ABFG. Moreover, we can combine ABFG with D, as they are independent realizations, allowing us to measure the theory error using a volume about 14 times the one of BOSS data. To do so, we combine for each cosmological parameter the 1D posterior of the shift of the mean with respect to the truth, as the product of two Gaussian distributions. The distance of the 1σ region to zero in each resulting 1D posterior gives a measure of the theory-systematic error for the combination ABFG+D. For each cosmological parameter, the error bar obtained on ABFG+D represents the smallest theorysystematic error which we can measure, which is between 0.3 · σ data and 0.5 · σ data , where σ data is the error bar obtained by fitting BOSS data. ω cdm h ln(10 10 A s ) n s w Ω m σ stat |σ sys σ stat |σ sys σ stat |σ sys σ stat |σ sys σ stat |σ sys σ stat |σ sys ABFG 0.007|0.000 0.027|0.000 0.11|0.04 0.044|0.000 0.139|0.000 0.021|0.000 D 0.006|0.000 0.018|0.000 0.11|0.04 0.039|0.000 0.093|0.000 0.014|0.000 ABFG+D 0.005|0.000 0.015|0.000 0.08|0.07 0.029|0.000 0.077|0.000 0.012|0.000 Table 1: 68%-confidence intervals σ stat and theory-systematic errors σ sys obtained fitting clustering quintessence to the lettered challenge simulations with a BBN prior.
In Fig. 5 and Tab. 1, we show the results obtained by fitting the lettered challenge simulations at scale cut k max = 0.23h Mpc −1 . We find for all cosmological parameters zero theory-systematic error, with the exception of ln(10 10 A s ), where we find a marginal theory-systematic error of 0.07, which is ∼ 0.4 · σ data 8 . These results show that we can confidently fit the data up to k max = 0.23h Mpc −1 on our high redshift (z eff = 0.57) sample CMASS. For LOWZ sample at z eff = 0.32, we rescale the scale cut as in [1] and fit up to k max = 0.2h Mpc −1 .

LSS constraints
In Fig. 6 and Tab. 2, we show the results obtained by fitting BOSS FS+BAO, and in combination with BAO measurements from 6DF/MGS and eBOSS, and with Pantheon SN, on clustering quintessence with a BBN prior. We see that all cosmological parameters can be measured. For all analyses performed, w is consistent with −1 at 1σ.   Physical considerations We now discuss why all cosmological parameters can be measured by analyzing the FS using the EFTofLSS, and how the addition of the SN measurements helps to obtain better constraints. Let us start with the contribution from the BAO information. The two angles corresponding to the BAO components perpendicular and parallel to the line of sight are given by: Here r d (z CMB ) is the sound horizon at the end of the baryon-drag epoch z CMB , and D A (z LSS ) and H(z LSS ) are the angular diameter distance and the Hubble parameter at the effective redshift of the survey z LSS . As discussed in e.g. [1,7], these angles carry information about h, Ω m and w. The dependence on parameters is the same as in wCDM, as the angles only depend on the background geometry [7]: where z Lyα = 2.35, z CMASS = 0.57, z LOWZ = 0.32 and z 6dF/MGS = 0.106. θ LSS, V is a combination of θ LSS, ⊥ and θ LSS, (see e.g. [7]). The dependences on the cosmological parameters above and in the rest of this section are obtained expanding around a fiducial cosmology (Ω m = 0.3, h = 0.7, w = −1). Furthermore, the relative amplitude of the BAO wiggles with respect to the smooth part instead gives a measurement of ∼ Ω m h 2 (though the amplitude is not part of the standard BAO analysis). Clearly, at least in principle, this information allows for a determination of w, Ω m and h. Notice however that the measurements for w and Ω m are strongly degenerate when using solely the BAO information from CMASS and LOWZ, and the breaking of the degeneracy by measuring both θ LSS,⊥ and θ LSS, is mild, insufficient to get strong constraints [7]. Of course, the situation is greatly ameliorated by the addition of the information from 6dF/MGS and eBOSS, but it is also ameliorated by the inclusion of the FS analysis.
In fact, the FS contains information not only through the BAO signal, but also by its shape and amplitude [1]. The shape depends on the equality scale, and therefore on Ω m h 2 . The amplitude and the anisotropy of the FS can be roughly summarized by the fact that the monopole and quadrupole mainly depend on the combinations b 1 (z) 2 is the amplitude of the linearly evolved power spectrum at the maximum wavenumber of our analysis, A (kmax) s ∼ (k/k 0 ) ns−1 (k eq /k max ) 2 A s , with k eq being the wavenumber that re-enters the horizon at equality and k 0 the pivot scale. D + and f + are respectively the growth factor and growth rate of the growing adiabatic mode. k max is the maximum wavenumber of our analysis, which is where the signal to noise is dominated. Given that there are two redshifts in BOSS, this clearly offers a way to measure both A s and n s , together with b 1 (z CMASS ) and b 1 (z LOWZ ). In this way, all cosmological parameters are, at least in principle, measured. However, we should keep in mind that the FS offers an independent measurement for each wavenumber, therefore, by combining the information from several k's, further information on w and Ω m is obtained. In fact, just by looking at the dependence at linear level of the monopole and quadrupole at z CMASS and z LOWZ , one can see that on top of b 1 and A s , one can measure the combination f (z CMASS )D(z CMASS ) , which, around the fiducial cosmology, goes as ∼ Ω −0.12 m |w| 0.44 . This can be seen by using the fitting functions for D + and f + as a function of Ω m and w given in [67], which read: where C(a) = 1 + (1 + w)Ω D (a)/Ω m (a). This is to be contrasted with the same ratio for the case of a smooth dark energy component, namely wCDM, around the same cosmology: ∼ Ω −0.12 m |w| 0.006 . We can see that the change in the dependence on w going from LOWZ to CMASS is stronger in the case of clustering quintessence compared to wCDM, physically originating from the fact that clustering quintessence contributes to the clustering. The mild degeneracy present for wCDM between Ω m and w is thus less pronounced in clustering quintessence when jointly fitting LOWZ and CMASS. Furthermore, these measurements give different correlations between Ω m and w with respect to the ones in θ LSS , thus further breaking the degeneracies. This can be seen in Fig. 7, where we compare the posteriors obtained fitting BOSS FS+BAO on clustering quintessence and wCDM. To summarize, Ω m , h, w, A s , n s and b 1 can be measured from the BAO angles in combination with the broadband signal. By looking at the same Fig. 7, one can also see that in wCDM there is a large degeneracy in lowering w and lowering A s . This can be explained by the fact that, in wCDM with w < −1 (which, we remind, is physically inconsistent at the quantum level but can still be analyzed as a model), matter domination lasts longer, so that structures grow more and therefore the power spectrum is left unchanged by lowering A s . In clustering quintessence, this degeneracy is broken by the fact that the adiabatic mode receives a contribution from clustering quintessence proportional to 1 + w. This can be see from solving the linear equations, which, at early times, give (see e.g. [53], eq. (4.15)): with a 0 the present epoch and a early a time early on during matter domination. This effect acts in a direction contrary to the extra growth that one gets from the extension of the epoch of matter domination for 1 + w < 0, in practice bounding the degeneracy between w and A s . Note that this discussion gives only rough estimates of the parameter dependence of the FS. In practice, there is no separation between the broadband and the other sources of information within the FS analysis as all the signal is analyzed up to the chosen scale cut. In particular, the loop provides additional information. For example, the growth function enters as D 4 + in the loop, providing yet another parametric dependence on w. In Fig. 7, we also show the posteriors obtained on clustering quintessence with the EdS approximation. The difference with the posteriors obtained with exact time dependence is clearly visible: most notably, about 0.2 σ for H 0 and Ω m , and 0.3 σ for w. At the level of the power spectrum in Fig. 4, the difference is somewhat larger in terms of error bars, but we should remember that in that figure the EFT parameters are fixed. In particular, the large deviation that can be seen in the monopole of Fig. 4 can be partially absorbed below the error bars with a small offset in the shot noise c ,0 /n g of ∼ 0.1. The difference we see between the EdS evaluation and the exact-time one can be traced to the time functions, as for example G 2 , in some loop terms when evaluated with exact time dependence: G(z LOWZ ) 2 ∼ |w| 0.42 and G(z CMASS ) 2 ∼ |w| 0.27 . Because of this, the EdS approximation leads to significant shifts in the posteriors for clustering quintessence.
Finally, the distance-redshift relation of SN data from Pantheon brings evidently more constraints. Approximately, the line degeneracy of the luminosity distance D L = (1 + z) 2 D A is D L (z = 0.25) = Ω −0.05 m |w| 0.1 , which further helps break the degeneracy between Ω m and w when fitting jointly with the FS and BAO.

CMB+LSS constraints
In Fig. 8 and Tab. 3, we show the results obtained fitting clustering quintessence with Planck data in combination with BOSS FS+BAO, BAO measurements from 6DF/MGS and eBOSS and with Pantheon SN.
As expected and apparent from the posteriors, we can see that Planck gives precise measurements on ω b , ω cdm , ln(10 10 A s ) and n s , while constraints on H 0 or Ω m are obtained by the combination with late-time probes, that break the degeneracy in the H 0 − Ω m plane present in the CMB. As discussed in the previous subsection, w is mainly measured thanks to low-redshift measurements. However, the constraints on w are better when adding Planck since the precise measurements of the other cosmological parameters by Planck helps to further break the degeneracies.

wCDM with w ≥ −1
From an effective field theory point of view, there is no known theory, at least to us, that can realize w < −1 with c 2 s → 1. As discussed in previous sections, such theory has a negative kinetic term. For a theory with no Lorentz-violating UV cutoff, the scalar perturbations are unstable, and the vacuum decays into gravitons at an infinite rate [69]. Therefore, w < −1 would either need some other, physical, motivation or one can posit that w is not allowed to be smaller than −1 in wCDM. By doing so, we get the results depicted in Fig. 7 obtained by fitting BOSS data on wCDM with a BBN prior and a flat prior w ≥ −1. We see that the results differ substantially from the ones obtained without a prior on w. In particular, the degeneracy line w − H 0 , open when allowing w to vary below −1, can not be exploited to lift H 0 to higher values than the one found in ΛCDM analyzing CMB or LSS data.
In Fig. 9 and Table 4     where Θ(a −ã) is the Heaviside step function, W (ã) is the Wronskian of D + and D − : and we impose the boundary conditions At second order, the resulting time-dependent functions are given by for σ = 1, 2. At third order we have The degeneracies pointed out in (39) result from the following identities: where again σ ∈ {1, 2}. One can derive these relations using (65)- (70) and the fact that Furthermore, for the derivation of some of the flow terms in Appendix C it is important to use the following relations:
First, we expand the overdensity and velocity divergence perturbatively. Apart from δ (2) , the only second-order term is in the first line, which is given by Next, we take this same term with δ at second order and v at first order. This gives where the expression for clustering quintessence takes the same form as for wCDM, and we used (74).
In the second and third lines of (94) we can take all fields at linear order. We have a da a c δ (a, a ) 1 2 ∂ i ∂ j δ(x, a ) For completeness, the flow terms from δ 2 and s 2 read lm ∂ i (s lm ) (1) ∂ i ∂ 2 θ (1) ] k (a) .