The two-body potential of Vainshtein screened theories

TL;DR

This paper investigates the two-body potential in Vainshtein screened theories, revealing that nonlinear effects can violate the Equivalence Principle and impact celestial orbits.

Contribution

It provides the first analysis of two-body potentials in Vainshtein screening beyond spherical symmetry using an Effective One-Body approach.

Findings

Vainshtein screening suppresses fifth forces in dense environments.

Nonlinear dependence on mass ratio leads to Equivalence Principle violation.

Potential impact on lunar orbit dynamics.

Abstract

Adding a light scalar degree of freedom to General Relativity often induces a fifth force whose magnitude is strongly constrained by laboratory experiments and solar system tests. The Vainshtein screening mechanism ensures that the effects of this supplementary force are suppressed in dense environments. However, the field solution of theories exhibiting Vainshtein screening is only known in spherically symmetric situations. In this article we examine in different configurations the two-body potential energy of pointlike particles in a specific $P(X)$ theory with Vainshtein screening. We use ideas borrowed from the Effective One-Body approach of Buonanno and Damour in order to restrict the form of the solution. Our results indicate that, even if Vainshtein screening is also fully active in the equal-mass case, the nonlinear dependance of the two-body energy on the mass ratio implies a violation of the Equivalence Principle. We compute the contribution of this effect to the Moon orbit for generic theories equipped with Vainshtein screening.

Equations

$$S=\int d^{4}x\left[-\frac{(\partial \varphi {)}^2}{2}-\frac{1}{4{\Lambda }^4}(\partial \varphi {)}^4+\frac{\varphi T}{M_P}\right],$$

$$T=-m_1{\delta }^3(\mathbf{x}-{\mathbf{x}}_1)-m_2{\delta }^3(\mathbf{x}-{\mathbf{x}}_2).$$

$$\tilde{\varphi }=\frac{\varphi }{{\Lambda }^2},\tilde{m}=\frac{m}{4\pi M_P{\Lambda }^2},$$

$$S={\Lambda }^4\int dt{d}^{3}x\left[-\frac{1}{2}(\nabla \tilde{\varphi }{)}^2-\frac{1}{4}(\nabla \tilde{\varphi }{)}^4+\tilde{\varphi }\tilde{T}\right],$$

$${\partial }_i\left({\partial }_i\varphi +(\partial \varphi {)}^2{\partial }_i\varphi \right)=-T.$$

$${\varphi }^{\prime }+({\varphi }^{\prime }{)}^3=\frac{M}{r^2}.$$

$${\varphi }_M(r)=-\frac{M}{r}{}_3{F}_2\left(\frac{1}{4},\frac{1}{3},\frac{2}{3};\frac{5}{4},\frac{3}{2};-\frac{27M^2}{4r^4}\right),$$

$$\begin{array}{rl}{\varphi }_M& =-\frac{M}{r}+\frac{M^3}{5r^5}+\dots ,r>r_{\ast },\\ & =C+3{\left(Mr\right)}^{1/3}+\dots ,r<r_{\ast },\end{array}$$

$$r_{\ast }={\left(\frac{27}{4}\right)}^{1/4}\sqrt{M},$$

$$e^{i{S}_{\mathrm{cl}}[{\mathbf{x}}_1,{\mathbf{x}}_2]}=\int \mathcal{D}[\varphi ]e^{iS[{\mathbf{x}}_1,{\mathbf{x}}_2,\varphi ]}.$$

$$\frac{E}{4\pi {\Lambda }^4}=-\frac{m_1{m}_2}{r}+\frac{m_1{m}_2(m_1^2+m_2^2)}{5r^5}+\dots ,$$

$$\tilde{E}=\frac{E}{{\Lambda }^4},$$

$$⟨T{\varphi }_{\mathbf{k}}(t){\varphi }_{\mathbf{q}}(t^{\prime })⟩=-i(2\pi {)}^3{\delta }^3(\mathbf{k}+\mathbf{q})\delta (t-t^{\prime })\frac{1}{k^2}.$$

$$\frac{E_{\mathrm{Newt}}}{4\pi }=-\frac{m_1{m}_2}{r}.$$

$$\begin{array}{rl}\mathrm{Fig}\text{ref}subfig:correctio{n}_{k}mouflag{e}_{4}a& =i\frac{(4\pi {)}^4}{3}{m}_1^3{m}_2{\int }_{\mathbf{K},{\mathbf{k}}_1,{\mathbf{k}}_2}{e}^{i\mathbf{K}\cdot \mathbf{r}}\frac{1}{{\mathbf{k}}_1^2{\mathbf{k}}_2^2{\mathbf{K}}^2(\mathbf{K}+{\mathbf{k}}_1+{\mathbf{k}}_2{)}^2}\\ & \times \left({\mathbf{k}}_1\cdot {\mathbf{k}}_2(\mathbf{K}+{\mathbf{k}}_1+{\mathbf{k}}_2)\cdot \mathbf{K}+2\mathrm{perm}\right),\end{array}$$

$${\int }_{\mathbf{K}}{\mathbf{K}}^2{e}^{i\mathbf{K}\cdot \mathbf{r}},$$

$$\mathrm{Fig}\text{ref}subfig:correctio{n}_{k}mouflag{e}_{4}a=-4\pi i\frac{m_1^3{m}_2}{5r^5}.$$

$$\frac{E_{\mathrm{pp}}}{4\pi }=-\frac{m_1{m}_2}{r}+\frac{m_1{m}_2^3}{5r^5}+\dots$$

$$P=2N+1,$$

$${\left(\frac{m_1}{m_2}\right)}^{q-1},$$

$$\frac{E_{pp}}{4\pi }=m_1\sum _{N\ge 0}{\alpha }_N\frac{m_2^{2N+1}}{r^{4N+1}},$$

$$\frac{E}{4\pi }=\sum _{N\ge 0}{\beta }_N\frac{m_1{m}_2(m_1^{2N}+m_2^{2N})}{r^{4N+1}},$$

$$\frac{E}{4\pi }=\frac{m_1{m}_2}{r}+m_1{\varphi }_{m_2}(r)+m_2{\varphi }_{m_1}(r).$$

$$\begin{array}{rl}\mu & =\frac{m_1{m}_2}{m_1+m_2},\\ M& =m_1+m_2,\\ x& =\frac{m_1}{m_1+m_2}.\end{array}$$

$$\frac{1}{2}{m}_1{v}_1^2+\frac{1}{2}{m}_2{v}_2^2=\frac{1}{2}\mu v^2,$$

$$\frac{E}{4\pi }=\mu \sum _{N\ge 0}{\beta }_N\frac{M^{2N+1}\left(x^{2N}+(1-x{)}^{2N}\right)}{r^{4N+1}}.$$

$$\begin{array}{rl}{\tilde{\alpha }}_N& ={\alpha }_N\left(x^{2N}+(1-x{)}^{2N}\right),N>0,\\ {\tilde{\alpha }}_0& ={\alpha }_0.\end{array}$$

$$r_{\ast }={\left(\frac{27}{4}\right)}^{1/4}\sqrt{M},$$

