Cosmological Constraints on Horndeski Gravity in Light of GW170817

The discovery of the electromagnetic counterpart to GW170817 severely constrains the tensor mode propagation speed, eliminating a large model space of Horndeski theory. We use the cosmic microwave background data from Planck and the joint analysis of the BICEP2/Keck Array and Planck, galaxy clustering data from the SDSS LRG survey, BOSS baryon acoustic oscillation data, and redshift space distortion measurements to place constraints on the remaining Horndeski parameters. We evolve the Horndeski parameters as power laws with both the amplitude and power law index free. We find a 95% CL upper bound on the present-day coefficient of the Hubble friction term in the cosmological propagation of gravitational waves is 2.38, whereas General Relativity gives 2 at all times. While an enhanced friction suppresses the amplitude of the reionization bump of the primordial B-mode power spectrum at $\ell<10$, our result limits the suppression to be less than 0.8%. This constraint is primarily due to the scalar integrated Sachs-Wolfe effect in temperature fluctuations at low multipoles.


Introduction
The detection of gravitational waves opens a new window into constraining gravity. In general relativity (GR), the line element for scalar mode perturbations in the Newtonian gauge is given by ds 2 = a 2 (τ ) − (1 + 2Ψ) dτ 2 + (1 − 2Φ) dx 2 . The line element for tensor mode perturbations is ds 2 = a 2 (τ ) −dτ 2 + (δ ij + h ij ) dx i dx j , where h ij = ±h + , h × are small perturbations and are the divergenceless, traceless component of the metric. The linearized Einstein equation without a source takes the form of the wave equation h ij = 0, where is the D'Alembertian. Then, going to Fourier space, the gravitational waves can be described byḧ ij + 2ȧ aḣ ij + k 2 h ij = 0, where k is the wavenumber. Dots throughout the paper denote derivatives with respect to conformal time.
Changing the gravitational theory modifies the propagation of gravitational waves. In Horndeski theory, the most general tensor-scalar theory in which the equations of motion are second order [1], the tensor mode propagation equation becomes h ij + [2 + α M (a)] Hḣ ij + c 2 T (a) k 2 h ij = 0, (1.1) using the parameterization developed by [2] and where H =ȧ/a. The Planck mass run rate α M describes how the Planck mass evolves over time and contributes to the gravitational wave friction. It is defined as -1 -

JCAP12(2018)030
where M 2 * is the effective Planck mass. If not constant, α M creates anisotropic stress in the Jordan frame.
The tensor speed c T is given by c 2 T = 1 + α T , where α T is the tensor speed excess that quantifies how much the gravitational wave speed deviates from that of light. Recent observations of GW170817 and its electromagnetic counterpart have placed the stringent bounds This new constraint effectively eliminates all Horndeski theories with α T 0 = 0 and alters the physically allowed values of the other Horndeski parameters [3]. See e.g. [4] for an example of a non-trivial theory compatible with c T = c, [5,6] for discussions of how gravitational wave detections impact Horndeski theory, [7][8][9] for good theoretical discussions about the impact of GW170817, and [10] for observational constraints on Horndeski theory using only the simultaneous detection of GW170817 and GRB170817A. In this paper we use measurements of the cosmic microwave background (CMB), galaxy clustering, baryon acoustic oscillations (BAO), and redshift space distortions (RSD) to constrain the remaining Horndeski parameters in light of this discovery.
The remaining viable tensor mode propagation equation is Two additional parameters, the kineticity α K and braiding α B , complete the parameterized Horndeski theory. Kineticity describes the scalar perturbations' kinetic energy. Large values reduce the scalar sound speed. The kinetic braiding parameter describes how the scalar and metric kinetic terms mix [see e.g. 11-13, for good discussions]. A nonzero value indicates the clustering of dark energy. All four parameters are independent of each other and the background, which we choose to be ΛCDM. ΛCDM GR cosmology is regained when all α i = 0. In section 2 we discuss our parameterization of the Horndeski parameters and what data we use to constrain them, and in section 3 we define the model's stability constraints. We present our results in section 4 and discuss the source of constraining power in section 5. In section 6 we compare our models and discuss our results. We summarize in section 7. In appendix A we provide additional background information about Horndeski's theory and its implementation. Appendix B details how kineticity affects the model's stability constraints and how this influences our results, and in appendix C we discuss our evolution of the Horndeski parameters.

