Vainshtein regime in Scalar-Tensor gravity: constraints on DHOST theories

TL;DR

This paper investigates the Vainshtein screening mechanism in general scalar-tensor theories compatible with gravitational wave observations, deriving constraints on model parameters from Solar System and pulsar tests.

Contribution

It derives the fully nonlinear effective Lagrangian for these theories, identifies key parameters, and constrains them using astrophysical observations.

Findings

For $eta_1 geq 0$, the potentials do not behave as in standard gravity.

Solar System light bending constrains $eta_1 o 10^{-5}$.

Pulsar timing constrains $eta_1 o 10^{-2}$.

Abstract

We study the screening mechanism in the most general scalar-tensor theories that leave gravitational waves unaffected and are thus compatible with recent LIGO/Virgo observations. Using the effective field theory of dark energy approach, we consider the general action for perturbations beyond linear order, focussing on the quasi-static limit. When restricting to the subclass of theories that satisfy the gravitational wave constraints, the fully nonlinear effective Lagrangian contains only three independent parameters. One of these, $\beta_1$, is uniquely present in degenerate higher-order theories. We compute the two gravitational potentials for a spherically symmetric matter source and we find that for $\beta_1 \ge 0$ they decrease as the inverse of the distance, as in standard gravity, while the case $\beta_1 < 0$ is ruled out. For $\beta_1 > 0$, the two potentials differ and their gravitational constants are not the same on the inside and outside of the body. Generically, the bound on anomalous light bending in the Solar System constrains $\beta_1 \lesssim 10^{-5}$. Standard gravity can be recovered outside the body by tuning the parameters of the model, in which case $\beta_1 \lesssim 10^{-2}$ from the Hulse-Taylor pulsar.

Equations

$$\displaystyle\begin{split}S_{\rm DHOST}=\ &\int d^{4}x\sqrt{-g}\Big{[}P(\phi,X)+Q(\phi,X)\Box\phi\\ &+f(\phi,X){}^{(4)}\!R+\sum_{I=1}^{5}a_{I}(\phi,X)L_{I}(\phi,\phi_{;\nu},\phi_{;\rho\sigma})\Big{]}\;,\end{split}$$

$$L_1$$

$$L_4$$

$$\begin{split}&S_{\rm EFT}=\int d^{4}x\sqrt{h}\frac{M^{2}}{2}\big{[}-(1+\delta N)\delta{\cal K}_{2}+c_{T}^{2}{}^{(3)}\!R+4H\alpha_{\text{B}}\delta K\delta N\\ &+(1+\alpha_{\text{H}}){}^{(3)}\!R\delta N+4\beta_{1}\delta KV+\beta_{2}V^{2}+{\beta_{3}}a_{i}a^{i}+\alpha_{\text{V}}\delta N\delta{\cal K}_{2}\big{]}\;,\end{split}$$

$$\begin{array}{rl}{\alpha }_{\text{B}}& ={\alpha }_{\text{V}}-3{\beta }_1+\dot{\varphi }(f_{\varphi }+2X{f}_{,\varphi X}+X{Q}_{,X})/(M^{2}H),\\ c_T^2& =2f/M^2,\\ {\alpha }_{\text{H}}& =4X(a_2-f_{,X})/M^2,\\ {\beta }_1& =2X(f_{,X}-a_2+a_{3}X)/M^2,\\ {\beta }_2& =-8X^2\left(a_3+a_4-2a_{5}X\right)/M^2,\\ {\beta }_3& =-8X(f_{,X}-a_2-a_{4}X)/M^2,\\ {\alpha }_{\text{V}}& =4X(f_{,X}-2a_2-2X{a}_{2,X})/M^2.\end{array}$$

$${\beta }_2=-6{\beta }_1^2,{\beta }_3=-2{\beta }_1\left[2(1+{\alpha }_{\text{H}})+{\beta }_1{c}_T^2\right],$$

$${\alpha }_{\text{M}}\equiv \frac{d\ln M^2}{d\ln a}.$$

$$d{s}^2=-(1+2\Phi )d{t}^2+a^2(t)(1-2\Psi )d{\vec{x}}^2.$$

$$S_{\mathrm{EFT}}=\int d^{4}x\left[\frac{M^{2}a}{2}\left({\mathcal{L}}_2+\frac{1}{a^2}{\mathcal{L}}_3+\frac{1}{a^4}{\mathcal{L}}_4\right)+{\mathcal{L}}_{\mathrm{m}}\right],$$

$$\begin{array}{rl}{\mathcal{L}}_2=& \left(c_1\Phi +c_2\Psi +c_3\pi \right){\partial }^2\pi +c_4\Psi {\partial }^2\Phi +c_5\Psi {\partial }^2\Psi \\ & +c_6\Phi {\partial }^2\Phi +(c_7\dot{\Psi }+c_8\dot{\Phi }+c_9\ddot{\pi }){\partial }^2\pi ,\\ {\mathcal{L}}_3=& -\frac{b_1}{2}(\partial \pi {)}^2{\partial }^2\pi +(b_2\Phi +b_3\Psi )Q_2[\pi ,\pi ]\\ & -\frac{1}{2}(\partial \pi {)}^2(b_4{\partial }^2\Psi +b_5{\partial }^2\Phi +b_6{\partial }^2\dot{\pi }),\\ {\mathcal{L}}_4=& -\frac{d_1}{2}(\partial \pi {)}^2{Q}_2[\pi ,\pi ]+\frac{d_2}{4}{\partial }_k(\partial \pi {)}^2{\partial }_k(\partial \pi {)}^2,\end{array}$$

$${\mathcal{L}}_{\mathrm{m}}=-a^3\Phi \delta {\rho }_{\mathrm{m}},$$

$${\mathcal{E}}_{{\phi }_a}\equiv \frac{1}{2a^2{M}^2}\frac{\delta S_{\mathrm{EFT}}}{\delta {\phi }_a}=0,$$

$$c_{\mathrm{T}}^2=1,\text{and}{\alpha }_{\mathrm{V}}=-{\alpha }_{\mathrm{H}},\text{(speed of gravity)}.$$

$${\alpha }_{\mathrm{H}}+2{\beta }_1=0,\text{(no decay)}.$$

$$\begin{split}S=\ &\int d^{4}x\sqrt{-g}\Big{[}P+Q\Box\phi+f{}^{(4)}\!R+\frac{3f_{,X}^{2}}{2f}\phi^{;\mu}\phi_{;\mu\nu}\phi^{;\nu\rho}\phi_{;\rho}\Big{]}\;.\end{split}$$

$$\begin{array}{rl}{\beta }_1& =X{f}_{,X}/f,{\alpha }_{\text{M}}=f_{,\varphi }\dot{\varphi }/(fH),\\ {\alpha }_{\text{B}}& =-X{f}_{,X}/f+\dot{\varphi }(f_{,\varphi }+2X{f}_{,\varphi X}+X{Q}_{,X})/(2fH).\end{array}$$