$$\begin{array}{rl}\frac{E_{\mathrm{tm}}}{4\pi }& =\mu {\varphi }_M(r)\\ & =\mu \sum _{N\ge 0}{\alpha }_N\frac{M^{2N+1}}{r^{4N+1}},\end{array}$$

$$\begin{array}{rl}\frac{E}{E_{\mathrm{tm}}}& =f\left(\frac{E_{\mathrm{tm}}}{E_N}-1\right)\\ & =a_0+a_1\left(\frac{E_{\mathrm{tm}}}{E_N}-1\right)+a_2{\left(\frac{E_{\mathrm{tm}}}{E_N}-1\right)}^2+\dots \end{array}$$

Full text

I Introduction

The possibility of adding a new degree of freedom to General Relativity is often motivated by the cosmic acceleration. Horndeski Horndeski (1974) and its recent extensions GLPV and DHOST theories Gleyzes et al. (2015); Langlois and Noui (2016) are prime examples of such alternative explanations to the cosmological constant. However, these theories rely on a screening mechanism in order to hide the fifth force effects of the scalar field on solar system scales. Among these mechanisms, Vainshtein or K-mouflage screening Vainshtein (1972); Babichev et al. (2009); Babichev and Deffayet (2013) is due to the fact that classical nonlinearities are dominant on scales smaller than the nonlinear Vainshtein radius (the quantum dominated regime, where loop corrections should be taken into account, happens at a much smaller scale de Rham and Ribeiro (2014)).

In this nonlinear regime, the exact field generated by massive bodies is known only in the spherically symmetric case. Generically, it is found that the field goes as $\phi\sim r^{n}$ with $n>-1$ so that it is subdominant compared to the Newtonian gravitational potential $\Phi\sim r^{-1}$. However, tests of GR on solar system scales are very precise and provide useful bounds on these theories using the anomalous perihelion precession that they predict even if the field is screened. Moreover, the recent detection of a binary neutron star merger GW170817 both in the gravitational and electromagnetic channels Collaboration and Collaboration (2017) has constrained the speed of gravitational waves to be that of light, killing a large set of modified gravity models Creminelli and Vernizzi (2017); Ezquiaga and Zumalacarregui (2017); Baker et al. (2017) and constraining the Vainshtein mechanism Sakstein and Jain (2017); Crisostomi and Koyama (2018); Langlois et al. (2018).

Focussing on the solar system for now, the simplest approximation that one can do is to treat planets as test-masses in the (screened) field generated by the sun. This implies that there is no observable effect on the precession of the perihelion of Mercury for the quartic galileon but that the cubic Galileon is already constrained by observations Andrews et al. (2013); Dvali et al. (2003); Iorio (2012). However, for comparable mass bodies the test-mass approximation is not good and this could affect the motion of bodies in a way that depend of their composition, thus violating the Equivalence Principle Belikov and Hu (2013); Hiramatsu et al. (2013). In this case, a analytic form of the two-body potential in the screened regime is missing; it could be particularly relevant for the Earth-moon system, where the Earth-Moon distance is measured with millimeter accuracy by Lunar Laser Ranging experiments Williams et al. (2012); Battat et al. (2009).

Another important physical situation where the two-body potential could be of interest (at least in principle) is the case of two gravitationally bound compact objects (neutron stars or black holes) of comparable mass whose gravitational wave emission is now directly detected at the LIGO/Virgo interferometers. A simple order-of-magnitude estimate shows that the nonlinear Vainshtein radius for solar mass objects is $r_{*,\;1}\sim 10^{15}$ meters and $r_{*,\;2}\sim 10^{19}$ meters for two typical theories that incorporate Vainshtein screening, namely K-Essence Armendariz-Picon et al. (2000, 2001) and Galileons Nicolis et al. (2009); Deffayet et al. (2009) respectively. As the correction to any quantity like the Newtonian energy is suppressed by powers of $r/r_{*}$ where $r$ is the typical size of the system, this makes the effect of the scalar field on such a binary system suppressed by powers of respectively $10^{-10}$ and $10^{-14}$ for the two theories considered. As the gravitational wave phase is measured with an accuracy of $10^{-4}$ in detectors, this hides efficiently the effect of the scalar field from observations. However, the comparable mass case has not been studied in detail yet.

Futhermore, recent work on the radiation of binary systems in Galileon theories, both theoretical Chu and Trodden (2013); de Rham et al. (2013a, b) and numerical Dar et al. (2018), have shown that the radiation itself is screened with powers of $\lambda/r_{*}$, where $\lambda$ is the wavelength of the emitted radiation which is greater by an amount $1/v$ than the size of the system $r$ (here $v$ is the typical velocity of the bodies). The theoretical calculations, which are performed assuming a background field generated by a central mass $M=m_{1}+m_{2}$, are better suited for a small mass ratio inspiral, and it would be interesting to investigate on the case of nearly equal mass ratio. While we will only consider the simpler case of the conservative dynamics in this article, we find that it would be an interesting direction to continue this work.

Here we present a first step towards the study of comparable mass objects in Vainshtein screened theories, namely an approximation to the non-relativistic energy of two bodies. We first approximate the two-body energy outside of the screening radius where traditional calculations based on an expansion in terms of the interactions are reliable. We then try to obtain information on the energy in the screened region based on the knowledge that we gained from the outside. In order to restrict the form of the energy inside the screening region, we use ideas similar to the Effective One-Body (EOB) approach of Damour and Buonanno Buonanno and Damour (1999). The EOB approach is the relativistic generalization of the well-known fact that the two-body problem in Newtonian mechanics is solved by relating it to the motion of a test particle of mass $\mu=\frac{m_{1}m_{2}}{m_{1}+m_{2}}$ in an external field generated by a body of mass $M=m_{1}+m_{2}$. The generalization to GR is that the two-body motion can be approximated by the motion of a test particle of mass $\mu$ in an external deformed Schwartzchild metric, the deformation parameter being the symmetric mass ratio $\nu=\frac{\mu}{M}$. This allows to provide an accurate description of the dynamics of binary systems in regimes that are inaccessible to other analytical approaches, such as Post-Newtonian or Gravitational Self-Force computations Damour (2014). Concerning modified gravity, the EOB formalism has been generalized to scalar-tensor theories in Julie (2018); Julie and Deruelle (2017).

We propose to apply this idea to a two-body system in a Vainshtein screened theory. Since it maps the two-body problem into a spherically symmetric one-body problem, it allows to recover the energy as a function of the known exact field solution. The difficulty is then to identify the parameters of this deformed exact solution. Outside the screening radius, a perturbative (Feynman) expansion allows us to identify the EOB parameters. However, the perturbative expansion breaks down for objects separated by less than the screening radius.