Modeling & data sets
We use EFTCAMB 1 [14,15], a modified version of the Boltzmann code CAMB [16] that does not use a quasi-static approximation, as well as the complimentary CosmoMC version [17], to study the effects of modified gravity on perturbations. See appendix A for further detail on how the parameters are handled in the code. We set R − 1 0.03, where R is the Gelman-Rubin diagnostic [18], as our criterion to obtain convergence in the chains.
To explore the time evolution of the remaining Horndeski parameters and achieve a sufficiently large stable parameter space to perform an MCMC analysis (see appendix B), we 1 http://www.eftcamb.org/codes/download.html, version 2.0.
where M 2 * /m 2 0 = 1 + M in the context of EFTCAMB, and m 2 0 is the Planck mass. M 0 is the fractional deviation of M 2 * from m 2 0 today. α i 0 denotes the parameter value today, and a is the scale factor. This parameterization is similar to that used by the Planck collaboration for the Horndeski parameters [19]. While we do not know the exact functional form of the α i 's, ref. [20] discusses the challenges of parameterizing them as a function of Ω DE and ref. [21] finds that evolving the α i 's as power laws, with both the amplitude and power law index free, is preferred over evolving them with simpler or more complex models.
The authors of [12] parameterized the Horndeski parameters as a function of the dark energy density and found that α K could not be well constrained [see also erratum of 12]. To limit the number of additional degrees of freedom in our analysis, we fix the evolution of α K with α K 0 = 0.1 and κ = 3. We also choose to explore the case of α B 0 = 0 to probe a theory in which the primary modification is due to gravitational waves. Fixing α B 0 = 0 yields a perfect-fluid model that includes anisotropic stress [2]. If both α K 0 = 0 and α B 0 = 0, the scalar propagation speed diverges. For this reason we fix the kineticity at a nonzero value throughout our analysis. See appendix B for a discussion of how kineticity affects the stable parameter space given the imposed stability conditions defined in section 3. We vary the standard cosmological parameters Ω b h 2 , Ω c h 2 , θ, τ , ln(10 10 A s ), n s , as well as the tensorto-scalar ratio r, Planck calibration y cal , the dust power ( = 80, ν = 353 GHz) A B,dust , and the dust frequency scaling parameter β B,dust . We choose not to vary w 0 and w a to minimize the number of free parameters in our analysis and to focus on the propagation of gravitational waves. In table 1 we list our adopted prior cutoffs for the different parameters. All parameters have uniform priors except for y cal and β B,dust , which have the Gaussian priors y cal = 1.0000 ± 0.0025 and β B,dust = 1.59 ± 0.11, respectively.
Given the direct effect the Horndeski parameters have on tensor perturbations [see eq. (1.4) and e.g. [22][23][24], we include in our analysis the B-mode data from the joint analysis of the BICEP2/Keck Array and Planck [25]. Because the Horndeski parameters influence the scalar perturbations, as well, we use the 2015 Planck low-CMB temperature and polarization data, high-temperature data, and the lensing potential measurements [26,27].
It is reasonable to posit the CMB places optimal constraints on the remaining Horndeski parameters at the epoch of decoupling. We studied the optimal pivot scale to measure α B and α M and find that they are uncorrelated at present and are best constrained at late times rather than during recombination. See appendix C for a further discussion. The constraints primarily come from the late integrated Sachs-Wolfe (ISW) effect (see section 5). The pivot scale is model dependent, however, and should be re-examined when the evolution of the Horndeski parameters is better determined. To take advantage of this late time sensitivity, we include the SDSS LRG DR4 [mpk, 28] and BOSS BAO and RSD data sets in our MCMC analysis [29][30][31][32].

Stability constraints
Several viability priors can be set by EFTCAMB to ensure the parameter space yields a stable theory. We enforce the following constraints [see 2, 14, 33, for a full description]: 1. Physical stability: the theory must have a positive effective Newton's constant (i.e., 1 + Ω > 0), and avoid ghost and gradient instabilities. Ghost instabilities arise when the kinetic energy becomes negative, and gradient instabilities occur when the squared speed of sound for perturbations becomes negative [see section 3.3 of 2, 14].
2. Mathematical stability: neither the coefficient ofḧ ij nor the coefficient ofπ may equal 0, ensuring the tensor perturbation and π field equations are well defined, respectively [14,33,34]. In this work we exclude all exponential growth of the π field perturbations, including those due to tachyon instabilities. 2 For details of the mathematical stability conditions, see Equations 40 and 52 in [33], as well as the corresponding π field equation discussion in section IV D and viability condition discussion in section IV F of [33].
3. We require a positive matter density and dark energy density with ω DE ≤ −1/3 at all times.
We do not restrict our analysis to regions of parameter space where c 2 s ≤ 1 or m 2 π ≥ 0. Subluminal sound speeds are required for the theory to be UV complete through standard methods [2,36]. The authors of [37] have shown, however, that Horndeski theories will always have regions of parameter space that include superluminal sound speeds [see also 2, for a brief JCAP12(2018)030 discussion]. If the scalar field couples to matter, its speed and mass are more complicated to compute than if the field were in a vacuum. They cannot be directly read from the π propagation equation because the scalar degree of freedom of the theory is a combination of π and the matter fields. Enforcing c 2 s ≤ 1 in EFTCAMB would not put a limit on the true scalar sound speed in the non-minimal coupling scenario [38]. We note that restrictions on c s can have a severe impact on the stable parameter space. The authors of [39] have shown that restricting the scalar field to propagate subluminally yields a stable parameter space that is a very small subsection of the parameter space allowed when the scalar field is free to have any real sound speed.

Results
In figure 1 and figure 2 we show the resulting posterior probability distributions for the cases with α B 0 = 0 and α B 0 = 0, respectively. CMB denotes the low-TEB, high-TT [26], and BKP [25] data set combination. LSS denotes combining the lensing [27], mpk [28], BAO [29,30], and RSD [31,32] data sets. Constraints on the friction and braiding parameters are quoted in table 2 and table 3. See [40] for constraints on α B from galaxy cluster observations in light of the observation of the electromagnetic counterpart to GW170817.
One may be surprised to see that the bounds on M 0 and α M 0 exclude their GR values at 95% confidence in the case of α B 0 = 0. However, these lower bounds are driven largely by stability constraints on the model when Both ξ, the power law index for the braiding parameter, and β, the power law index related to friction, are relatively unconstrained by the data. In fact, ξ cannot be constrained at the 95% CL for the α B 0 = 0 case with CMB data alone. At 68% confidence, ξ > 2.82 for this case. Note that the nonzero probability for ξ = 0 in figure 1 is due to smoothing artifacts. The constraints on β are almost purely from stability constraints (see appendix B).
There is a small stable region for β near 0 if ξ is near 0, but the data did not prefer this region. The β posterior bounds in figure 1 and figure 2 correspond to the boundaries of the remaining portion of the stable parameter space. ξ has a large stable parameter space, so the best-fit appears driven by the data. However, the plateau in the posterior for large ξ is due to the data being unable to further constrain the parameter. Larger values of β are preferred for larger M 0 . This suggests that, for a power law evolution, the data prefer to minimize the deviation from m 2 0 at early times. The ξ − α B 0 contour in figure 1 shows a similar trend for ξ 2 as α B 0 deviates farther from 0.

Sources of constraints
Where do the constraints on the Horndeski parameters come from? To identify the source of dominating constraining power, we compute power spectra derivatives at our fiducial cosmology for the temperature, E-mode polarization, B-mode polarization, lensing potential, and matter power spectra in figures 3-6. Our fiducial values for the Horndeski parameters,   tensor-to-scalar ratio, and tensor tilt are listed in table 4, and we use the best-fit Planck TT+lowP+lensing+ext ΛCDM parameter values [41]. The fiducial Horndeski values were chosen based on the posteriors in figure 1 and figure 2. Computing the relative difference between the Horndeski and fiducial (rather than GR) spectra ensures all the spectra have a stable theory while only varying a single parameter. The power spectra are most affected at large scales by the Horndeski parameters, and the CMB scalar modes and large scale structure (LSS) are significantly more sensitive to these parameters than the CMB tensor modes.

Scalar perturbations
In figure 3 we show the sensitivity of the unlensed temperature and polarization anisotropies to changes in the Horndeski parameters in the neighborhood of the fiducial parameters, which is a good proxy for the constraining power near the best-fit values. In general, the CMB anisotropies are most sensitive to the Horndeski parameters at 10, though measurements at such scales are limited by cosmic variance. The perturbations are most sensitive to the friction parameter M 0 , and M 0 has a more dramatic effect on the scalar temperature anisotropies than on polarization, lensing, and matter clustering.
The sensitivity of the scalar TT power spectrum to the Horndeski parameters at lowsuggests the ISW effect may be the source of sensitivity. In figure 4 we show power spectra derivatives with and without the ISW effect. Removing the ISW effect erases any sensitivity the scalar TT power spectrum had to the Horndeski parameters, including braiding, indicating the ISW effect is the primary source of constraint for the temperature anistropies. Because a nonzero Planck-mass run rate creates anisotropic stress, the evolution of the Bardeen potentials changes over time [35,42]. The anistropic stress constraint from the spatial traceless component of the Einstein equations makes this clear [2]: where v X = − aδφ φ is the scalar velocity potential, and H = H/a. Friction alters the relationship between the Bardeen potentials, leading to a change in the ISW and, thus, the temperature anisotropies. As seen in figure 4, increasing braiding decreases the ISW effect, reducing the power of the temperature anisotropies at low-[see 42, for a further discussion of α M , α B , and the ISW effect].
Friction's effect on the temperature anisotropies via the ISW also changes the local temperature quadrupole seen by electrons during reionization, altering the low-E-mode. Scattering of the reionization electrons off the quadrupole produces additional polarization at the scale that enters the horizon during reionization. Increasing M 0 makes the reionization bump peakier for the scalar E-modes, boosting the power for < 10 and damping the power near = 10.

Tensor perturbations
Braiding only affects the scalar modes. Increasing the friction term M 0 dampens the tensor perturbations, decreasing the B-mode and tensor E-mode amplitudes (see middle and bottom panels in figure 3). This is expected since the form of eq. (1.4) is that of a damped driven oscillator [see also 24,43]. For the tensor mode polarization, increasing M 0 decreases the  reionization bump. See e.g. [22] and [24] for a further discussion on how friction affects the reionization peak. Only gravitational waves outside the horizon at recombination affect the temperature anisotropies since the gravitational waves decay and oscillate as soon as they enter the horizon. Since friction is a damping term for the gravitational waves, the tensor temperature power spectrum is damped for low-, as well. Nonetheless, friction has a more dramatic effect on the scalar modes, so the scalar temperature and E-mode polarization dominate the CMB constraints.

Large scale structure
Increasing friction damps the lensing potential at all scales, but most significantly at large angular scales (see figure 5). This occurs because |Φ + Ψ| is damped from increasing the friction. The lensing potential is almost as sensitive to the friction parameters as the scalar TT power spectrum at large scales, but present lensing potential measurements do not cover multipoles L < 40 and so the lensing potential does not have strong constraining power.
Both α B and α M alter the growth rate [2]. In the neighborhood of our fiducial cosmology, increasing M 0 boosts clustering for k 10 −2 h Mpc −1 (see figure 6). Increasing α B 0 so that it is closer to GR reduces clustering on similar scales. The matter power spectrum shows weak sensitivity to the Horndeski parameters for the scales directly probed by the LSS measurements. The dominating constraining power on the Horndeski parameters from these measurements comes through their constraints on σ 8 .

Discussion
To compare the two models we can use the Akaike information criterion (AIC) [44], defined as: where χ 2 = χ 2 BKP + χ 2 high− TT + χ 2 low− TEB for the CMB data set combination and χ 2 = χ 2 BKP + χ 2 high− TT + χ 2 low− TEB + χ 2 lens + χ 2 mpk + χ 2 BAO + χ 2 RSD for the CMB+LSS data set combination. Then, where L is the maximum-likelihood and k is the number of fit parameters, yielding ∆AIC = 1.97 for the CMB data sets and ∆AIC = 3.72 after including large scale structure measurements. The AIC takes into account how well the model fits the data while incorporating a penalty proportional to the number of parameters fit. When comparing two models, the lower AIC corresponds to the preferred model. In principle, the Bayes factor should be used to compare the models and the AIC proves only an approximation to it [see, e.g., the introduction of 45, for the Bayes factor].
In table 5 we list the individual ∆AIC values for each data set used in the CMB and CMB+LSS combinations, as well as the total ∆AIC value for both combinations. All ∆AIC values are positive, indicating the α B = 0 model is preferred for all data sets. For this case the data are consistent with GR.
Incorporating the LSS data leads to a lower preferred value and upper bound on M 0 (see table 2 and table 3). Increasing M 0 boosts the power of the matter power spectrum for large scales, increasing σ 8 (see figure 6). A weak correlation between M 0 and σ 8 is visible in the lower boundary of the relevant contour in figure 1. The RSD measurement's preference for a lower σ 8 helps shrink the contour and leads to a slightly lower preferred M 0 .
A primary goal of current and future CMB polarization experiments is to make the first detection of B-mode polarization from primordial gravitational waves. With our constraints we can investigate the impact of the Horndeski parameters on such experiments. In figure 7 we plot the primordial B-modes while varying α M 0 in the 95% CL range allowed by the CMB+LSS data set combination for α B 0 = 0 to show the range of B-mode spectra consistent with the data. We use r = 0.05, n t = 0, the best-fit Planck TT+lowP+lensing+ext ΛCDM parameter values [41], and our best-fit α B 0 and ξ, the power law index for braiding, values for the CMB+LSS data set combination. Even for the largest α M 0 allowed by the data, the deviations from GR are less than 1% at the reionization bump at < 10, which is smaller than the cosmic variance. Thus, even without a detection of primordial B-modes to date, we can conclude that the effect of the Horndeski parameters on the primordial B-mode polarization is constrained to be insignificant.

Summary
With a fixed kineticity of α K = 0.1a 3 we have shown: • The friction α M 0 has a 95% CL upper limit of 0.38 and 0.41 when α B 0 = 0 and α B 0 = 0, respectively, when using both CMB and LSS data.
• The lower 95% CL limit for α M 0 excludes GR for the α B 0 = 0 case but not for the α B 0 = 0 case. We believe this is primarily due to tachyon instabilities (defined in section 3) imposed by fixing α K 0 = 0.1.
• The effects of Horndeski theory on primordial B-modes are constrained by CMB and LSS data to be insignificant with 95% confidence.
It is important to remember that even when using Horndeski theory, making different assumptions about the α i parameters can yield dramatically different results. See e.g. the tight α M 0 constraint by [19], the negative values allowed by [12,24], and [12, erratum of], and the constraints found by [39]. We stress that choice of kineticity has a non-negligible -13 -

JCAP12(2018)030
impact on the constraints for the other Horndeski parameters due to its effects on the viable parameter space.
The observation of the electromagnetic counterpart to GW170817 has constrained −6 × 10 −15 ≤ α T 0 ≤ 1.4 × 10 −15 , which seems to eliminate all Horndeski theories with α T = 0. Typically, the α i 's are parameterized so that they are 0 in the matter dominated era and only have late time effects [see e.g. 20,21]. However, the evolution of α T could take a form such that α T → 0 as z → 0, but where α T = 0 in the past. In this case, eq. (1.1) is still viable. The power law evolution probed in this paper does not permit this behavior, but other functional forms can. It would be interesting to explore the parameter space and constraints from such a theory.

Acknowledgments
CK would like to thank the referee for their very helpful comments, as well as Marco Raveri, Paul Steinhardt, Anna Ijjas, Alexandre Barreira, Kris Pardo, and Alice Pisani for useful discussions. This work was partially funded through a DAAD (German Academic Exchange Service) Research Grant and the National Science Foundation Graduate Research Fellowship, NSF award number DGE1656466.
where a 2 δg 00 = a 2 g 00 + 1, are the perturbations to the time-time metric component, extrinsic curvature, curvature trace, and the three dimensional spatial Ricci scalar, respectively. n µ is the normal to surfaces of constant time. We adopt EFTCAMB's definitions for clarity. Ω, Λ, and c are the background evolution equations, as a function of conformal time [14,15]. The background -14 -

JCAP12(2018)030
equations Λ and c can be written as a function of Ω, the coupling to gravity: where ρ m and P m are the matter energy density and pressure, respectively [14,15]. EFT-CAMB multiplies R by 1 + Ω rather than Ω for numerical accuracy. S m [χ i , g µν ] is the action for all matter fields χ i . The Stückelberg technique makes the scalar perturbations explicit in unitary gauge. The conformal time is perturbed by a scalar field π, known as the Stückelberg field. All equations are now functions of τ + π, and the perturbation operators transform as [14,15] δg 00 → δg 00 − 2π The Horndeski theory of gravity is the most general tensor-scalar theory in which the equations of motion are second order [1]. However, the authors of [4] and [49] have presented an extended Horndeski theory in which the equations of motion have higher order derivatives. The equations of motion that describe the propagating degrees of freedom reduce to second order equations and, thus, avoid Ostrogradski instabilities [see also 50,51]. In this study we restrict ourselves to ordinary Horndeski theory in which operators contain at most two derivatives. The authors of [52] detail the derivatives introduced by the perturbation operators that act on the metric and scalar field perturbations. They note both (K µ µ ) 2 and δK µ ν δK ν µ contain terms with four spatial derivatives on the scalar perturbations and with one time and two spatial derivatives. Cancelling the two operators removes the four spatial derivatives, while δg 00 δR (3) can cancel with the mixed time and spatial derivative term. The ∂ µ a 2 g 00 ∂ ν a 2 g 00 term also contains higher order derivatives that cannot cancel with any other term, so it must be removed. The coefficient relationships required for these cancellations to occur are [52]  The authors of [2] have formulated a physically motivated parameterization of the coefficients in the action displayed in eq. (A.1) for Horndeski theory. The following four parameters are independent of both themselves and the background:

B Effects of kineticity on stability constraints
Kineticity's effects on the observables are closely related to the accuracy of the quasi-static approximation. The authors of [53] have shown that when analyzing the CMB, the quasi-static limit should be used neither when dark energy has a non-negligible effect at recombination nor when modeling the integrated Sachs-Wolfe (ISW) effect. If the dark energy sound speed, also known as the π field sound speed, is less than 0.1, the quasi-static limit is not valid for CMB lensing. When the approximation is valid, kineticity does not enter the equations of motion and is not well constrained by observations. As mentioned in section 2, ref. [12] found that α K could not be well constrained with their parameterization. They then presented constraints on α T , α M , and α B for a few fixed values of kineticity. To limit the number of additional degrees of freedom in our analysis, we also fix α K 0 in our analysis. With α T 0 = 0, we find that evolving the remaining parameters as constants, with α K = 0.1, yields a stable parameter space too small to explore with an MCMC analysis. Evolving the Horndeski parameters as power laws (see section 2 and appendix C) enlarges the stable parameter space and provides the opportunity to probe the time evolution of the parameters.
Although the authors of [12] found the remaining Horndeski parameters have a weak dependence on kineticity, we note that fixing the kineticity has non-trivial effects on our viable parameter space. We believe this is primarily due to the tachyon instabilities defined in section 3. In figure 8 we explore the stable parameter space for α K 0 = 0.1. We fix κ = 3.0, which we found provides a large range of likely values for the other Horndeski parameters. The stable region for β, the power law index related to friction, has a gap from 2 β 3 due to mathematical stability conditions. Without using a nested sampling method, the MCMC cannot reach both stable β regions. We choose to only explore the smaller friction exponent values β 2 since large values will drive the friction parameter α M close to 0. Having M 0 ≈ 0 effectively achieves the same result if the data prefer α M ≈ 0.
Stability requires the scalar propagation speed to satisfy [2] using the definition for α B of [33]. Because both kineticity and braiding are in the denominator, it is not stable to have α K 0 = 0 and α B 0 = 0 simultaneously. Thus, to explore the α B 0 = 0 case and compare with the α B 0 = 0 case in a self-consistent manner, we must fix α K 0 at a nonzero value. We analyzed the stable parameter space for α K 0 = 0 and α K 0 = 1.0, as well. Increasing α K 0 from 0.1 to 1.0 significantly shrinks the stable parameter space to the extent it cannot be explored with an MCMC analysis. We fix α K 0 = 0.1 throughout our analysis due to its larger stable parameter space and closeness to GR.
In figure 9 we compare the stable parameter space for α K 0 = 0 and α K 0 = 0.  kineticity values is dramatic, indicating the importance of understanding how all imposed stability conditions affect the viable parameter space to correctly interpret the parameter posteriors. When α K 0 = 0.1, the viable parameter space when 0 M 0 0.1 is small, whereas for α K 0 = 0 the same M 0 values provide a larger stable parameter space. The tachyon stability conditions solely drive these differences and those in the other stability contours. For the α K 0 = 0.1 case, the viable M 0 − α B 0 space converges to a single point as both α B 0 and M 0 approach 0 + , i.e. their GR value from the right. Thus, the choice of kineticity affects the M 0 posterior due to the tachyon stability constraints (see section 4). Indeed, if we remove the mathematical stability conditions the viable parameter spaces for α K 0 = 0 and α K 0 = 0.1 are identical.
We also note that while the stability contours for α K 0 = 0 make it appear α B 0 = 0 is stable for this kineticity, there are no stable points when both α K 0 = 0 and α B 0 = 0 exactly. The appearance of the plot is an artifact of discretely sampling the parameter space.

C Parameter evolution
To better constrain the Horndeski parameters we explore the pivot point a * at which to measure the parameters. The parameters of interest to measure then become parameters of interest as well as σ 8 as a function of pivot redshift z * . α B * is relatively uncorrelated with the other parameters at decoupling, which could be an artifact of the power law evolution. α M * and σ 8 are positively correlated at present. This is consistent with the effects seen in the matter power spectrum when increasing M 0 (see section 5). Large scales entered the horizon at low redshifts when α M * and σ 8 were positively correlated. Hence, increasing M 0 creates an amplifying effect on large scales. α B * and α M * are positively correlated with their respective exponents at the present epoch. However, the only epoch during which α B * and α M * are uncorrelated is at present. α B * has the smallest relative uncertainty at present. The relative uncertainty for α M * has a small minimum near z * ≈ 10 − 20, similar redshifts to where α M * is uncorrelated with β and σ 8 . This suggests a pivot point earlier than the present may be preferable for α M . However, the relative uncertainty for α M * at present is almost as small. Thus, because α B * and α M * are uncorrelated at present and the relative uncertainties for each are minimized or close to the minimum, we choose to keep the pivot point at present.