$$\begin{split}\alpha c_{s}^{2}\equiv\ &{2}(1+\alpha_{\text{B}}-\dot{\beta}_{1}/H)^{2}\bigg{[}\frac{1}{aM^{2}}\frac{d}{dt}\left(\frac{aM^{2}(1-\beta_{1})}{H(1+\alpha_{\text{B}})-\dot{\beta}_{1}}\right)-1\bigg{]}\\ &\qquad-\frac{\bar{\rho}_{\rm m}(1-\beta_{1})^{2}}{H^{2}M^{2}}\;,\end{split}$$

$${\partial }^2\pi =-\frac{a^2}{2M_{\mathrm{Pl}}{\nu }_2}\left({\nu }_4\delta {\rho }_{\mathrm{m}}+{\nu }_5\delta {\dot{\rho }}_{\mathrm{m}}\right),$$

$$\begin{array}{rl}{\nu }_2& \equiv \frac{M^2{H}^2\alpha c_s^2}{2M_{\mathrm{Pl}}(1-{\beta }_1)},\\ {\nu }_4& \equiv \frac{-H[{\alpha }_{\mathrm{B}}-{\alpha }_{\mathrm{M}}(1-{\beta }_1)+{\beta }_1(4-3{\beta }_1)]}{1-{\beta }_1},\\ {\nu }_5& \equiv -{\beta }_1.\end{array}$$

$$\begin{array}{rl}{\partial }^2\Phi & ={\mu }_{\Phi }\delta {\rho }_{\mathrm{m}}+{\nu }_{\Phi }{\delta \dot{\rho }}_{\mathrm{m}}+{\sigma }_{\Phi }{\delta \ddot{\rho }}_{\mathrm{m}},\\ {\partial }^2\Psi & ={\mu }_{\Psi }\delta {\rho }_{\mathrm{m}}+{\nu }_{\Psi }{\delta \dot{\rho }}_{\mathrm{m}}+{\sigma }_{\Psi }{\delta \ddot{\rho }}_{\mathrm{m}},\end{array}$$

$$x\equiv \frac{1}{{\Lambda }^3}\frac{{\pi }^{\prime }}{a^{2}r},y\equiv \frac{1}{{\Lambda }^3}\frac{{\Phi }^{\prime }}{a^{2}r},z\equiv \frac{1}{{\Lambda }^3}\frac{{\Psi }^{\prime }}{a^{2}r},$$

$$\mathcal{A}\equiv \frac{1}{8\pi M_{\mathrm{Pl}}{\Lambda }^3}\frac{m}{a^3{r}^3},m\equiv 4\pi {\int }_0^{r}d\tilde{r}{\tilde{r}}^{2}a(t{)}^3\delta {\rho }_{\mathrm{m}}(t,\tilde{r}),$$

$${\nu }_1{x}^2+\left[{\nu }_2+{\nu }_3{r}^{-2}(r^3\mathcal{A}{)}^{\prime }\right]x+{\nu }_4\mathcal{A}+{\nu }_5\dot{\mathcal{A}}=0,$$

$$\begin{array}{rl}{\nu }_1& =\xi M_{\mathrm{Pl}}^{-1}{M}^{2}H{\Lambda }^3,{\nu }_3={\Lambda }^3{\beta }_1,\end{array}$$

$$\xi \equiv 2{\alpha }_{\text{B}}-{\alpha }_{\text{M}}(1-{\beta }_1)+2{\beta }_1-2{\dot{\beta }}_1/H.$$

$$r_V\sim {\left(\frac{m}{8\pi M_{\mathrm{Pl}}}\right)}^{1/3}\frac{1}{a\Lambda }.$$

$$\begin{array}{rl}x& =-\frac{{\nu }_2+{\nu }_3\kappa \mathcal{A}-\sqrt{{\left({\nu }_2+{\nu }_3\kappa \mathcal{A}\right)}^2-4{\nu }_1({\nu }_4\mathcal{A}+{\nu }_5\dot{\mathcal{A}})}}{2{\nu }_1},\end{array}$$

$$\kappa (t,r)\equiv \frac{\partial \ln m(t,r)}{\partial \ln r}\ge 0.$$

$$\mathcal{A}|{\beta }_1|\kappa \gg \frac{{\beta }_1{\nu }_2}{{\nu }_3}\sim \frac{\alpha c_s^2}{2(1-{\beta }_1)},$$

$$\mathcal{A}{\left({\beta }_1\kappa \right)}^2\gg \frac{4{\beta }_1^2{\nu }_1({\nu }_4-3H{\nu }_5)}{{\nu }_3^2}\sim \frac{4\xi ({\nu }_4-3H{\nu }_5)}{H},$$

$$\mathcal{A}\sim 10^{30}\frac{m}{M_{\odot }}{\left(\frac{R_{\odot }}{r}\right)}^3,$$

Full text

I Introduction

Scalar-tensor theories are currently used to extend gravity beyond General Relativity on cosmological scales. The higher-derivate terms that characterize Horndeski Horndeski:1974wa ; Deffayet:2011gz and beyond-Horndeski theories Zumalacarregui:2013pma ; Gleyzes:2014dya ; Langlois:2015cwa ; Crisostomi:2016czh ; BenAchour:2016fzp (see also Langlois:2018dxi ; Kobayashi:2019hrl for recent reviews), while possibly providing an origin for the observed accelerated expansion of the Universe, are crucial to suppress modifications of gravity on Solar-System, where very stringent tests apply Will:2014xja .

Indeed, close to matter sources, the nonlinearities of the scalar field suppress its effects through the so-called Vainshtein screening mechanism Vainshtein:1972sx ; Babichev:2013usa , allowing General Relativity to be recovered. Screening is essential for the theories under consideration to be observationally viable.

However, recent observations of gravitational waves have put severe constraints on higher-derivative operators. The simultaneous observation of gravitational waves and gamma ray bursts from the GW170817 event has constrained very precisely the relative speed between gravitons and photons TheLIGOScientific:2017qsa , eliminating part of the derivative couplings of the scalar field to the curvature Creminelli:2017sry ; Sakstein:2017xjx ; Ezquiaga:2017ekz ; Baker:2017hug .111As discussed in deRham:2018red , this conclusion can be evaded if new physics appears at a scale parametrically smaller than the observed LIGO/Virgo frequencies. More recently, it has been pointed out that the theories beyond Horndeski belonging to the Gleyzes-Langlois-Piazza-Vernizzi (GLPV) class Gleyzes:2014dya ; Gleyzes:2014qga display a cubic graviton-scalar-scalar interaction that can mediate an observable decay of the graviton into dark energy particles Creminelli:2018xsv . The same vertex is responsible for an anomalous gravitational wave dispersion, when the speeds of scalar and gravitational wave differ. These two effects may become important at frequencies relevant for LIGO/Virgo observations, leading to a very tight constraint on these theories.222Recently, positivity bounds for scalar theories coupled to gravity Bellazzini:2019xts have decreased significantly the regime of validity of the corresponding EFTs (e.g. for the cubic galileon Nicolis:2008in ). It would be interesting to study how these theoretical constraints affect general theories with screening.