It is however easy to identify the EOB parameters using a simple numerical simulation, and the rest of this article will be devoted to this task. Numerical simulations of screened theories in the quasi-static limit have mostly focused on N-body simulations for cosmological applications Barreira et al. (2013); Li et al. (2012, 2013); Khoury and Wyman (2009); Schmidt (2009, 2009), which is not our concern here. We thus implement our own code using the Finite element solver Fenics Alnaes et al. (2015). We first obtain the energy outside the screening radius in order to confirm our preliminary results. We then concentrate on the two-body energy in the fully screened situation which is relevant for astrophysical systems, and for arbitrary mass ratios. The most important result of this paper is the final mass ratio dependance of the two-body energy, and it is shown in Figure 6. We will finally derive the effect on the orbit of the Moon around the Earth of such a violation of the Equivalence Principle in Sec. V.

Let us now be more precise about the screening theory that we consider. We will focus on a typical theory that incorporate this effect, namely a specific form of K-Mouflage or $P(X)$ Babichev et al. (2009)

[TABLE]

with two static point-like sources,

[TABLE]

Note that we have concentrated on the scalar part of the action, the gravitational action being the one of standard GR. It would be easy to generalize ours results to other well-known theories that incorporate Vainshtein screening, e.g the cubic Galileon Nicolis et al. (2009).

We stress that the chosen form for the interaction in the K-Mouflage case (with a negative sign in front of the quartic term) has a speed of sound around a cosmological background greater than one. Ref Barreira et al. (2015) discusses the conditions needed in order to have a viable k-Mouflage theory. Nonetheless, we chose the theory (1) to test our method because it is one of the simplest settings in which there is Vainshtein screening and the numerical implementation with finite elements is easier since its (weak form) equation of motion involves only first derivatives.

This paper is organized as follows. In Sec. II we compute the two-body energy in the weakly coupled regime where standard perturbative techniques apply. After resumming a class of Feynman diagrams, we recast the problem into an EOB framework in Section III and discuss the form of the EOB energy map in the strongly coupled regime. Section IV then presents the numerical solution to the two-body problem in different regimes. We analyze the consequences of the violation of the Equivalence Principle that it implies on the Earth-Moon-Sun system in Sec. V, before concluding in Sec. VI. We use natural units where $\hbar=c=1$.

II Two-body energy outside the screening radius

In this Section we will concentrate on the action (1). We will only consider static point-sources, so we will ignore time from now on and focus only on the potential energy $E$. In order to have simpler expression in what follows and also to compare directly our results to the numerical simulation, we will introduce rescaled variables

[TABLE]

so that the action of the system writes as

[TABLE]

where $\tilde{T}=-4\pi\tilde{m}_{1}\delta^{3}(\mathbf{x}-\mathbf{x}_{1})-4\pi\tilde{m}_{2}\delta^{3}(\mathbf{x}-\mathbf{x}_{2})$, and we will ignore the tilde from now on. Be careful that now the dimension of $m$ is (in natural units) $-2$ and the dimension of $\phi$ is $-1$.

II.1 Spherically symmetric case

Consider the action (4) with a single point-like source $T=-4\pi M\delta^{3}(\mathbf{x})$. The equation of motion that follows from it reads

[TABLE]

Using spherical symmetry and integrating over a sphere, one can reduce it to a single ordinary differential equation

[TABLE]

This can be solved exactly and writes, for a field vanishing at infinity,

[TABLE]

where ${}_{p}F_{q}$ is the generalized hypergeometric function.

The solution has two regimes separated by the nonlinear scale $r_{*}$

[TABLE]

where $C\simeq-3.7\sqrt{M}$ is a constant of integration (we chose the constant such that the field vanishes at infinity, so it cannot also vanish in zero), and the nonlinear scale is given by

[TABLE]

and corresponds to the radius of convergence of the two series written above.

These two regimes can be seen as expansions in $r_{*}/r$ and $r/r_{*}$ respectively (in fact, it is even possible to reformulate the initial action with additional fields in order to make the screened regime appear from the beginning, see Gabadadze et al. (2012)). We have expanded up to next-to-leading order in the $r>r_{*}$ case because it will prove useful in the following. Choosing $\Lambda$ to be responsible for the cosmic acceleration, $\Lambda^{2}=M_{P}H_{0}$ Barreira et al. (2015), one has $r_{*}\sim 0.1$ Pc for the sun, thus rendering the next order in the $r<r_{*}$ case very subdominant concerning Solar System experiments.

II.2 Two-body problem

Let’s now take a two-body source $T=-4\pi m_{1}\delta^{3}(\mathbf{x}-\mathbf{x}_{1})-4\pi m_{2}\delta^{3}(\mathbf{x}-\mathbf{x}_{2})$. The salient feature of the two-body problem in screened theories is that one can not compute the energy by superposing one-body energies, as one usually does in Newtonian gravity. Said differently, one can not recast the problem in terms of a simple ODE, and instead one should solve a nonlinear PDE of 2 variables, thus explaining the lack of analytical results. Nonetheless, the problem is well formulated outside the screening radius where the nonlinear term can be treated as an interaction. Let us now imagine to have two bodies outside of their respective screening radius and calculate the first nonlinear correction to the potential energy of the two objects. To this aim, we will use Feynman rules derived from the action (4). The potential energy of the system is $\int dtE=-S_{cl}$ where the classical action is obtained as the saddle-point of the path integral,

[TABLE]

Before we start, let us give directly the main result of this Section : the gravitational energy of the two bodies is

[TABLE]

where $r=|\mathbf{x}_{1}-\mathbf{x}_{2}|$. We have factored out a multiplicative coefficient of $\Lambda^{4}$ coming directly from eq. 4: for convenience we will also use a rescaled energy defined by

[TABLE]

and forget the tilde from now on.

Let us now derive this result. Introducing the Fourier modes of the field $\phi=\int_{\mathbf{k}}\phi_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{r}}$ (where $\int_{\mathbf{k}}$ means $\int\frac{d^{3}k}{(2\pi)^{3}}$), the non-relativistic propagator of the field is given by

[TABLE]

To obtain $S_{cl}$, we have to calculate all connected (once we remove the source) Feynman diagrams without internal loops, treating $(\partial\phi)^{4}$ and $T$ as interactions (loops would represent quantum corrections to this potential and are suppressed by powers of $\hbar/L$ where $L$ is the total angular momentum of the system Goldberger and Rothstein (2006), so we will completely ignore them). The full Feynman rules of a simple conformally coupled scalar field in the presence of gravity are given in Ref. Kuntz et al. (2019). At lowest order, there is the one-scalar exchange of Fig 1 that gives rise to the Newtonian potential between the two sources,

[TABLE]

In order to compute the first nonlinear correction to it, we have to calculate the two diagrams of Figure 2. The first one can be put in the following form

[TABLE]

where $\mathbf{r}=\mathbf{x}_{1}-\mathbf{x}_{2}$. The remaining integrals over momentum should be computed with dimensional regularization. To this aim a set of useful integrals in $d=3-2\epsilon$ are given in Appendix A. One can simplify the calculation with the following observation : after having integrated $k_{1}$ and $k_{2}$, one is left with an integral over $K$ which, for dimensional reasons (the first correction to the two-body potential is proportional to $r^{-5}$, see (8)), is of the form

[TABLE]