In this article we study the subclass of scalar-tensor theories that leave gravitational waves unaffected and satisfy these constraints exactly. This study is the natural continuation of Ref. Kobayashi:2014ida , which examined the Vainshtein mechanism in GLPV theories, and of Refs. Crisostomi:2017lbg ; Langlois:2017dyl ; Dima:2017pwp , where this investigation was extended to Degenerate Higher-Order Scalar-Tensor (DHOST) theories Langlois:2015cwa ; Crisostomi:2016czh ; BenAchour:2016fzp , in the case where the speed of gravitational waves equals the one of light. These references showed that for the subclass of these theories extending Horndeski, the two gravitational potentials differ inside matter and depend on the density profile, signalling the breaking of the Vainshtein screening. (See for instance also Bartolo:2017ibw ; Ganz:2018vzg ; Babichev:2018rfj ; Kase:2018iwp ; Crisostomi:2017pjs ; Crisostomi:2018bsp ; Frusciante:2018tvu for other studies on the surviving theories and Kase:2018aps ; Ezquiaga:2018btd for reviews.)

As in these references, here we consider the solutions for the gravitational potentials in these theories near matter sources, i.e. in the Vainshtein regime, and compare them to observational data. We will rely on the quasi-static approximation, which is valid on scales much smaller than the Hubble radius when we restrict to non-relativistic sources.

In the next section we review DHOST theories and their effective description, and we expand the action in the metric potentials and scalar field perturbations. We then discuss the subset of theories that evade the gravitational wave constraints. We briefly study the linear theory in Sec. III and the Vainshtein regime around a spherically symmetric source in Sec. IV. Here, we also discuss our results in the more familiar Horndeski frame. The previous discussions do not apply to theories that are related to General Relativity by a conformal transformation, which are the subject of Sec. V. Finally, constraints on the parameters are derived in Sec. VI, and conclusions are left to Sec. VII.

To simplify the main text, many of the coefficients and equations have been postponed to App. A. In App. B we also present some extra astrophysical constraints, and we show that they are weaker than the ones obtained from the Solar System and the Hulse-Taylor pulsar.

II Action and perturbation equations

II.1 DHOST theories

We denote by the semicolon a covariant derivative and we define $X\equiv-\phi_{;\mu}\phi^{;\mu}/2$. The action for DHOST theories includes all possible quadratic combinations up to second derivatives of the field $\phi$ and reads Langlois:2015cwa

[TABLE]

where ${}^{(4)}\!R$ is the 4D Ricci scalar and the $L_{I}$ are defined by

[TABLE]

The functions $P$ and $Q$ do not affect the degeneracy character of the theory. Instead, for DHOST theories, the functions appearing in the second line of eq. (1) must satisfy three degeneracy conditions Langlois:2015cwa that fix three of these functions in terms of the others. Here we are going to focus on the subclass that satisfy $a_{1}+a_{2}=0$ (and two other degeneracy conditions). Other subclasses have been shown to display a linear instability, either in the scalar or in the tensor sector deRham:2016wji ; Langlois:2017mxy .

II.2 Effective description and constraints

It is convenient to discuss observational constraints on these theories in terms of the EFT of dark energy parameters, which for DHOST theories have been introduced in Langlois:2017mxy and extended to nonlinear order in the perturbations in Dima:2017pwp .

Specifically, in the presence of a preferred slicing induced by a time-dependent scalar field, we can choose the time as to coincide with the uniform field hypersurfaces. In this gauge, called the unitary gauge, and using the ADM metric decomposition with line element $ds^{2}=-N^{2}dt^{2}+h_{ij}(dx^{i}+N^{i}dt)(dx^{j}+N^{j}dt)$, cosmological perturbations around an FRW solution $ds^{2}=-dt^{2}+a^{2}(t)d\vec{x}^{2}$ are governed by the action

[TABLE]

where we have written only the operators with the highest number of spatial derivatives, which are relevant in the quasi-static limit. Here $H\equiv\dot{a}/a$ (a dot denotes the time derivative), $\delta N\equiv N-1$, $\delta K_{i}^{j}\equiv K_{i}^{j}-H\delta_{i}^{j}$ is the perturbation of the extrinsic curvature of the time hypersurfaces, $\delta K$ its trace, and ${}^{(3)}\!R$ is the 3D Ricci scalar of these hypersurfaces. Moreover, $\delta{\cal K}_{2}\equiv\delta K^{2}-\delta K_{i}^{j}\delta K^{i}_{j}$, $V\equiv(\dot{N}-N^{i}\partial_{i}N)/N$ and $a_{i}\equiv\partial_{i}N/N$.

The time-dependent functions in this action are related to the free functions in eq. (1). One finds that the effective Planck mass that normalizes the graviton kinetic energy is given by $M^{2}=2(f-2a_{2}X)$. The other parameters read Dima:2017pwp

[TABLE]

The function $\alpha_{\text{B}}$ Bellini:2014fua measures the kinetic mixing between metric and scalar fluctuations Creminelli:2008wc and $c_{T}^{2}$ the fractional difference between the speed of gravitons and photons. The function $\alpha_{\text{H}}$ measures the kinetic mixing between matter and the scalar fluctuations Gleyzes:2014dya ; Gleyzes:2014qga ; DAmico:2016ntq and vanishes for Horndeski theories, while the function $\alpha_{\text{V}}$ parameterizes the only operator that starts cubic in the perturbations Cusin:2017mzw .

Finally, the functions $\beta_{1}$, $\beta_{2}$, and $\beta_{3}$ parameterize the presence of higher-order operators. The degeneracy conditions, which we assume hereafter, read Langlois:2017mxy

[TABLE]

so that only $\beta_{1}$ is independent. Another function that we define here for later convenience is Bellini:2014fua

[TABLE]

II.3 Action in Newtonian gauge

To study scalar linear and higher-order perturbations we will exit the unitary gauge and work in the Newtonian gauge, where the metric is written as

[TABLE]

The scalar field $\pi$ is introduced by performing a space-time dependent shift in the time $t\rightarrow t+\pi(t,\vec{x})$. Then, we expand the action eq. (LABEL:EFTaction) in terms of the metric and scalar field perturbations. We keep only terms with the highest number of derivatives per field, which are those relevant in the quasi-static limit, and we find

[TABLE]

with

[TABLE]

where we have defined $Q_{2}[\varphi_{a},\varphi_{b}]\equiv\varepsilon^{ikm}\varepsilon^{jlm}\partial_{i}\partial_{j}\varphi_{a}\partial_{k}\partial_{l}\varphi_{b}$ for $\varphi_{a}\equiv\{\Phi,\Psi,\pi\}$, and $c_{1},\ldots,c_{9}$, $b_{1},\ldots,b_{6}$ and $d_{1},d_{2}$ are time-dependent coefficients, reported in App. A.1.

The last term in the bracket is the matter Lagrangian. If we define by $\rho_{\rm m}$ the matter energy density and by $\bar{\rho}_{\rm m}$ its mean cosmological value, this is given by

[TABLE]

where $\delta\rho_{\rm m}\equiv\rho_{\rm m}-\bar{\rho}_{\rm m}$.

The field equations can be derived straightforwardly by varying the action eq. (8) with respect to $\Phi$, $\Psi$ and $\pi$. We have

[TABLE]

where the explicit expressions of ${\cal E}_{\varphi_{a}}$ can be found in App. A.1 (see eqs. (70)–(72)).

The above Lagrangians and equations are valid for general DHOST theories. We will now consider the subclass of theories leaving the gravitational wave unaffected.

II.4 Gravitational wave constraints

We now focus on the subset of theories that evade the gravitational wave constraints. In particualr, we demand that gravity and light travel at the same speed, i.e. (see e.g. Creminelli:2017sry ),

[TABLE]

Moreover, we require that gravitons do not decay into dark energy by setting Creminelli:2018xsv

[TABLE]

Unless otherwise stated, in the following we impose these two equations and replace $\alpha_{\text{H}}$ and $\alpha_{\text{V}}$ in terms of $\beta_{1}$.

In terms of the functions $a_{I}$ in the action eq. (1), these requirements, together with the degeneracy conditions, read $a_{1}=a_{2}=a_{3}=a_{5}=0$ and $a_{4}=3f_{,X}^{2}/({2f})$, for any $X$, i.e., Creminelli:2018xsv ,

[TABLE]

From eq. (4), in this subset of theories the remaing free parameters are thus given by

[TABLE]

III Linear regime

We briefly discuss how matter inhomogeneities source $\pi$ and the gravitational potentials at linear order. For convenience, we first define the combination

[TABLE]

where $c_{s}^{2}$ is the effective sound speed of dark energy fluctuations Langlois:2017mxy , which must be positive to avoid gradient instabilities, while the time-dependent function $\alpha$ is the coefficient in front of the time kinetic term of scalar fluctuations (see Langlois:2017mxy ), which must be positive to avoid ghosts. In practice, we do not need their explicit expressions because in the quasi-static limit these two parameters always appear in the combination $\alpha c_{s}^{2}>0$.

By solving for $\Phi$ and $\Psi$ the linear equations obtained by varying the action with respect to the two metric potentials, and replacing these solutions in the scalar field equation, we obtain

[TABLE]

where $M_{\rm Pl}$ is the value of the Planck mass measured today and the parameters $\nu_{2}$, $\nu_{4}$ and $\nu_{5}$ above are defined as

[TABLE]

This choice of definitions will become clearer when we consider the full nonlinear equation for $\pi$, in Sec. IV.

For completeness, we also provide the solutions for the metric potentials, which can be obtained by replacing eq. (17) back into the equations for $\Phi$ and $\Psi$. One finds Crisostomi:2017pjs ; Hirano:2019nkz

[TABLE]

where the coefficients on the right-hand side are given in App. A.2. Observational constraints on the above coefficients of the linear equations are discussed in Hirano:2019nkz . We are now ready to discuss the Vainshtein regime.

IV Vainshtein regime

In this section we want to study the Vainshtein mechanism around a spherically symmetric body, such as for instance a non-relativistic star (see Kobayashi:2018xvr for a study of relativistic stars in DHOST theories). This is at play close to the body, where non-linearities of the scalar field become important suppressing the scalar force. Far from the body, the linear solutions affected by the fifth force (see eq. (19)) are recovered. This can allow the theory to be compatible with stringent Solar System observations while at the same time modifying gravity on large scales.

To derive the equations relevant for spherically symmetric solutions, we assume that all of the fields depend only on time and the radial variable, $\varphi_{a}(t,\vec{x})=\varphi_{a}(t,r)$, where $r\equiv|\vec{x}|$. Then, we integrate the field equations eq. (11) over the radial variable and use Stoke’s theorem. Following Kobayashi:2014ida , we use the following notation,

[TABLE]

where a prime denotes a derivative with respect to $r$ and $\Lambda$ is a mass scale of order $\Lambda^{3}\sim H_{0}^{2}M_{\rm Pl}$, where $H_{0}$ is the Hubble rate today. Moreover, we define

[TABLE]

where $m$ is the physical mass of an overdensity contained in a spherical ball of physical radius $a(t)r$. The quantity ${\cal A}(t,r)$ represents the comoving mass density contrast in this ball.

The explicit expressions of the field equations in spherical symmetry are reported in App. A.3. Following Crisostomi:2017lbg ; Langlois:2017dyl ; Dima:2017pwp , we can solve the first two equations for $y$ and $z$ and plug the solutions into the third equation. After imposing the degeneracy conditions eq. (5), one obtains a polynomial equation for $x$ only. This equation was previously studied after imposing that gravitons and photons propagate with the same speed, eq. (12), in Crisostomi:2017lbg ; Langlois:2017dyl and in the general case in Dima:2017pwp and a cubic equation was obtained. Here we further restrict to theories where the gravitational waves do not decay, i.e. eq. (13). In this case one finds instead a quadratic equation,

[TABLE]

where the time-dependent functions $\nu_{2}$, $\nu_{4}$ and $\nu_{5}$ already appeared in the linear equations and are defined in eq. (18). The other functions are given by

[TABLE]

with

[TABLE]

If $\nu_{1}\neq 0$, the term proportional to $x^{2}$ becomes important close to the matter source, while it is negligible far away from it. The transition happens at the so-called Vainshtein radius, roughly corresponding to $\mathcal{A}\sim 1$, i.e.,

[TABLE]

For example, for a star like the sun, the Vainshtein radius is about a tenth of the size of the Milky Way.

For $\nu_{1}=0$, the term quadratic in $x$ vanishes and the usual Vainshtein screening cannot be supported. As shown in Sec. IV.3, this case corresponds to scalar-tensor theories that are conformally related to General Relativity, i.e., that can be described by the Einstein-Hilbert term plus conformally coupled matter. Here we will assume that $\nu_{1}$ (and thus $\xi$) does not vanish. We discuss the case $\xi=\nu_{1}=0$ in Sec. V.

As discussed in Sec. III, $\alpha c_{s}^{2}>0$. Thus, for $\beta_{1}<1$, $\nu_{2}$ defined in eq. (18) is also positive and the solution to eq. (22) that matches the linear regime, i.e., that has the correct behavior for $\mathcal{A}\ll 1$, is

[TABLE]

where for convenience we have defined the positive function,

[TABLE]

The terms proportional to $\nu_{3}$ and $\nu_{5}$ in eq. (22) contain respectively the radial and time derivative of the comoving mass of the body. Thus, when the mass of the central overdensity is constant in time and space, such as at a radial distance larger than the body size, we have $(r^{3}\mathcal{A})^{\prime}=0$ and $\dot{\mathcal{A}}=-3H\mathcal{A}$. Therefore, we consider two different cases, depending on whether we are outside or inside the object.