(there could also be a factor of $\epsilon$ in the power of $\mathbf{K}^{2}$, but it does not change the validity of the argument). Then using formula (71), one can see that there is a pole of the Gamma function in the denominator. Consequently, this integral will vanish in dim. reg (it will be proportional to $\epsilon$) unless there is also a pole of the Gamma function in the numerator. Keeping the only pole that appear in the numerator, and repeatedly using the formulaes given in Appendix A, one can find that

[TABLE]

One can use the same machinery to calculate the second diagram 2b. However, in this case there is no pole of the Gamma function at the numerator, and consequently this diagram vanishes in dim. reg. Including the symmetric counterpart of 2a, we finally obtain the formula for the first correction to the two-body energy given in eq. (11).

II.3 Resumming the test-mass diagrams

The result of Feynman integrals concerning diagram 2a should come as no surprise. Indeed, as we will now explain, it should come back to the energy of a point-particle mass in case of one of the masses goes to zero. In the following, we will assume without loss of generality that $m_{1}<m_{2}$.

Consider the limit $m_{1}\rightarrow 0$. Then the two-body energy should reduce to the energy of a point-particle $m_{1}$ in the external field generated by $m_{2}$, which is (the energy being obtained from the field $\phi_{\mathrm{cl}}$ of eq. (8) via $E_{\mathrm{pp}}=-\int d^{3}x\phi_{\mathrm{cl}}T$)

[TABLE]

From there we see that the result of diagram 2a, which ultimately gives the numerical prefactor in front of the first-order correction to the two-body energy, could not have been otherwise. Had the diagram of figure 2b been nonzero, its value would not have been fixed by this observation, because it is proportional to $m_{1}^{2}m_{2}^{2}$, which vanish (compared to the test-mass diagram proportional to $m_{1}$) in the test-mass limit.

Based on this observation, we propose to resum a particular class of diagram which share the same property at any order in the nonlinear expansion, which we call test-mass diagrams. They are given by a single coupling to the mass $m_{1}$ and any number of couplings to $m_{2}$, as illustrated in Figure 3. Denoting by $N$ the expansion order (i.e, the number of vertex corresponding to the insertion of the nonlinear operator in the diagrams) and by $P$ the number of $m_{2}$ mass insertions, we have that there are $1+4N+P$ field insertions in this diagram. Since there are no loops and the diagram should be connected once we remove the particles wordlines, there are $N-1$ ’internal’ propagators, i.e propagators that connect two nonlinear vertices. Finally, the total number of propagators is $P+1+(N-1)$. Since the total number of field insertions is two times the number of propagators, we get the relation

[TABLE]

so that at order $N$ the mass coefficient of this graph (and its symmetric counterpart) is $m_{1}m_{2}(m_{1}^{2N}+m_{2}^{2N})$.

One can even argue that, at a given perturbation order $N$, this test-mass diagram is the leading one away from the test-mass limit. Indeed, the contribution of other kinds of diagrams, with more $m_{1}$ insertions (we found the only other second-order diagram to vanish, but we see no reason why it should be the case at higher orders), would be of the form $(m_{1}m_{2})^{q}(m_{1}^{2(N+1-q)}+m_{2}^{2(N+1-q)})$ where $2\leq q\leq N+1$. For $m_{2}>m_{1}$ and in the large $N$ limit, the ratio of this quantity to the test-mass diagram is

[TABLE]

which is less than one. However, this is not accounting for the fact that there can be a large number of these other diagrams, which can make them count as much as the test-mass one. Anyway, we shall content ourselves with having understood at least a part of the nonlinear energy.

Now the exact energy of a point-particle $m_{1}$ in the field generated by $m_{2}$ is

[TABLE]

where $\alpha_{N}$ are numerical coefficients that can be easily found by solving eq. (6). By writing in a identical way the contribution of the test-mass graphs to the two-body energy,

[TABLE]

we see that in order to have the good test-mass limit one should impose $\beta_{N}=\alpha_{N}$ for $N\geq 1$, and $\beta_{0}=\frac{\alpha_{0}}{2}=-\frac{1}{2}$. Denoting by $\phi_{m}$ the (exact) spherically symmetric field generated by a body of mass $m$ in eq. (8), one finally finds for the contribution of test-mass graphs to the two-body energy

[TABLE]

This result can be intuitively understood as being the symmetric sum of one-body energies, plus a compensating term that ensures that the Newtonian limit $r\rightarrow\infty$ (where $\phi_{m}\sim-m/r$) is correct. We will compare this analytical resummed energy to the numerical solution in Sec. IV.3.

III Effective One-Body approach

III.1 Energy map outside

Very much like in GR, where the motion of a two-body system (expanded in powers of $r_{s}/r$, where $r_{s}$ is the Schwartzchild radius of the combined mass $M=m_{1}+m_{2}$) can be recast in the motion of a test-mass in a modified external Schwartzchild metric Buonanno and Damour (1999), the above formula for the two-body energy can be expressed into the energy of such a point-particle in a modified external field. Of course, what will now play the role of the Schwartzchild radius is the nonlinear radius $r_{*}$. Let us define, in top of the reduced mass $\mu$ and the total mass $M$ which are the two masses naturally associated to the effective problem, the mass ratio $x$ as

[TABLE]

The kinetic energy of the two objects is easy to rewrite in terms of an effective kinetic energy since we have the well-known relation

[TABLE]

where $\mathbf{v}=\mathbf{v}_{1}-\mathbf{v}_{2}$, and we have set the center-of-mass to the origin of coordinates (this center-of-mass definition would be modified by relativistic corrections, but we do not consider these in this treatment).

As for the approximate potential energy that we obtained above (22), it can be rewritten in terms of the effective parameters as

[TABLE]

We now see that, outside the nonlinear radius, the motion of a two-body system separated by $r=|\mathbf{x}_{1}-\mathbf{x}_{2}|$ can be identified with the motion of a test-particle of mass $\mu$ (at a distance $r$ from the origin) in a modified external field created by $M$, whose modified coefficients in the nonlinear expansion are given by

[TABLE]

What is the nonlinear radius associated to the effective problem ? The perturbative expansion (27) that we wrote above breaks down at the nonlinear radius associated to the biggest of the two masses, i.e $r_{*}=M(1-x)$. However this result is really tied to the test-mass resummation, and calculating more precisely the two-body energy could change it. It is nonetheless natural to assume that the nonlinear radius of the full two-body problem is the one associated to the total mass $M$, i.e.,

[TABLE]

knowing that the real two-body Vainshtein radius (defined as the radius of convergence of the energy expansion for $r\rightarrow\infty$) could differ by numerical factors dependant on $x$.

Alternativaly, we can formulate the equivalent one-body problem in a different way which will prove useful when investigating the behavior inside the nonlinear radius. Using the energy of a test-mass $\mu$ in an external field generated by $M=m_{1}+m_{2}$,

[TABLE]

we can build an energy map between the real energy $E$ and the effective test-mass energy $E_{\mathrm{tm}}$ as follows

[TABLE]

Here $E_{N}=-\mu M/r$ is the Newtonian reference energy, and the function $f$ has been Taylor expanded for small values of the dimensionless ratio $E_{\mathrm{tm}}/E_{N}-1$. Indeed, from eq. (8) this ratio can be expanded outside the nonlinear radius as