IV.1 Inside of matter

Inside of matter, $(r^{3}\mathcal{A})^{\prime}$ is generally different from zero, in which case, when $\mathcal{A}\gg 1$, the term proportional to $\nu_{3}$ can dominate the square root in eq. (26). This regime is characterized by

[TABLE]

and

[TABLE]

where on the right-hand side we have used the definitions in eq. (18) and eq. (23) and that $M^{2}H^{2}\sim M_{\rm Pl}\Lambda^{3}$. The right-hand sides of these equalities are of order unity while, on the left-hand side, ${\cal A}\gg 1$. For instance, inside a star like the sun,333Here and in the rest of this paper, we use the symbol $\odot$ to denote the appropriate quantity for our sun.

[TABLE]

where we have used $\Lambda\sim(10^{3}\,\text{km})^{-1}$. These conditions are thus realized inside matter unless $\beta_{1}\kappa$ becomes $\mathcal{O}(10^{-15})$. Using a mass profile from the Lane-Emden equation (see eq. (83) in the appendix), we have verified that this happens only on a very thin region near the surface of the star.

Now, expanding the solution in eq. (26) for $\kappa\neq 0$ and $\mathcal{A}\gg 1$, we can distinguish between two cases. When $\beta_{1}>0$, the solution for $x$ is given by

[TABLE]

This solution can then be plugged back into the field equations for $y$ and $z$ (see eqs. (77) and (78)) to solve for these two variables. Since $x_{\rm in}\sim{\cal O}({\cal A}^{0})$, the equations for $y$ and $z$ are dominated by the usual matter term linear in ${\cal A}$ and terms both linear and quadratic in $x$ can be neglected.

The solutions for the potentials can be straightforwardly computed and read

[TABLE]

where $G_{*}\equiv 1/(8\pi M^{2})$ is the gravitational constant that canonically normalizes the graviton, and

[TABLE]

Thus, in this theory $\Phi\neq\Psi$ so that Vainshtein screening is broken inside the body. Notice that the breaking is different from the one found in Kobayashi:2014ida ; Crisostomi:2017lbg ; Langlois:2017dyl ; Dima:2017pwp , which depends on the radial derivatives of the object mass.

The case $\beta_{1}<0$ is instead ruled out. Indeed, in this case the solution eq. (26) reads

[TABLE]

Therefore, $\pi$ is not suppressed inside matter and nonlinear terms proportional to $x^{2}\sim{\cal O}({\cal A}^{2})$ dominate the field equations for $y$ and $z$. The solutions,

[TABLE]

are incompatible with the existence of stars or other bounded objects.

IV.2 Outside of matter

Using eq. (21) outside the object, where the physical mass is constant, we have $\dot{\mathcal{A}}=-3H\mathcal{A}$ and $(r^{3}\mathcal{A})^{\prime}=0$. Replacing this in eq. (26), near the source we obtain

[TABLE]

In contrast with the solution inside matter, eq. (31), where $x_{\rm in}\sim{\cal O}({\cal A}^{0})$, here $x_{\rm out}^{2}\sim{\cal O}({\cal A})$ so that terms quadratic in $x$ contribute to the $y$ and $z$ solutions while linear terms are negligible.

The solutions for the potentials become now

[TABLE]

with

[TABLE]

where for convenience we have defined a new quantity,

[TABLE]

From these expressions, it is clear that the Vainshtein mechanism is broken also outside of the matter source because the two gravitational potentials are different. Only for

[TABLE]

are the two potentials the same and the Vainshtein screening recovered outside of matter.

An example of the solutions eq. (32) (for $r<R_{\odot}$) and eq. (37) (for $r>R_{\odot}$) is presented in Fig. 1. Notice the suppressed value of $x$. For a large overdensity ${\cal A}$, such as in a star, the transition between the interior of the star, where eqs. (28) and (29) are satisfied, and the exterior, where $\kappa=0$, takes place abruptly, so that $x$, $y$ and $z$ visually display a discontinuity for $\beta_{1}\upsilon\neq 0$. Of course, the solutions are actually continuous, as follows for instance from eq. (26).

IV.3 The Horndeski frame

The above results can be understood by noticing that DHOST theories that we consider in this paper can be mapped to Horndeski theories via an invertible $X$-dependent conformal and disformal transformation Crisostomi:2016czh ; Achour:2016rkg .

Since we want to preserve the speed of propagation of gravitons, we focus on conformal transformations only, i.e.

[TABLE]

and use the tilde to denote quantities in the Horndeski frame. For convenience, we also introduce the dimensionless time-dependent parameter

[TABLE]

with the right-hand side of this equation evaluated on the background solution.

The relations to compute the transformation of the EFT parameters of the action eq. (LABEL:EFTaction) under the change of frame eq. (41) are given in Langlois:2017mxy (see for instance eq. (2.22) of that reference). Without loss of generality, we will assume that $C=1$ on the background solution, so that $\tilde{M}^{2}=M^{2}$ and $\tilde{a}=a$. Moreover, for the transformation above, one finds

[TABLE]

while $\tilde{c}_{T}=c_{T}$.

Using these relations, one can show that a DHOST theory with gravitational metric $g_{\mu\nu}$ satisfying the conditions eq. (12) and eq. (13) (see its covariant form in eq. (14)), is equivalent to a Horndeski theory with gravitational metric $\tilde{g}_{\mu\nu}$,444Since matter is coupled to $g_{\mu\nu}$, the coupling to matter in the new frame $\tilde{g}_{\mu\nu}$ is non-minimal; see below. with $\tilde{c}_{T}=1$ and $\tilde{\alpha}_{\rm V}=\tilde{\alpha}_{\rm H}=\tilde{\beta}_{1}=0$ and action

[TABLE]

provided that $\alpha_{\rm Y}=-\beta_{1}$. Additionally, we obtain Langlois:2017mxy (see eq. (C.14) of that reference)

[TABLE]

Using the above transformations we can check that $\xi$ given in eq. (24) can be written in terms of Horndeski-frame quantities as $\xi=(1-\beta_{1})(2\tilde{\alpha}_{\rm B}-\tilde{\alpha}_{\rm M})$. As mentioned below eq. (24), this vanishes for $\tilde{\alpha}_{\rm M}=2\tilde{\alpha}_{\rm B}$, i.e. for Brans-Dicke Brans:1961sx , $f(R)$ Carroll:2004de or other theories conformally related to General Relativity Gleyzes:2015pma .

In the Horndeski frame, the EFT action eq. (8) is given by

[TABLE]

Here the Lagrangians $\tilde{\cal L}_{2}$, $\tilde{\cal L}_{3}$ and $\tilde{\cal L}_{4}$ have analogous expressions as those in eq. (9), but now fields and parameters are tilded and $\tilde{\beta}_{1}=\tilde{\beta}_{2}=\tilde{\beta}_{3}=0$.

Using eq. (41) and focusing on the leading order in spatial derivatives to retain only terms relevant in the quasi-static limit, the potentials in the Horndeski frame are related to those in the DHOST frame by

[TABLE]

while it is straigforward to verify that the scalar field fluctuation does not change,

[TABLE]

While the metric transforms when changing frame, matter remains always minimally coupled to the gravitational metric in the DHOST frame, i.e. $g_{\mu\nu}$. Therefore, $\tilde{\cal L}_{\rm m}={\cal L}_{\rm m}$. Using eq. (10) and the relations above, the coupling in the Horndeski frame is thus

[TABLE]

We can now vary the action eq. (46) with respect to $\tilde{\Phi}$, $\tilde{\Psi}$ and $\tilde{\pi}$. Actually, using the above functional relationships between the Horndeski action $\tilde{S}_{\rm EFT}$ and the DHOST action $S_{\rm EFT}$, we can easily derive the relationship between the Horndeski frame equations of motion, $\tilde{\mathcal{E}}_{\rm\tilde{\varphi}_{a}}\equiv(2a^{2}M^{2})^{-1}\delta\tilde{S}_{\rm EFT}/\delta\tilde{\varphi}_{a}=0$ and the DHOST frame equations of motion eq. (11) using the chain rule,

[TABLE]

This gives

[TABLE]

Using the field equations $\tilde{\cal E}_{\tilde{\varphi}_{a}}=0$, we can solve the system assuming spherical symmetry around a body, similarly to what was done earlier in Sec. IV. The equations for $\tilde{y}$ and $\tilde{z}$ read

[TABLE]

(i.e. all nonlinear terms vanish), and the closed equation for $\tilde{x}$ remains the same as that of $x$, eq. (22), as expected. Since terms linear in $\tilde{x}$ can be neglected both inside and outside of matter, the solutions to the above equations become, in terms of the potentials in the Horndeski frame,

[TABLE]

both inside and outside of the source, as long as ${\cal A}\gg 1$.

With eq. (53), we can now verify what we found in the previous subsections. Inside matter, $\pi$ can be neglected in eq. (47) and eq. (32) is recovered. Outside matter, $\dot{\pi}$ can be neglected but $(\partial\pi)^{2}$ cannot. Replacing in eq. (47) the solution for $x$, eq. (36), we can recover the solution outside the body, eq. (37).

Moreover, the conformal transformation eq. (47) leaves the sum of the potentials invariant. One can verify directly that our expressions inside eq. (32) and outside eq. (37) of the source satisfy $y+z=\tilde{y}+\tilde{z}=2G_{*}m/(r^{3}\Lambda^{3}(1-\beta_{1}))$, where we have used eq. (53) for the last equality. Additionally, the fact that the solutions in the Horndeski frame eq. (53) are valid both inside or outside of the source means that $\epsilon_{\Phi}^{\rm out}+\epsilon_{\Psi}^{\rm out}=\epsilon_{\Phi}^{\rm in}+\epsilon_{\Psi}^{\rm in}$, which can be verified directly in eq. (38).

V Theories conformally related to General Relativity

In the previous section we focused on theories with $\xi\neq 0$. Here we consider the case $\xi=0$, which corresponds to theories related to General Relativity by the conformal transformation (41), as shown in Sec. IV.3. Their general action is given by eq. (14) with $Q=0$. Examples are Jordan-Brans-Dicke Brans:1961sx and $f(R)$ theories Carroll:2004de , but here we are interested in extentions of these theories where the conformal factor relating them to General Relativity depends also on $X$, in which case $\beta_{1}\neq 0$.

When $\xi=0=\nu_{1}$, the equation for $x$, eq. (22), becomes linear. Inside of matter, the term proportional to $\nu_{3}$ dominates the solution so that one obtains the same conclusions as in Sec. IV.1 (including the constraint $\beta_{1}\geq 0$).

Thus, we focus on the solutions outside of matter and near the source. In that case, the solution for $x$ is simply the linear solution,

[TABLE]

Solving the equations for $y$ and $z$ with this solution we have, in terms of $\Phi^{\prime}$ and $\Psi^{\prime}$,

[TABLE]

This is clearly incompatible with Solar System tests unless $\beta_{1}\upsilon^{2}$ is practically zero.

A possible way out with $\beta_{1}\neq 0$ is to consider the subclass of theories with $\upsilon=0$. Taking into account that $\xi=0$, this fixes all parameters as a function of $\beta_{1}$, i.e.,

[TABLE]

For this particular case the source term in eq. (22) is absent and the fifth force vanishes, $x=0$. Thus, the gravitational potentials outside the source are linear in the object mass and coincide with the solutions eq. (32), which were obtained ignoring the scalar field contribution.

VI Constraints

Let us discuss observational constraints on the parameters of the theories studied above, in particular on $\beta_{1}$. A first class of constraints comes from stellar physics. We have studied them in App. B and have shown that they lead to complementary (but weaker) constraints to those derived in the following.

As shown in Jimenez:2015bwa , for $c_{T}=1$ the decrease of the orbital period of binary stars is proportional to the ratio between the gravitational constant normalizing the gravitons, $G_{*}$, and the one entering the Kepler law, which here is taken to be the one outside the object, $G_{*}(1+\epsilon^{\rm out}_{\Phi})$. Thus, using the results from Weisberg:2010zz , the Hulse-Taylor pulsar (PSR B1913+16) Hulse:1974eb allows us to constrain $\epsilon_{\Phi}^{\rm out}$,

[TABLE]

Next, we move to the Cassini constraints. Measurements of the frequency shift of radio waves, sent to and from the Cassini spacecraft as they passed near the sun Bertotti:2003rm , constrain the post-Newtonian parameter $\gamma_{\rm PN}\equiv\Psi/\Phi$ to be $-0.2\times 10^{-5}<\gamma_{\rm PN}-1<5.5\times 10^{-5}$. Since eq. (57) says that $\epsilon_{\Phi}^{\rm out}$ is small, we can approximate $\gamma_{\rm PN}-1\approx\epsilon_{\Psi}^{\rm out}-\epsilon_{\Phi}^{\rm out}$ and use this measurement to constrain the relative difference between the gravitational potentials outside of matter,

[TABLE]

Using eq. (37) we can rewrite the above quantity as

[TABLE]

where we remind the reader that $\xi$ and $\upsilon$ are respectively defined in eqs. (24) and (39). For generic values of $\alpha_{\text{B}}$, $\alpha_{\text{M}}$ and $\dot{\beta}_{1}$, one expects that $2-\upsilon/\xi\sim\mathcal{O}(1)$, so that the above turns into a tight constraint on $\beta_{1}$,

[TABLE]

The bound on the left-hand side comes from the result of Sec. IV.1 that negative values of $\beta_{1}$ are not allowed.

A more conservative constraint on $\beta_{1}$, independent of $\alpha_{\text{B}}$, $\alpha_{\text{M}}$, and $\dot{\beta}_{1}$, comes from isolating $\beta_{1}$ by using eq. (38). One obtains