[TABLE]

To construct such an energy map, one can choose each value $a_{N}$ such that each coefficient in front of $(M^{2}/r^{4})^{N}$ of eq. (30) matches. Since the small ratio $E_{\mathrm{tm}}/E_{N}-1$ is chosen such that $a_{N}\left(E_{\mathrm{tm}}/E_{N}-1\right)^{N}$ contributes only at order $(M^{2}/r^{4})^{N}$ or higher, this procedure yields an unambiguous value for $a_{N}$ for all $N$.

This energy map proves very useful because it allows to resum the nonlinear behavior into the small parameter $E_{\mathrm{tm}}/E_{N}-1$. We refer the reader to Ref. Buonanno and Damour (1999) for its derivation in the context of the Post-Newtonian formalism.

III.2 Energy map inside

Having understood the behavior of the energy outside the nonlinear radius, we would like now to generalize to astrophysical situations of interest where the two bodies lie deep within their nonlinear radius. The above formula (23) for the two-body energy, even if it resums part of the nonlinear corrections, has no chance to be valid inside $r_{*}$ because it contains the Newtonian reference energy (14) valid only at large radius. However, one can still obtain the coefficients of an energy map from numerical simulations. In App. B we adopt an analytical approach that attempts to relate the two energy maps by a matching condition at the nonlinear radius. While we will argue that this part should yield a qualitative result concerning the two-body energy, we will see when comparing to the numerical simulation that even this qualitative result does not compare well to the real two-body energy. Further improvement is needed in order to obtain a sensible analytical result.

In the following, we will assume that the two-body energy for $r<r_{*}$ can be related, in a spirit similar to the one for $r>r_{*}$, to the energy of a test-mass $\mu$ in an external field generated by the total mass $M$ through an energy map. We have no possibility to calculate directly the modified coefficients as we did in the last Section, but we can obtain them numerically as we will do in Sec. IV.

Inside the nonlinear radius, the energy map should take the form

[TABLE]

A few comments are required here. First, in order to avoid the appearance of an unphysical (mass-dependant) constant in the energy (which we calculated before by assuming its vanishing at infinity, so it cannot also vanish in zero), we chose to write the energy map using the $r$-derivative of the energy $E^{\prime}$. Second, we use a reference energy level equal to the small-$r$ value of the energy since the Newtonian energy is irrelevant inside the Vainshtein radius,

[TABLE]

Finally, note that the small parameter $E^{\prime}_{\mathrm{tm}}/E^{\prime}_{\mathrm{ref}}-1$ is this time expanded as

[TABLE]

where the coefficients $\gamma_{N}$ can be found by solving exactly the nonlinear equation (6) (we will not need their exact expression here).

In the test-mass limit, the real two-body energy should be approximated by the point-particle energy, and consequently $b_{0}=1$ and $b_{1}$, $\dots$ $b_{N}=0$. Now, away from the test-mass limit, we are only interested by the behavior of $b_{0}$ as a function of the mass ratio $x$, since we recall that the next order is further Vainshtein suppressed and thus irrelevant for astrophysical systems. We will now turn to a numerical implementation that will allow us to obtain the coefficient $b_{0}$.

IV Numerical simulation

In this Section, we will directly solve the nonlinear PDE in the two-body case using a finite element solver and obtain the two-body energy in order to compare it to the analytical result.

IV.1 Setup

Using the action (4), one has the following equation of motion for the scalar field

[TABLE]

where we recall that

[TABLE]

To numerically solve this equation, a finite element (FEM) solver 111We use the Python FEM solver fenics Alnaes et al. (2015) is well adapted to the problem, since the PDE can easily be put into a weak form : for any test function $v$ that vanish on the boundary of the integration domain,

[TABLE]

FEM solvers solve the weak form equation by decomposing the unknown $\phi$ on some basis functions $\psi_{j}$, here chosen to be the continuous Lagrange polynomials of second order on the grid chosen to discretize the problem. The solver then finds the coefficients $c_{j}$ of this decomposition $\phi=\sum_{j}c_{j}\psi_{j}$ by evaluating the weak form equation (37) with the test function $v$ being one of the basis functions $\psi_{i}$. This produces a matrix equation for the unknown vector of coefficients $(c_{j})_{j\geq 0}$ that is solved by an efficient sparse LU decomposition. The non-linear term is dealt with Newton iterations, i.e. by setting $\phi=\phi_{0}+\phi_{1}$ with $\phi_{0}$ a function that is close to the solution sought after, linearizing over $\phi_{1}$ then solving for it, and finally iterate the procedure until a desired convergence threshold has been reached.

In cylindrical coordinates (so that the problem becomes effectively two-dimensional), we choose the two bodies to lie along the $z$ axis at positions $+a/2$ and $-a/2$. We regularize the delta-functions by replacing them with Gaussians,

[TABLE]

where we have taken care of the fact that the $r$ variable goes from [math] to $\infty$ while $z$ goes from $-\infty$ to $\infty$.

There are two scales involved in this problem, the separation between the two bodies that we denote by $a$ and the nonlinear scale $r_{*}=\sqrt{M}$. We choose for the domain of integration the half-disk defined by $r>0$, $R\leq R_{\mathrm{max}}$ where $R=\sqrt{r^{2}+z^{2}}$. The domain is automatically discretized by the FEM solver with a resolution of approximately $64\times 64$ points, with a manual refinement of the grid near the two bodies. The boundary conditions are chosen such that $\partial_{r}\phi(0,z)=0$ (as required by the symmetry of the problem) and $\phi(r,z)=-\frac{m_{1}+m_{2}}{R}$ for $R=R_{\mathrm{max}}$. The second boundary condition corresponds to recovering both spherical symmetry and the Newtonian-like behavior of the field far from the two bodies, where we know the exact solution which is given by eq. (8). For it to be consistent, we must also ensure $r_{*}\ll R_{\mathrm{max}}$.

Once we get the field solution, the energy can be computed as

[TABLE]

A remarkable fact to be noted is that the self-energy contribution $4\pi m_{1}\phi(\mathbf{x}_{1})+4\pi m_{2}\phi(\mathbf{x}_{2})$ is not divergent, contrary to the Newtonian case. This is due to the fact that the field goes to a constant as $\phi\sim|\mathbf{r}-\mathbf{x}_{\alpha}|^{1/3}$ close to the source $\alpha$, instead of diverging as $\frac{1}{|\mathbf{r}-\mathbf{x}_{\alpha}|}$. The same goes for the integral over all space. We consequently do not need to renormalize the energy.

Finally, after having compared the behavior of the energy for $r\gtrsim r_{*}$ to the theoretical predictions, we will need the energy in the realistic $r\ll r_{*}$ case. In this setup, we can ignore the quadratic term in the equation of motion (37). This corresponds to take an infinite Vainshtein radius. Correspondingly, we have to change the exterior boundary condition to $\phi(r,z)=3((m_{1}+m_{2})R)^{1/3}$ for $R=R_{\mathrm{max}}$ because the field is screened over all space. The energy is computed ignoring also the quadratic term, but there is a subtelty in its definition because the integral is formally divergent as $R_{\mathrm{max}}^{1/3}$ (a consequence of the good UV and bad IR behavior of the field). This divergence is the same as the one associated to a single particle of mass $M=m_{1}+m_{2}$ sitting at the origin, and so to renormalize the theory we susbtract to the energy the contribution of such a field which is