[TABLE]

Combined with the constraint $\beta_{1}>0$, this gives

[TABLE]

which also holds in the presence of the tuning $2\xi-\upsilon\approx 0$.

The previous discussion assumed $\xi\neq 0$. Theories that are conformally related to General Relativity have $\xi=0$. As discussed in Sec. V, these theories are ruled out unless they satisfy the condition (56). In this case, the post-Newtonian parameter $\gamma_{\rm PN}$ reads

[TABLE]

Combining $\beta_{1}\geq 0$ with the Cassini constraint, eq. (58), gives

[TABLE]

VII Conclusion

Vainshtein screening is crucial to pass the stringent tests of gravity available at local scales. We studied the Vainshtein regime in the most general subset of DHOST theories that evade all of the constraints from gravitational wave observations, i.e., for which gravitons propagate at the speed of light and do not decay.

For non-zero values of the parameter that characterizes the higher-order derivatives in DHOST theories, i.e. $\beta_{1}$, the screening mechanism is broken. Negative values of $\beta_{1}$ are ruled out because the gravitational potentials inside matter do not scale as the inverse of the distance. For positive values of $\beta_{1}$, the scalar field can be neglected inside matter while its non-lineartities contribute to the potentials outside of matter. Therefore, we find that the gravitational potentials scale as the inverse of the distance but have different gravitational constants inside and outside of the body, and between themselves.

For generic values of the other parameters, in particular of $\alpha_{\text{B}}$, $\alpha_{\text{M}}$ and $\dot{\beta}_{1}$, the measurement of the Shapiro time-delay with the Cassini spacecraft constraints the parameter $\beta_{1}$ to be smaller than $10^{-5}$. But there is a special region of the parameter space where the two gravitational potentials are equal outside of the body and the bound from Cassini is evaded. Here the constraint on $\beta_{1}$ weakens, $\beta_{1}<7.5\times 10^{-3}$ at $2\sigma$, and is obtained from the measurement of the orbital decay rate of the Hulse-Taylor pulsar.

Theories conformally related to General Relativity have been considered separately. Generically, they are ruled out unless $\beta_{1}=0$. For specific values of the parameters $\alpha_{\text{B}}$ and $\alpha_{\text{M}}$ they admit $0\leq\beta_{1}\lesssim 10^{-6}$.

In App. B we studied also other constraints, coming from the stellar structure, that can be imposed using the difference between the Newtonian potential inside and outside of a star. However, they all turned out to be weaker than the bounds obtained using the exterior solutions. These constraints can be improved with more sophisticated methods (see e.g. Saltas:2018mxc ) and with more accurate data that will be available in the future.

Note added.— Another article Hirano:2019scf , whose content overlaps with ours, appeared while finalizing this article.

Acknowledgements

We thank the authors of Hirano:2019scf for interesting correspondence. MC thanks Francesco Villante, Santiago Casas, and Virginia Ajani for useful discussions. MC is supported by the Labex P2IO and the Enhanced Eurotalents Fellowship. ML acknowledges financial support from the Enhanced Eurotalents fellowship, a Marie Sklodowska-Curie Actions Programme, and the European Research Council under ERC-STG-639729, preQFT: Strategic Predictions for Quantum Field Theories.

Appendix A Coefficients and equations

We report here the explicit expressions of the coefficients introduced in the text and of some of the long expressions for the field equations.

A.1 Action and field equations

The coefficients appearing in the action expanded in perturbations, eq. (9), are given explicitly by

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]

Let us define

[TABLE]

The variation of the action eq. (8) gives the field equations eq. (11), with

[TABLE]

[TABLE]

and

[TABLE]

where we have defined $Q_{3}[\varphi_{a},\varphi_{b},\varphi_{c}]\equiv\varepsilon^{ikm}\varepsilon^{jln}\partial_{i}\partial_{j}\varphi_{a}\partial_{k}\partial_{l}\varphi_{b}\partial_{m}\partial_{n}\varphi_{c}$. Note that these equations are general: no degeneracy conditions or observational constraints have been assumed.

A.2 Linear solutions

Here we give the coefficients relevant for the linear solutions in eq. (19). These can be written quite compactly in terms of the coefficients in the linear solution for $\pi$ eq. (17), which we define as

[TABLE]

Then, in terms of these, the coefficients in eq. (19) are,

[TABLE]

and

[TABLE]

where

[TABLE]

A.3 Field equations in spherical symmetry

In spherically symmetry, the field equations eq. (70)–eq. (72) reduce to

[TABLE]

[TABLE]

and

[TABLE]

Appendix B Stellar constraints

In this appendix, we use properties of stars to constrain the relative difference between the gravitational constant inside and outside of matter,

[TABLE]

We assume that masses of stars in the literature are measured using the value of the gravitational constant outside of a source, i.e. $G_{\Phi}^{\rm out}=G_{*}(1+\epsilon_{\Phi}^{\rm out})$ in eq. (37), and that this is the gravitational constant value measured on the Earth, $G_{\Phi}^{\rm out}=6.6276\times 10^{-8}\text{ dyne cm}^{2}/\text{gm}^{2}$.

Here, we show that the constraints involving the interior of the star are weaker than the ones presented in Sec. VI. The first constraint is an upper bound and comes from the Chandrasekhar mass limit. The second comes from the minimum main-sequence mass of brown dwarfs and gives a lower bound. We also present a third constraint, coming from fitting the mass-radius relation for white dwarfs. This gives the best bound, i.e.,

[TABLE]

Of course, with improved observational data or theoretical modeling, these bounds could be improved in the future.

B.1 Chandrasekhar mass limit

First, we consider the Chandrasekhar mass limit, which is the largest mass that a white dwarf can have Chandrasekhar:1931ftj ; Chandrasekhar:1935zz . This limit exists because if the white dwarf had a larger mass, the electron degeneracy pressure would not be able to support the star, and it would collapse into a neutron star or a black hole. (Reference Jain:2015edg considered a similar constraint in the context of the DHOST theories which have a modification of the Newtonian potential eq. (32) proportional to $m^{\prime\prime}$.)

The Chandrasekhar mass limit $M_{\rm Ch}$ for white dwarfs is proportional to $(G^{\rm in}_{\Phi})^{-3/2}$ Chandrasekhar:1935zz , and the current theoretical value of this limit is $M_{\rm Ch}\approx 1.44M_{\odot}$. Taking into account the difference between the Newton constants inside and outside of the star, this limit becomes $M_{\rm Ch}\approx 1.44M_{\odot}(1+\epsilon_{G})^{-3/2}$. The largest white dwarf that we have seen has a mass of $m_{\rm tot}=(1.37\pm 0.01)M_{\odot}$. Because the Chandrasekhar limit cannot be below the heaviest white dwarf that we have seen, we must have $m_{\rm tot}\lesssim M_{\rm Ch}$, which in our theory translates to

[TABLE]

B.2 Brown dwarfs

Next, we move on to consider constraints coming from the burning process in brown dwarfs and red dwarfs Burrows:1992fg . This bound is based on the fact that the luminosity generated in the interior of the star from hydrogen burning, $L_{\rm HB}$, must equal the luminosity emitted by the star from the photosphere, $L_{e}$. This analysis gives a smallest mass, the minimum main-sequence mass $M_{\rm MS}$, that is consistent with having a stable burning process in the body of the star. In this work, we use the results of Burrows:1992fg , but keep track of the dependence on $G^{\rm in}_{\Phi}$: the interested reader can find many more details in Burrows:1992fg . (Reference Sakstein:2015zoa has used a similar argument to constrain modification of the Newtonian potential proportional to $m^{\prime\prime}$.)

To find the luminosity in hydrogen burning, we start with the equation for hydrostatic equilibrium $dP/dr=-G^{\rm in}_{\Phi}m(r)\rho(r)/r^{2}$, where $P$ is the pressure, $\rho$ is the mass density, and $m(r)$ is defined in eq. (21) (although we neglect the expansion of the universe on these small scales). Then, we assume an equation of state $P=K(\eta)\rho^{5/3}$, where $\eta\equiv c_{\eta}\rho^{2/3}/T$ (for $c_{\eta}$ a constant) is a measure of the degeneracy of the electron gas and is constant throughout the star. The numerical values of both $K(\eta)$ and $c_{\eta}$ can be determined from fundamental parameters such as $\hbar$, the electron mass, the hydrogen mass, and the number of baryons per electron. This allows us to find an equation for the density profile, which we write as $\rho(r)=\rho_{c}\theta(r)^{3/2}$. After also defining $\tau=r/r_{*}$ and $r_{*}^{2}=5\rho_{c}^{-1/3}K(\eta)/(8\pi G_{\Phi}^{\rm in})$, we obtain the Lane-Emden equation

[TABLE]

where at $\tau=0$, we have both $\theta=1$ and $d\theta/d\tau=0$.

For the total luminosity in hydrogen burning, we obtain

[TABLE]

where $\alpha=4.82$, and for the luminosity emitted from the surface, we obtain

[TABLE]

where $\kappa_{R}\sim 10^{-2}$ (in units of $\text{ cm}^{2}/\text{gm}$) is the Rosseland mean opacity, which determines the optical depth of the star.

The condition to have stable burning is that $L_{\rm HB}=L_{e}$, which means that,

[TABLE]

where $I(\eta)\equiv\eta^{0.184}(1+\alpha/\eta)^{1.509}$. As discussed in Burrows:1992fg , the function $I(\eta)$ has a minimum value of 2.337 at $\eta=34.7$, so that we have the bound

[TABLE]

Now, the smallest red dwarf that has been measured has a mass of $m_{\rm tot}=(0.093\pm 0.0008)M_{\odot}$ Segransan:2000jq . Thus, we finally have

[TABLE]

Notice that the bound depends on the mean opacity $\kappa_{R}$. Previous results in the literature Burrows:1992fg ; Sakstein:2015zoa have set $\kappa_{R}=10^{-2}\text{ cm}^{2}/\text{gm}$ and noted the relatively weak dependence as a justification. However, even though the dependence is quite weak, it can have a fairly large impact on the final bound. For example, in Freedman:2007cm , we see that for $\rho_{e}\sim 10^{-5}-10^{-4}\text{ gm}/\text{cm}^{2}$ and $T_{e}\sim 2000\text{ K}$, that $10^{-2}\lesssim\kappa_{R}\lesssim 10^{-1}$. Thus, we obtain a range of bounds,

[TABLE]

B.3 White dwarfs

Finally, we move on to bounds set by the mass-radius relation of white dwarfs by comparing to the catalogue of stars in Holberg:2012pu . To simplify our analysis, we assume that the stars are each made of a low temperature, completely degenerate Fermi gas, and follow shapiro_teukolsky_1983 . Practically speaking, this assumption means that the profile of the star is not significantly affected by the temperature, an assumption which is reasonable for $T\lesssim 30,000\text{ K}$ Holberg:2012pu . For this reason, we omit the data point with $T=49,000\text{ K}$ in Holberg:2012pu . See Saltas:2018mxc for an example of a more involved analysis which includes these non-trivial temperature effects.

Under these assumptions, the equation of state of the star is given by,

[TABLE]

where $P_{0}=1.4218\times 10^{25}\text{dyne}\,\,\text{cm}^{-2}$, $\rho_{0}=1.9479\times 10^{6}\text{gm}\,\,\text{cm}^{-3}$, and $s\equiv p_{\rm F}/(m_{e}c)$ is the unitless Fermi momentum of an electron, which is in general a function of the distance $r$ from the center of the star, i.e. $s=s(r)$. The equation for hydrostatic equilibrium is the same as that mentioned above eq. (84), and combined with the mass conservation equation, $m^{\prime}(r)=4\pi r^{2}\rho(r)$, one obtains a system of first order differential equations for $s(r)$ and $m(r)$.

After defining $q\equiv r/R_{\odot}$, $\vartheta\equiv m/m_{0}$, $m_{0}\equiv(4\pi R_{\odot}^{3}\rho_{0}/3)$, $\gamma\equiv P_{0}/(4\pi G_{\Phi}^{\rm in}\rho_{0}^{2}R_{\odot}^{2}/3)$ and $y(q)\equiv s(q)^{2}$, this system can be written as

[TABLE]

The initial conditions are $\vartheta(0)=0$ and $y(0)=y_{0}$ (which determines the total mass of the star), and the radius of the star $R$ is given by $y(R)=0$. In particular, we have $\gamma=0.00277$, and $M_{\odot}/m_{0}=7.24\times 10^{-7}$.

To do the fit to the data, we first find the points $\tilde{R}_{i}$ on the theoretical curve that match best with each data point $(R_{i},M_{i})$. In particular, we define

[TABLE]

where $M_{i}$, $\sigma_{M,i}$, $R_{i}$, and $\sigma_{R,i}$ are respectively the mass, error bar for the mass, radius, and error bar for the radius of the $i$th star, and we minimize each $\chi^{2}_{i}$ to find the $\tilde{R}_{i}$ where it is minimum. The total $\chi^{2}$ is then given by

[TABLE]

where for us $N=12$. We then minimize eq. (93) with respect to $\epsilon_{G}$ to find our constraints.

Calling $\tilde{\epsilon}_{G}$ the minimum, we then determine the $n\sigma$ region by finding where $\Delta\chi^{2}(\epsilon_{G})\equiv\chi^{2}(\epsilon_{G})-\chi^{2}(\tilde{\epsilon}_{G})=n^{2}$. In particular, we find a best fit value of $\epsilon_{G}=-0.017$, a $1\sigma$ range of $-0.039\lesssim\epsilon_{G}\lesssim 0.0066$, and a $2\sigma$ range of $-0.060\lesssim\epsilon_{G}\lesssim 0.031$. We present our results in Fig. 2.