[TABLE]

IV.2 Numerical tests

In order to assess the validity of our numerical scheme, we have performed two numerical tests. Before presenting them, let us first discuss the choice of the numerical parameters $R_{\mathrm{max}}$ and $\sigma$. A finite choice of $R_{\mathrm{max}}$ brings a correction of the order $a/R_{\mathrm{max}}$ to the field (where we recall that $a$ is the separation between the two bodies), and so we choose $R_{\mathrm{max}}=50a$ in our simulation.

Concerning the effect of $\sigma$, let us consider a single point-particle of mass $m_{2}$ in the $r\ll r_{*}$ case with boundary condition $\phi(r,z)=3(m_{2}R)^{1/3}$ for $R=R_{\mathrm{max}}$. In this case, if the particle was really pointlike, the field would vanish at the origin. This is no longer the case if the particle has a finite extent $\sigma$. Rather, the central value of the field is

[TABLE]

simply by dimensional analysis. But this could potentially be problematic in the calculation of the two-body energy which requires the evaluation of the field at the particle location, see eq. (39). The contribution of this term in the total two-body energy is

[TABLE]

this equality being true up to a numerical factor.

In order to avoid an unphysical dependance on $\sigma$ in the energy, we have to tune $\sigma$ in order that the term in eq. (41) is much smaller than the term in eq. (42). This gives the following bound:

[TABLE]

where we recall that $x=m_{1}/(m_{1}+m_{2})$ is the mass ratio. We choose $\sigma=10^{-6}a$, which is consistent with the minimal value of $x$, $x_{\mathrm{min}}=0.03$ for $a=1$, that we will use.

Having chosen a value for the numerical parameters, we have first checked that, for a single particle at the origin of the coordinates, the field solution corresponds to the exact solution given in eq. (7). The result is plotted in Figure 4, where we see that there is perfect agreement between the theoretical and numerical values.

The second nontrivial check that we have performed is to verify that the two-body energy in the fully screened regime (i.e, neglecting the quadratic operator as explained in Sec. IV.1) indeed varies as $E\propto a^{1/3}$. In Figure 4 we have plotted the numerical value of the two-body energy in the $r\ll r_{*}$ case in logarithmic plot, from which we immediately confirm the expected behavior.

IV.3 Results

Figure 5 presents the numerical two-body in the $r\gtrsim r_{*}$ regime against different theoretical predictions for equal masses. The value of the energy when the charges are taken infinitely far apart is not zero: it can be easily calculated as two times the value of the energy when we plug the exact one-body field solution (7) into the action. This yields

[TABLE]

We can see that the numerical value of the energy for large separation matches very well this result, thus providing again a strong check of the validity of our code. The resummed energy provides an improvement over the simple Newtonian potential, although the gain is $\sim 50$ % at the nonlinear radius.

Concerning the fully screened situation $r\ll r_{*}$, Figure 6 presents the effective parameter $b_{0}$ as a function of the mass ratio $x$, where $b_{0}$ is defined as

[TABLE]

where we chose $a=1$ to obtain our results ($b_{0}$ does not depend on $a$, as emphasised in Sec. IV.2) and the numerical energy $E$ is computed by keeping only the nonlinear term in the action, as explained in Sec. IV.1. This is the most important result of this work, since it presents the two-body energy in astrophysically relevant situations and for arbitrary mass ratios. We can see that $b_{0}\simeq 0.75$ in the equal-mass case, which means that screening is a bit more efficient than when there is a large mass hierarchy (for which $b_{0}=1$). For convenience, a fit to $b_{0}$ with a sixth-order polynomial gives

[TABLE]

The dependance of the two-body energy on the mass ratio implies a direct violation of the Weak Equivalence Principle (WEP). In the next Section we will explore the consequences of this result on the orbit of the Moon, which is measured with a great accuracy by Lunar Laser Ranging.

V Violation of the Weak Equivalence Principle

In this Section we will show how our results imply an Equivalence Principle violation that would be visible on the Moon orbit. We will make heavy use of the following mass ratios

[TABLE]

where $m_{\odot}$ is the Sun mass, $m_{\oplus}$ is the Earth mass and $m_{M}$ is the Moon mass.

V.1 Three-body system, and finite size corrections

In this part we will examine the applicability of our results to a three body system like the one formed by the Sun, the Earth and the Moon (herafter, SEM system). Fifth forces generated by the scalar interaction are expected to induce a supplementary perihelion precession that would be visible on planetary orbits Talmadge et al. (1988). But the lunar perihelion precession cannot be computed by ignoring the scalar field generated by the Sun, and using a perturbative treatment Ref. Andrews et al. (2013) reached the conclusion that the perturbations blows up once distances hierarchies are taken into account. We will recover their result using a different approach.

To start with, let us reformulate the action of a two-body system in a different way. Labelling the two objects as $1$ and $2$, we split the scalar field according to

[TABLE]

where $r_{2}=|\mathbf{x}-\mathbf{x}_{2}|$ is the distance to the body $2$, and $\phi_{2}$ is the spherically symmetric field (8) generated by the same body. In the fully screened situation that is of interest to us, this fields writes as

[TABLE]

If $m_{2}$ is much greater than the other mass, we can think of $\psi$ being a fluctuation on top of the dominant field generated by the mass $m_{2}$, but for now let’s keep the discussion general and not assume any mass hierarchy (we still assume $m_{1}\leq m_{2}$ by symmetry). Then by inserting this decomposition in the action (4) (where we ignored the quadratic term that is subdominant on small scales), we get the following action

[TABLE]

where $S[\phi_{2}]$ is the original action (4) applied to $\phi_{2}$, and $T_{1}=-4\pi m_{1}\delta^{3}(\mathbf{x}-\mathbf{x}_{1})$ is the source term corresponding to the object $1$. The term linear in $\psi$ vanishes because of the equations of motion for $\phi_{2}$.

We now ask the question : close to the object $1$, what is the behavior of the fluctuation $\psi$ ? In other words, is there an operator that dominates the action for $\psi$ in (50) ? It seems natural to assume that close enough to the source $1$, we recover the behavior $\psi\sim(m_{1}r_{1})^{1/3}$, which means that the last nonlinear operator in eq. (50) dominates. Let us assume this is the case and derive the condition on $r_{1}$ for this to be true.

If the term $(\nabla\psi)^{4}$ dominates in the action, then we recover the same spherically symmetric action as for the second object, and we consequently obtain

[TABLE]

Let us now make the ratio between the operator $(\nabla\psi)^{4}$ and another one in the action, say the term cubic in $\psi$. Using Cauchy-Schwarz inequality, we get

[TABLE]

Using the same scaling, one can show that the term quadratic in $\psi$ is similarly suppressed with respect to the term cubic in $\psi$. This result is quite important. It means that around the first body, the field can be well approximated simply by taking the linear superposition $\phi_{2}+\psi$ of two spherically symmetric solutions. This is true up to the maximal distance to the first body

[TABLE]

which depends on the mass ratio. Since the Vainshtein radius of a massive body is $r_{*,\alpha}=\sqrt{m_{\alpha}}$, this equation can also be interpreted as

[TABLE]

which simply means that the distance to each body is measured in units of its Vainshtein radius.

This has several consequences. First, it means that we can ignore the finite-size of the bodies and treat them as point-particles as long as their radius is less than $r_{1}^{\mathrm{max}}$. In this particular theory, this is true both for the Sun-Earth system and for the Earth-Moon system if we take the numerical values of their respective radius and masses.

Second, it means that we can simply add the fifth forces felt by a satellite that orbits sufficiently close to a planet itslef orbiting around its star. Unfortunately, for the SEM system and the particular screening theory considered here this is not true, as

[TABLE]

where we have taken $r_{1}$ to be the Earth-Moon distance and $r_{2}$ to be the Sun-Earth distance. The end result is that we expect the Lunar perihelion precession rate to be corrected by an amount depending on the masses and distances hierarchies.

For another type of Vainshtein screening such as the Galileon-3 Nicolis et al. (2009), the same reasoning shows that we can superpose the nonlinear solutions provided

[TABLE]

which is true at the 10% level for the SEM system. The same conclusion was also reached in Nicolis and Rattazzi (2004) with similar arguments. This shows that in the case of a Galileon-3, the Lunar perihelion precession can be calculated by simply ignoring the Sun 222Note that we can also apply this reasoning to the Sun embedded inside our own galaxy. In this case, $m_{2}/m_{1}=m_{\mathrm{Gal}}/m_{\odot}\simeq 10^{12}$, $r_{1}/r_{2}=r_{\mathrm{SE}}/r_{\mathrm{Gal}}\simeq 6\times 10^{-10}$ where $m_{\mathrm{Gal}}$ is the galactic mass, $r_{\mathrm{Gal}}$ is the distance to the Galactic centre, and we have made the simplifying assumption that the mass of the Galaxy is concentrated at its centre. Eqs 52 and 56 then show that, for both Galileons and $P(X)$, we can indeed neglect the background field generated by the galaxy for the inner planets of the Solar System. I am grateful to Philippe Brax for pointing this out. Since the Earth-Moon mass ratio is $x\simeq 10^{-2}$, the calculation can be carried out in the test-mass approximation and yields an interesting constraint on the size of a Galileon-3 operator Dvali et al. (2003).

V.2 Weak Equivalence Principle violation

The fact that the two-body energy of two massive particles is not the one of a reduced mass $\mu$ in an external field implies a violation of the Weak Equivalence Principle, as we will now show. In this Section we will assume a general Vainshtein screening mechanism that gives rise to a fifth force (not necessarily the specific $P(X)$ theory that we studied in the rest of the article). If we first neglect the moon, the total (gravitational and scalar) interacting Lagrangian of the Earth and the Sun writes as

[TABLE]

Here $r_{SE}=|\mathbf{y}_{S}-\mathbf{y}_{E}|$ is the Earth-Sun distance, $x_{SE}=m_{\oplus}/(m_{\odot}+m_{\oplus})$ is the Earth-Sun mass ratio, $r_{*}$ is the Vainshtein radius of the Sun, $n$ is an exponent that depends on the type of screening considered, and $\alpha$ is an unknown function of the coefficient $x_{SE}$ that can be found numerically as in Sec. IV. In order to avoid confusion between positions and mass ratios, we denote in this Section the positions by the letter $\mathbf{y}$. This type of Lagrangian is common to all theories endorsed with Vainshtein screening, with the expressions of $n$, $r_{*}$ and $\alpha$ different among theories. For example, $n=4/3$ for the $P(X)$ theory considered above in this article, and $n=3/2$ for a Galileon-3.

We will make the supplementary assumption that we can get the total force felt by the Moon by simply adding the fifth forces of the Earth and the Sun. As we showed in Sec. V.1, this is the case for a Galileon-3 but not for the particular $P(X)$ example examined in the rest of the article. A complete treatment of this case would necessitate further work. With this assumption we can write the total (non-relativistic) interaction Lagrangian of this three-body system (to get the total Lagrangian from this, one should also add the kinetic energies from the three bodies)

[TABLE]

where $S$ designates the Sun, $E$ the Earth, $M$ the Moon, and each line of this equation is the two-body Lagrangian of eq. (57) adapted to each pair of bodies. We will now derive the fifth force incidence on the Lunar motion along the lines of Ref. Damour and Vokrouhlicky (1996).

The first line of eq. (58) gives rise to the Earth anomalous perihelion precession, and the third line to the Lunar anomalous perihelion precession. These effects are already discussed in other References Dvali et al. (2003); Iorio (2012); Andrews et al. (2013) and we will not comment on them. From now on, we will ignore the last line of eq. (58) which does not give rise to the leading order equivalence principle violation that we are going to derive. Let us expand the distances to the Sun around the Earth-Moon center-of-mass, which is defined with the usual expression

[TABLE]

where $\mathbf{y}_{E}$ and $\mathbf{y}_{M}$ are the positions of the Earth and the Moon respectively. Then the distances to the Sun can be expressed by

[TABLE]

where $\mathbf{y}_{S}$ is the Sun position. Expanding the total Lagrangian to first order in $r_{EM}$, one finds

[TABLE]

where $r=|\mathbf{y}_{S}-\mathbf{Y}|$ is the center-of-mass distance to the Sun, and as discussed above we have dropped the third line of eq. (58). Upon introducing the reduced mass $\mu_{EM}=m_{\oplus}m_{M}/(m_{\oplus}+m_{M})$, one finds the following expression

[TABLE]

As discussed in Ref. Damour (2014), there are two physical effects that stem from this Lagrangian 333In fact Ref. Damour (2014) discusses three physical consequences, but the third one is not apparent at the expansion order we work with, and is expected not to be measureable.. The first line of eq. (62) implies that the gravitational constant involved in the motion of the Earth-Moon system around the Sun is not the same than the one involved in the motion of this system around its barycenter (third line of eq. (58)). However, this effect has practically no observable consequences.

The second physical effect, on which we will concentrate, is the perturbation on the lunar orbit implied by the second line of eq. (62). Nordtvedt Nordtvedt (1968) showed in 1968, in the context of Scalar-Tensor theories, that this term implies a modulation of the Lunar orbit with amplitude

[TABLE]

In scalar-tensor theories, $\delta_{AB}$ is the fractional variation of Newton’s constant due to the gravitational self-energy of each body ; said equivalently, Newton’s constant appearing in the gravitational attraction between two bodies is

[TABLE]

Although the physical origin is quite different, the equivalence principle violation considered here gives rise to the same term in the second line of the Lagrangian (62) than the equivalence principle violation of Scalar-Tensor theories. The parameter $\delta_{\oplus}-\delta_{M}$ is replaced by the following quantity

[TABLE]

In the second line we have expanded $\alpha$ to first order in the mass ratios, $\alpha(x)\simeq\alpha_{0}+\alpha_{1}x$ and neglected the Moon mass ratio compared to the Earth mass ratio, $x_{SM}\ll x_{SE}$.

Current LLR data give the constraint $|\delta_{\oplus}-\delta_{M}|\lesssim 10^{-13}$ Williams et al. (2012), which by ignoring the $\mathcal{O}(1)$ factor of $(1-n)$ yield the constraint

[TABLE]

Let us contrast this with the anomalous perihelion constraint. For the earth, the precession is measured with an accuracy of $10^{-11}$ (for a precise constraint on a Galileon-3, see Ref. Iorio (2012)), and the anomalous perihelion precession in typical Vainshtein screened theories is proportional to the ratio of scalar and gravitational potentials Dvali et al. (2003), which gives the constraint

[TABLE]

On the one hand, we gain two orders of magnitude by using Lunar Laser Ranging data, but on the other hand the equivalence principle violation is suppressed by the mass ratio $x_{SE}\simeq 10^{-6}$ with respect to the perihelion bound, resulting in a looser constraint if we assume that $\alpha_{0}\sim\alpha_{1}\sim\mathcal{O}(1)$.

VI Conclusions

In this paper we have analyzed for the first time the two-body potential energy of pointlike objects in Vainshtein screened theories for arbitrary mass ratios. One the one hand, from outside the nonlinear radius, the problem is amenable to a perturbative treatment which we use to resum a class of Feynman graphs. We derive an Effective One-Body energy map which relates the two-body energy to the one of a test particle in an external field. On the other hand, we conjecture the existence of such an energy map inside the nonlinear radius where the nonlinear screening term dominates the action. We have tried to get the analytical behavior of this expansion by a matching procedure at the nonlinear scale. Improving the accuracy of this matching procedure would necessitate to know the exterior potential energy with a higher accuracy. This could be done by calculating Feynman diagrams with two insertions of the nonlinear operator, instead of only one as we did in Sec. II.2. This would necessitate more involved computations that we leave for future work.

At the same time, we have performed a numerical simulation in order to get the most relevant effective coefficient $b_{0}$ corresponding to the ratio between the real two-body energy and the energy of a test-mass in an external field. In the test-mass limit, $b_{0}\simeq 1$ while in the equal-mass limit, we find $b_{0}\simeq 0.75$ in the particular $P(X)$ theory that we considered. This means that Vainshtein screening is active even in the fully nonlinear situation where we do not assume one mass to be smaller than the other, with departure from simple order-of-magnitude estimates encoded into a simple coefficient $b_{0}$ which can be found with a numerical simulation. While we have focused on a particular model exhibiting nonlinearities, we expect such a feature to be valid in any model endorsed with Vainshtein screening such as Galileons.

As we showed in Sec. V, the fact that the two-body energy differs from the test-mass one implies a violation of the weak equivalence principle, as the Earth and the Moon would not fall the same way towards the Sun. We used this fact to bound the size of the coupling parameters of any theory relying on Vainshtein screening to hide the effects of a fifth force in the Solar system. Although the final constraint is looser than the one obtained from perihelion precession for a Galileon-3, the methodology employed is quite general and we intend to use it to investigate on the case of two neutron stars or black holes in their inspiral phase.

Acknowledgements.

I would like to thank Federico Piazza and Julien Bel for insightful comments on the draft.

Appendix A Useful integrals

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix B Matching interior and exterior solutions

In this Appendix, we will attempt to obtain a qualitative analytical behavior of the energy as a function of the mass ratio inside the nonlinear radius. We can use the formula for the energy outside the screening radius that we found in Sec. II.3 as a limiting boundary condition for the energy map inside the screening radius. Since the energy map partly resums the nonlinear behavior, it should allow to get a qualitative description of the coefficient $b_{0}$.

Using the energy derived in eq. (23), one finds that the left-hand side of the energy map (32) writes as

[TABLE]

for $r$ close to $r_{*}$. We recall that $x$ is the mass ratio $x=m_{1}/M\leq 1/2$ and that $\phi_{m}$ is the spherically symmetric field (8) generated by a mass $m$. Note that, as $\frac{\phi^{\prime}_{Mx}}{x}\rightarrow\frac{M}{r^{2}}$ as $x\rightarrow 0$, this expression has the good test-mass limit, and in this limit we recover that $b_{0}=1$ and $b_{1},\;\dots\;b_{N}=0$.

If the energy map really resums a part of the nonlinear dynamics, then it should give an asymptotic expansion near to the Vainshtein radius, as suggested by the relative smallness of the expansion parameter at $r_{*}$,

[TABLE]

where we recall that $r_{*}=\left(27/4\right)^{1/4}\sqrt{M}$.

Consequently, we can try to find a first approximation for the value of $b_{0}$ by neglecting higher-order contributions in the right-hand side of eq. (32). This gives for $b_{0}$

[TABLE]

This first prediction for $b_{0}$ is plotted in Figure 7.

One can try to improve the result by matching not only the numerical value of the energy map (32) at $r=r_{*}$, but also its derivatives with respect to $r$. This will allow us to extract $b_{1}$, $b_{2}$… while at the same time refining our prediction on $b_{0}$. More precisely, matching up to the $k^{\mathrm{th}}$ derivative will allow us to determine $k+1$ coefficients in the energy map, say $b_{0}^{(k)},\;b_{1}^{(k)},\;\dots\;b_{k}^{(k)}$.

When should we stop this procedure ? Obviously it cannot be pushed to arbitrary order because the boundary condition on the energy is not sufficiently precise. A sensible criterion can be found by thinking in terms of asymptotic series. Let us write the $k^{\mathrm{th}}$ prediction for $b_{0}$ as

[TABLE]

where we have defined $b_{0}^{(-1)}=0$.

In this equation we have highlighted the fact that each step of the iteration brings out an additional factor of $b_{0}^{(l)}-b_{0}^{(l-1)}$ to $b_{0}$. Then the asymptotic series criterion requires us to cut the sum at the term at which $|b_{0}^{(l)}-b_{0}^{(l-1)}|$ is minimized in order to get the best precision on $b_{0}$. In Figure 8 we have plotted the ratio

[TABLE]

for $l=1,\;2,\;3$. The minimal term of the series is the one for which $r_{l}<1$ and $r_{l+1}>1$. We can see directly that, while the two first iterations of the procedure (i.e., matching up to the second derivative) seem to bring an improvement in the determination of $b_{0}$, the third one gives a large correction.

The end result of this Section is shown in Figure 7, where we plot the first four predictions $b_{0}^{(0)},\;b_{0}^{(1)},\;b_{0}^{(2)}$ and $b_{0}^{(3)}$ for $b_{0}$. Of course, an immediate drawback of this method is that we are not able to assess quantitatively the precision of the estimation of this coefficient ; moreover, it is apparent from Figure 7 that even the qualitative behavior is not correct. A refinement of the method is needed but we leave this for future work.