On the four-loop static contribution to the gravitational interaction potential of two point masses

TL;DR

This paper calculates specific four-loop contributions to the gravitational potential between two point masses at the fourth post-Newtonian order, revealing a rational coefficient and confirming previous theoretical results.

Contribution

It provides an analytical computation of a four-loop scalar integral and corrects prior effective-field-theory calculations, enhancing understanding of the 4PN gravitational interaction.

Findings

The total four-loop contribution has a rational coefficient due to cancellation effects.

The generalized Riesz formula enables analytical evaluation of complex integrals.

The results confirm the current understanding of the 4PN effective gravitational action.

Abstract

We compute a subset of three, velocity-independent four-loop (and fourth post-Newtonian) contributions to the harmonic-coordinates effective action of a gravitationally interacting system of two point-masses. We find that, after summing the three terms, the coefficient of the total contribution is rational, due to a remarkable cancellation between the various occurrences of $\pi^2$. This result, obtained by a classical field-theory calculation, corrects the recent effective-field-theory-based calculation by Foffa et al. [arXiv:1612.00482]. Besides showing the usefulness of the saddle-point approach to the evaluation of the effective action, and of x-space computations, our result brings a further confirmation of the current knowledge of the fourth post-Newtonian effective action. We also show how the use of the generalized Riesz formula [Phys. Rev. D 57, 7274 (1998)] allows one to analytically compute a certain four-loop scalar master integral (represented by a four-spoked wheel diagram) which was, so far, only numerically computed.

Equations

$$M\equiv m_1+m_2;\mu \equiv \frac{m_1{m}_2}{m_1+m_2};\nu \equiv \frac{\mu }{M}=\frac{m_1{m}_2}{(m_1+m_2{)}^2}.$$

$${\hat{H}}_{\text{4PN}}\left(\frac{{\mathbf{r}}_{12}}{GM},\frac{\mathbf{p}}{\mu }\right)$$

$$+{\nu }^3{\hat{H}}_3^{\text{4PN}}+{\nu }^4{\hat{H}}_4^{\text{4PN}}.$$

$$L_{\le \text{4PN}}^{\mathrm{int}}[{\mathbf{x}}_a,{\mathbf{v}}_a,{\dot{\mathbf{v}}}_a,\cdots ]$$

$$\sim \frac{G{m}_1{m}_2}{r_{12}}\sum _{n\le 4}{\left(\frac{Gm}{r{c}^2}+\frac{v^2}{c^2}+\frac{r_{12}\dot{v}}{c^2}+\cdots \right)}^n$$

$$\sim \frac{G{m}_1{m}_2}{r_{12}}\sum {\left(\frac{Gm}{r{c}^2}\right)}^{n_1}{\left(\frac{v^2}{c^2}\right)}^{n_2}{\left(\frac{r_{12}\dot{v}}{c^2}\right)}^{n_3}\cdots .$$

$$L_{{\text{4PN}}}^{\rm int}\big{|}^{O(G^{5})}_{{\rm linear}\ {\rm in}\ \dot{\bm{v}}_{a}}\subseteq L_{\leq{\text{4PN}}}^{\rm int}\big{|}^{(n=4)}+\sum_{n=0,1,2}L_{\leq{\text{4PN}}}^{\rm int}\big{|}^{(n)}_{\dot{v}^{2}},$$

$$\frac{105}{32}{\pi }^2\frac{G^5(m_1^4{m}_2^2+m_1^2{m}_2^4)}{c^8{r}_{12}^5}-\frac{71}{16}{\pi }^2\frac{G^5{m}_1^3{m}_2^3}{c^8{r}_{12}^5};$$

$$L_{33}^{\text{FMSS}}$$

$$L_{49}^{\text{FMSS}}$$

$$L_{50}^{\text{FMSS}}$$

$$-\frac{32}{3}{\pi }^2\frac{G^5{m}_1^3{m}_2^3}{c^8{r}_{12}^5}$$

$$S_{\mathrm{eff}}^{\mathrm{class}}[x_a(s_a)]$$

$$+\frac{1}{2}\sum _{a,b}{e}_a{e}_b\iint \mathrm{d}{x}_a^{\mu }\mathrm{d}{x}_{b\mu }\delta \left((x_a-x_b{)}^2\right).$$

$$e^{\frac{i}{\hslash }{S}_{\mathrm{eff}}^{\mathrm{quant}}}=\int D{A}_{\mu }{e}^{\frac{i}{\hslash }(S_{\mathrm{particle}}+S_{\mathrm{field}})}.$$

$$S_{\mathrm{eff}}^{\mathrm{class}}[x_a(s_a)]={\left[S_{\mathrm{pm}}+S_{\mathrm{EH}}+S_{\mathrm{gf}}\right]}_{g_{\mu \nu }(x)\to g_{\mu \nu }^{\mathrm{gf}}[x_a(s_a)]},$$

$$e^{\frac{i}{\hslash }{S}_{\mathrm{eff}}^{\mathrm{quant}}}=\int D{g}_{\mu \nu }{e}^{\frac{i}{\hslash }(S_{\mathrm{pm}}+S_{\mathrm{EH}}+S_{\mathrm{gf}})}.$$

$$S_{\mathrm{eff}}^{\mathrm{quant}}[x_a(s_a)]=S_{\mathrm{eff}}^{\mathrm{class}}[x_a(s_a)]+O(\hslash ).$$

$$S_{\mathrm{eff}}$$

$$+{\left[\frac{1}{3}ϵ+\frac{1}{2}Z+F+\frac{1}{2}H+\frac{1}{4}X\right]}_{G^3{m}^4}+\cdots .$$

$$\mathcal{G}(t,\mathbf{x};t^{\prime },{\mathbf{x}}^{\prime })\equiv -4\pi {\left(\Delta -\frac{1}{c^2}{\partial }_t^2\right)}^{-1}$$

$$=-4\pi \left({\Delta }^{-1}+\frac{1}{c^2}{\Delta }^{-2}{\partial }_t^2+\frac{1}{c^4}{\Delta }^{-3}{\partial }_t^4+\dots \right)\delta (t-t^{\prime }),$$

$$e^{\frac{i}{\hslash }{S}_{\mathrm{eff}}}=\int D\phi e^{\frac{i}{\hslash }(\int [\frac{1}{2}\phi \mathcal{K}\phi +\phi s+g{\phi }^3+\cdots ]}.$$

$$⟨\phi (x)\phi (y)⟩=\int D\phi e^{\frac{i}{\hslash }\int [\frac{1}{2}\phi \mathcal{K}\phi ]}\phi (x)\phi (y)=i\hslash {\mathcal{K}}_{x,y}^{-1},$$

$$\int D\phi e^{\frac{i}{\hslash }\int [\frac{1}{2}\phi \mathcal{K}\phi ]}\sum _n\frac{(i/\hslash {)}^n}{n!}{\left(\int (\phi s+g{\phi }^3+\cdots )\right)}^n$$

$$\mathcal{K}\mathcal{G}(x,y)=-\delta (x-y).$$

$$Z[s]=e^{\frac{1}{ϵ}{S}_{\mathrm{eff}}}=\int D\phi e^{\frac{1}{ϵ}(\int [\frac{1}{2}\phi \mathcal{K}\phi +\phi s+g{\phi }^3+\cdots ])},$$

$$Z[s]\approx e^{\int [\frac{1}{2}{\phi }_{\ast }\mathcal{K}{\phi }_{\ast }+{\phi }_{\ast }s+g{\phi }_{\ast }^3+\cdots ]},$$

$$0$$

$$=\mathcal{K}\phi +s+3g{\phi }^2+\cdots .$$

$${\phi }_{\ast }=\mathcal{G}s+\mathcal{G}\left(3g(\mathcal{G}s{)}^2\right)+\cdots ,$$

$$S_{\mathrm{eff}}[s]\approx \ln Z[s{]}^{\mathrm{saddle}}\approx \frac{1}{2}s\mathcal{G}s+g(\mathcal{G}s{)}^3+\cdots .$$

$$Z[s]=\int D\phi e^{\frac{1}{ϵ}(\int \frac{1}{2}\phi \mathcal{K}\phi )}\sum _n\frac{1}{n!}{\left(\int [\phi s+g{\phi }^3+\cdots ]\right)}^n,$$

$$⟨\phi (x)\phi (y)⟩=-{\mathcal{K}}_{x,y}^{-1}=\mathcal{G}(x,y).$$

$$S=S_{\mathrm{pm}}+S_{\mathrm{field}},$$

Full text

I Introduction

The analytical study, to ever-increasing accuracy, of the motion and radiation of two compact bodies (with comparable masses) in General Relativity has been vigorously pursued over the last decades, with the aim of helping the construction of accurate templates for the data-analysis pipeline of the network of ground-based interferometric gravitational-wave detectors. And indeed, the bank of 250 000 templates used in the matched-filter searches and data-analyses of the first observing run of advanced LIGO TheLIGOScientific:2016pea have been defined Taracchini:2013rva within the analytical effective one-body (EOB) formalism Buonanno:1998gg ; Buonanno:2000ef ; Damour:2000we ; Damour:2001tu ; Damour:2008gu . The EOB formalism combines, in a suitably resummed format, perturbative, analytical [post-Newtonian (PN)] results on the motion and radiation of compact binaries, with some non-perturbative information extracted from numerical simulations of coalescing black-hole binaries.

In this work we focus on the conservative dynamics of two spinless bodies. The current level of accuracy on the analytical knowledge of this problem is the fourth post-Newtonian (4PN) accuracy. The 4PN Hamiltonian [in Arnowitt-Deser-Misner (ADM) coordinates] of two mass points111It was shown long ago Damour1983 that the extension effects of compact bodies show up only at the 5PN level, so that they can be modelled by point masses below the 5PN accuracy. is non-local in time, and was first obtained in complete form in Ref.  Damour:2014jta , based on the computation of the local contributions in Ref. Jaranowski:2015lha . (Earlier, partial results were obtained in Refs. Foffa:2012rn ; Jaranowski:2012eb ; Jaranowski:2013lca ; Bini:2013zaa .) The non-local action of Ref.  Damour:2014jta was reduced to a local Hamiltonian in Ref. Damour:2015isa . (This “local reduction” was obtained by using an expansion in powers of the eccentricity, together with suitable redefinitions of the phase-space variables, as detailed in Damour:2016abl .) Since then, the only other attempt to derive the complete 4PN dynamics has been the harmonic-coordinates calculation of Ref. Bernard:2015njp . Most of the terms in the action of Ref. Bernard:2015njp agree with the results of Refs. Damour:2014jta ; Damour:2015isa , except a couple of them.

To discuss the discrepancies between the harmonic-coordinates result of Ref. Bernard:2015njp and the ADM-coordinates one of Refs. Damour:2014jta ; Damour:2015isa ; Jaranowski:2015lha , it is convenient to order the various contributions to the interaction Hamiltonian (which starts by the Newtonian one $-Gm_{1}m_{2}/r_{12}$) by means of the powers of the symmetric mass ratio $\nu$. Our notation (besides using $G$ for Newton’s gravitational constant) is

[TABLE]

We denote the two masses of the binary system as $m_{1}$ and $m_{2}$, while $r_{12}=|{\bm{r}}_{12}|$ (where ${\bm{r}}_{12}\equiv{\bm{x}}_{1}-{\bm{x}}_{2}$) denotes the relative distance. We work here in the center-of-mass system; when doing so in a Hamiltonian framework, one considers the ratio ${\bm{p}}/\mu$ (where ${\bm{p}}={\bm{p}}_{1}=-{\bm{p}}_{2}$) as fixed. The $\mu$-reduced Hamiltonian ${\hat{H}}_{\text{4PN}}\equiv{\hat{H}}_{\text{4PN}}/\mu$ can then be decomposed as

[TABLE]

Here, ${\hat{H}}^{\text{4PN}}_{0}$ describes the 4PN-level contribution to the dynamics of a test mass moving around a central body of mass $M=m_{1}+m_{2}$, while $\nu{\hat{H}}^{\text{4PN}}_{1}$ describes the first self-force (1SF) correction to the latter test-mass dynamics, $\nu^{2}{\hat{H}}^{\text{4PN}}_{2}$ the second self-force (2SF) correction, etc. [In diagrammatic language, computing 1SF effects on the “small mass” $m_{1}$ (say) corresponds to computing one gravitational loop in the external gravitational field of a black hole of mass $m_{2}\gg m_{1}$.]

It was shown in Ref. Bernard:2015njp that all the terms that are non-linear in $\nu$ [i.e. $\nu^{2}{\hat{H}}^{\text{4PN}}_{2}+\nu^{3}{\hat{H}}^{\text{4PN}}_{3}+\nu^{4}{\hat{H}}^{\text{4PN}}_{4}$ in Eq. (I)] in their harmonic-coordinate result agree (modulo a suitable contact transformation) with the ADM action of Ref. Damour:2014jta . The discrepancies are limited to the $\nu$-linear (1SF-level) contribution $\nu{\hat{H}}^{\text{4PN}}_{1}$. It was later shown in Ref. Damour:2016abl that the $\nu$-linear terms in the local reduction Damour:2015isa of the ADM non-local action were in full agreement with several different (analytical and numerical) gravitational self-force computations (combined with results from EOB theory, and from the first law of binary dynamics LeTiec:2011ab ; Blanchet:2012at ; Tiec:2015cxa ), and it was concluded that several claims, and results, of Ref. Bernard:2015njp were incorrect, and must be corrected both by evaluating the energy in keeping with Refs. Damour:2014jta ; Damour:2015isa , and by the addition of a couple of ambiguity parameters linked to subtleties in the regularization of infrared and ultraviolet divergences. The values of the needed additional ambiguity parameters (denoted there $\Delta a$ and $\Delta b$, when using the “gauge” $c=0$) were determined in Damour:2016abl to be $\Delta a^{\rm tot}=\Delta a-11\,\frac{16}{15}\,\Delta C=\frac{2179}{315}$ and $\Delta b=+12\,\frac{16}{15}\,\Delta C=-\frac{192}{35}$ [inserting Eqs. (6.1), (7.4) of Damour:2016abl in Eq. (6.3) there]. Recently, Ref. Bernard:2016wrg confirmed all those conclusions, and notably the values of the ambiguity parameters (which they denote $-\delta_{1}\equiv\Delta a^{\rm tot}$ and $-\delta_{2}\equiv\Delta b$) that must be added to the harmonic-coordinates Hamiltonian to correct it.

Very recently, Ref. Foffa:2016rgu applied the so-called effective field theory (EFT) method Goldberger:2004jt to the computation of a subset of the contributions to the harmonic-coordinates Lagrangian $L$. Given some specified gauge-fixing additional contribution to the Einstein-Hilbert action [here the standard harmonic-coordinates gauge-fixing term $S_{gf}=(16\pi G)^{-1}\int\mathrm{d}^{D}x\sqrt{g}(-\frac{1}{2}g_{\mu\nu}\Gamma^{\mu}\Gamma^{\nu})$ with $\Gamma^{\mu}\equiv g^{\rho\sigma}\Gamma^{\mu}_{\rho\sigma}$], the effective222The reduced action (obtained by “integrating out” the mediating field) describing the conservative dynamics of some particles is called by various names: Fokker action, reduced action, effective action, … . Here, we shall use the name “effective action” to avoid confusion with the “order-reduced” local action Damour:2015isa which replaces the original non-local-in-time 4PN action by an equivalent local-in-time one. action, $S_{\rm eff}=\int\mathrm{d}tL$, describing the conservative dynamics of the binary system can be decomposed in powers of $G$ and of the velocities ${\bm{v}}_{a}$, $a=1,2$ (together with their various time derivatives $\dot{\bm{v}}_{a}$, $\ddot{\bm{v}}_{a}$, $\cdots$)333Here, we formally consider the non-local-in-time piece of the (interaction) action as a functional of the infinite set of time derivatives of ${\bm{v}}_{a}$.. In particular, the structure of the interaction Lagrangian (say up to the 4PN level) is roughly described by expanding the $n$-th power (with $n\leq 4$) on the first rhs of the following sketchy formula (where $c$ denotes the velocity of light):

[TABLE]

In the multiple sum on the last rhs the sum of the powers $n=n_{1}+n_{2}+n_{3}+\cdots$ must be $\leq 4$. As will be described in more detail below, the various contributions in the fully expanded form of $L_{\leq{\textrm{4PN}}}^{\rm int}$ can be described in terms of Feynman diagrams. Here, following Ref. Foffa:2016rgu , we shall focus on the contributions having the highest possible power of $G$, i.e. $n_{1}=4$, $0=n_{2}=n_{3}=\cdots$ in Eq. (I), corresponding to a purely “static” term, quintic in $G$, without effects linked to velocities, or derivatives of velocities.

It was understood long ago Damour:1985mt ; Damour:1990jh that any term that is non-linear in the derivatives of velocities can be eliminated from a higher-order Lagrangian $L(x,v,\dot{v},\cdots)$ by adding suitable “double-zero” terms [quadratic in $\dot{v}-(\dot{v})^{\text{on-shell}}$], thereby allowing one to replace a general higher-order Lagrangian by an equivalent simpler one that is linear in accelerations. (A further reduction, involving a redefinition of the particle variables allows one to eliminate the accelerations Schafer:1984mr ; Damour:1985mt ; Damour:1990jh .) The procedure of reduction of terms quadratic (or more) in accelerations to a linear dependence in accelerations involves the on-shell equations of motion $(\dot{v})^{\text{on-shell}}\sim Gmr_{12}^{-2}(1+O(1/c^{2}))$, and thereby introduces a mixing between the various powers of $G$ in the expanded Lagrangian Eq. (I). In particular, after reduction to a $\dot{v}$-linear form (as was done in Bernard:2015njp ), the contribution proportional to $G^{5}$ is given by a sum of terms coming from terms $\sim G^{1+n}$ in Eq. (I) having $n\leq 4$. More precisely, as terms quadratic in accelerations contain at least two powers of $1/c^{2}$, we have

[TABLE]

with values $n=0,1,2$.

Foffa et al.  pointed out Foffa:2016rgu three facts: (i) the terms non-linear in accelerations coming from $n=0$ and $n=1$ on the rhs of Eq. (4) only contribute rational coefficients to the lhs; (ii) the terms quadratic in accelerations coming from $n=2$ on the rhs contribute the following $\pi^{2}$-dependent 4PN $O(G^{5})$ terms to L_{{\text{4PN}}}^{\rm int}\big{|}^{O(G^{5})}_{{\rm linear}\ {\rm in}\ \dot{\bm{v}}_{a}}

[TABLE]

and, (iii) the $\pi^{2}$-dependent terms (5) [coming from the $O(G^{3}\dot{v}^{2}$) action] coincide with the ($v$- and $\dot{v}$-independent) $\pi^{2}$-dependent terms present in the full, linear-in-acceleration harmonic-coordinates 4PN Lagrangian derived in Ref. Bernard:2015njp [see Eq. (5.6f) there].

As the latter contributions in the harmonic-coordinates Lagrangian of Bernard:2015njp agree with corresponding contributions in the ADM Hamiltonian of Damour:2014jta , one would then conclude (barring a coincidental agreement between two incorrect results) from Eq. (4) that the coefficients entering the $n=4$ [i.e. $O(G^{5})$] contribution to the original (non-linear in derivatives of $v$) 4PN effective Lagrangian L_{\leq{\text{4PN}}}^{\rm int}\big{|}^{(n=4)} should not contain any $\pi^{2}$, i.e. should be a rational number. In other words, there should be no new, genuine $\pi^{2}$ at the $O(G^{5})$ level.

However, Foffa et al. Foffa:2016rgu have recently reported the computation, within the EFT approach, of the 50 Feynman diagrams contributing to the $n=4$ [i.e. $O(G^{5})$] contribution to $L_{\leq{\text{4PN}}}^{\rm int}[{\bm{x}}_{a},{\bm{v}}_{a},\dot{\bm{v}}_{a},\cdots]$ in Eq. (I). Their results comprise three contributions with $\pi^{2}$-dependent coefficients, namely

[TABLE]

Note that we cited here twice the quantitities respectively denoted ${\cal L}_{33}$, ${\cal L}_{49}$ and ${\cal L}_{50}$ in Foffa:2016rgu because it seems that they implicitly assume that the $m_{1}$-$m_{2}$ symmetric Lagrangian contributions ${\cal L}_{33}\sim{\cal L}_{49}\sim{\cal L}_{50}\sim m_{1}^{3}m_{2}^{3}$ should be augmented by their $1\leftrightarrow 2$ images, and thereby doubled.

The sum of the three contributions (6) contains the $\pi^{2}$-dependent term

[TABLE]

which disagrees with the result of Bernard:2015njp (which is derived with the use of the same, harmonic gauge-fixing term). In terms of the $\mu$-reduced Hamiltonian (I), this discrepancy is proportional to $\nu^{2}$, and therefore at the 2SF level. All the contributions $O(\nu^{2})$ to the $\mu$-reduced 4PN action agreed (modulo a contact transformation) between the two existing complete 4PN calculations Damour:2014jta and Bernard:2015njp .

The main aim of the present paper is to perform a new, independent calculation of the three contentious Lagrangian contributions ${L}_{33}$, ${L}_{49}$ and ${L}_{50}$ to decide whether there were subtle, hidden errors in Bernard:2015njp and Damour:2014jta that coincidentally agree, or whether there is an error in the EFT-theory evaluation of the corresponding Feynman integrals. A secondary aim of the present paper concerns the analytical computation of a certain $d$-dimensional, four-loop “master” Feynman integral, denoted ${\cal M}_{3,6}$ in Ref. Foffa:2016rgu . This master integral contributes to the values of both ${L}_{33}$, and ${L}_{50}$. Though they employed some of the most advanced Feynman-integral computation techniques, Foffa et al. did not succeed in analytically evaluating the $d$-dimensional, four-loop integral ${\cal M}_{3,6}$, and had to resort to a many-digit numerical evaluation of the coefficients of the Laurent expansion of ${\cal M}_{3,6}(d)$ in powers of $\varepsilon\equiv d-3$. This evaluation gave very solid numerical evidence for the presence of $\pi^{2}$ at the $\varepsilon^{0}$ level, and this has been assumed to be exactly true in the computation of the results Eqs. (6).

The two main results of the present paper will be: (i) to show that one can analytically evaluate (by notably using the generalized Riesz formula derived in Jaranowski:1997ky , which was also crucial to the computation of the local ADM 4PN Hamiltonian computation Jaranowski:2015lha ) the relevant first three terms in the $\varepsilon$ expansion of the four-loop master integral ${\cal M}_{3,6}(d=3+\varepsilon)$, and, in particular, rigorously prove the presence of $\pi^{2}$ at the $\varepsilon^{0}$ level ; and (ii) explain away the seeming contradiction following from the presence of $\pi^{2}$ in the EFT evaluation (6) of the four-loop integrals ${L}_{33}$, ${L}_{49}$ and ${L}_{50}$, by showing that a new, independent calculation of these integrals (using, instead of the EFT technique of Foffa:2016rgu , the alternative, diagrammatic “field theory” approach to the effective action introduced long ago by Damour and Esposito-Farèse Damour:1995kt , together with ${\bm{x}}$-space techniques, and the use of the generalized Riesz formula), leads to results that crucially differ from the ones cited above in that $\pi^{2}$ simply cancels out in the sum ${L}_{33}+{L}_{49}+{L}_{50}$.

II Various approaches to the effective action for gravitationally interacting point masses

The introduction of a classical “variational principle that takes account of the mutual interaction of multiple particles without introducing fields” dates back to Fokker’s 1929 definition Fokker1929 of the following relativistic functional of several worldlines (labelled by $a,b=1,\ldots,N$) describing $N$ electromagnetically interacting charged point masses (here we use $c=1$, and all the quantitites are defined in a Minkowski spacetime of signature mostly plus)

[TABLE]

The action (II) is obtained by classically “integrating out” the electromagnetic field $A_{\mu}(x)$ in the usual total relativistic action for the particles and the field, i.e. by replacing the (time-symmetric, Lorenz-gauge) solution of the equation of motion of $A_{\mu}(x)$ in presence of given worldlines (say $A_{\mu}^{\rm Lorenz}[x;x_{a}(s_{a})]$ ) in the original particle $+$ field action.

The action (II) played a central role in the 1949 work of Wheeler and Feynman Wheeler:1949hn . Let us also note that Fokker’s original paper features spacetime diagrams of worldlines interacting via time-symmetric propagators. It is therefore probable that the introduction of quantum interaction diagrams (or Feynman diagrams) by Feynman around the same time was partly motivated by Fokker’s classical interaction diagrams. Clear evidence for this is the 1950 paper of Feynman Feynman:1950ir in which he introduces the (complex) quantum effective action for charged particles defined (in modern notation) through taking the logarithm of a functional integral over the field (in presence of given classical charged worldlines)

[TABLE]

He then explicitly shows that $S_{\rm eff}^{\rm quant}$, (9), only differs from its classical counterpart, (II), by the replacement of the (real) time-symmetric propagator $\delta((x_{a}-x_{b})^{2})$ by the (complex) (Stückelberg-)Feynman propagator $\delta_{+}((x_{a}-x_{b})^{2})$, with $\delta_{+}(x)=\frac{i}{\pi(x^{2}+i0)}=\delta(x^{2})+PP\frac{i}{\pi x^{2}}$.

The gravitational analog of the above classical, effective action for the general relativistic interaction of point masses reads InfeldPlebanski ,

[TABLE]

where $S_{\rm pm}=-\sum_{a}\int m_{a}\sqrt{-g_{\mu\nu}(x_{a})\,\mathrm{d}x_{a}^{\mu}\,\mathrm{d}x_{a}^{\nu}}$ denotes the point-mass action, $S_{\rm EH}$ the Einstein-Hilbert action, and $S_{\rm gf}$ a gauge-fixing term, and where $g_{\mu\nu}^{\rm gf}[x_{a}(s_{a})]$ denotes the gauge-fixed solution of Einstein’s equations in presence of given worldlines. The gravitational analog of the above (formally) quantum, effective action reads444Note that this definition is misprinted in Refs. Foffa:2011ub ; Foffa:2016rgu , where the lhs of Eq. (11) is simply written as $\frac{1}{\hbar}S_{\rm eff}^{\rm quant}$, without the exponential, and without the imaginary unit; these omissions being later corrected by considering connected diagrams and by multiplying the rhs by $-i$.

[TABLE]

In the classical limit, one can evaluate the (formal) path integral (11) by the saddle point (or stationary phase) approximation. As the extrema with respect to $g_{\mu\nu}(x)$ of the exponent in (11) are simply classical solutions of the gauge-fixed Einstein equations, one immediately sees that (formally)

[TABLE]

We recalled the above rather well-known facts to clarify that the so-called EFT method [formally based on (11)] computes (in the classical limit, and when considering the conservative555However, as discussed in Goldberger:2004jt and several subsequent papers, the imaginary part of $S_{\rm eff}^{\rm quant}[x_{a}(s_{a})]$ gives useful information about radiation-damping effects. dynamics) exactly the same quantity as the classical, Fokker (or, for that matter) ADM, reduction method (10).

However, the two different definitions of the effective action suggest different technical methods for computing it, and this is where there is a real practical difference in the traditional PN (or post-Minkowskian) computations of $S_{\rm eff}$, and in the EFT-inspired one. First, let us recall that long before the EFT method was set up Goldberger:2004jt , an alternative, diagrammatic “field theory” approach to the (classically defined) effective action was introduced in Ref. Damour:1995kt . It was explicitly shown in Damour:1995kt how the pertubative, post-Minkowskian way of solving the gauge-fixed Einstein’s equations (say in harmonic gauge) leads to a (classical, Feynman-like) diagrammatic expansion of the effective action for the particles of the form (with the normalizations chosen there)

[TABLE]

Here, each letter $I,V,T,\cdots$ is chosen to evoque a correspondingly shaped diagram, when representing the source by, say, $\circ$. For instance, the diagram $I$ denotes the vertical concatenation of two sources, $\circ$ and $\circ$, located at the end points of the $I$, via an intermediate (time-symmetric) gravitational propagator $\arrowvert$, say $\arrowvert_{\!\!\circ}^{\!\!\circ}$. (Here, the propagator is defined as minus the inverse of the kinetic term, see more discussion of this choice below.) In other words, $I$ denotes the one-graviton-exchange diagram of the gravitationally interacting source. [When decomposing the material, two-body source according to the masses, say $\circ=m_{1}\circ_{1}+\,m_{2}\,\circ_{2}$, the $I$ diagram gives three contributions: two self-gravity ones, $O(Gm_{1}^{2})$ and $O(Gm_{2}^{2})$, and a relativistic Newtonian interaction one: $Gm_{1}m_{2}\,\arrowvert_{{\!\!\circ}_{1}}^{{\!\!\circ}_{2}}$.] Similarly, $T$ denotes a diagram where three sources (located at the end points of the $T$) are connected via three gravitational propagators that meet at a cubic vertex in the middle of the upper branch of the $T$. In addition, Ref. Damour:1995kt gave explicit rules for computing the numerical coefficients to be put in front of each diagram to correctly evaluate the effective action666The explicit coefficients shown in Eq. (II) above follow from the specific $S_{n}[\varphi]=\frac{1}{n}V_{n}[\varphi^{n}]$ vertex normalization chosen in Damour:1995kt . When absorbing the conventional prefactor $\frac{1}{n}$ in the definition of the vertex $V_{n}$ many of the factors in the effective action (II) become unity, and the remaining ones are usual symmetry factors.. It is sometimes convenient (to better exhibit the physics contained in the effective action) to draw each individual source $m_{1}\circ_{1}$ or $m_{2}\circ_{2}$ as a spacetime worldline. Then each diagram in the post-Minkowskian expansion (II) becomes made of concatenated propagators, with some propagators starting on the worldlines, and intermediate propagators joining either a gravitational vertex, or a worldline. We shall later give explicit examples of such spacetime representations of effective-action diagrams, which generalize the representation used by Fokker himself back in 1929.

When further taking the PN expansion of the time-symmetric (scalar) propagator, say

[TABLE]

each post-Minkowskian diagram in the expansion (II) will generate a sequence of PN-type diagrams (involving inverse powers of the Laplacian, together with time-derivatives, as propagators). These are now three-dimensional (or $d$-dimensional) diagrams made of PN-propagators $\Delta^{-n}$ connecting the two point masses $m_{1}\delta({\bm{x}}-{\bm{x}}_{1})$ and $m_{2}\delta({\bm{x}}-{\bm{x}}_{2})$ via some intermediate field points that are integrated over. It has been known for a long time that the computation of the effective action at the $n$PN level involves diagrams whose topology features $\leq n$ loops. The topological loops can be recognized either on the spacetime diagrams, or on the projections as $d$-dimensional diagrams. For instance, Fig. 1 in Ref. Damour:2001bu represents a spatial, two-point, three-loop diagram representing a 3PN-level contribution $O(G^{4}m_{1}^{3}m_{2}^{3})$ to the effective ADM action of two point masses. Below, we shall give examples (with four loops, at the 4PN level) of such spatial diagrams.

Summarizing: the usual, Fokker-like computation of the PN-expanded gravitational action (using either harmonic coordinates, or ADM coordinates, and using either traditional methods or the field-theory-diagrammatic technique of Damour:1995kt ) leads to a sum of ${\bm{x}}$-space integrals involving the concatenation of PN-propagators $\Delta^{-n}\partial_{t}^{n+1}\delta(t-t^{\prime})$ and their joining at intermediate spatial points, with vertices involving two derivatives (because of the structure of the gravitational action $\partial\partial hh+\partial\partial hhh+\cdots$).

The main points we wanted to emphasize here about the traditional Fokker-like computation of the effective action are: (i) all the contributions ot the effective action are explicitly real; (ii) all the integrals are in ${\bm{x}}$-space; (iii) all the integrations by parts used to reduce integrals to some “master” integrals are done in ${\bm{x}}$-space; (iv) at each stage of the calculation one keeps track of the numerical coefficients multiplying each integral, because they are directly furnished by the replacement $g_{\mu\nu}(x)\to g_{\mu\nu}^{\rm gf}[x_{a}(s_{a})]$ of the gauge-fixed solution in (essentially) the Einstein-Hilbert Lagrangian (be it in harmonic guise, or in the ADM one).

By contrast, the EFT approach to the effective action is based on expanding functional integrals of the type (here written, for pedagogical purposes, as a scalar toy-model, with a source $s(x)$, taken simply as a linear coupling here),

[TABLE]

Instead of expanding around the saddle point of the exponent (as done in the usual Fokker approach) one expands the functional integral around the Gaussian approximation defined by the free term with kinetic operator ${\cal K}$, and with elementary contraction given by

[TABLE]

where ${\cal K}^{-1}_{x,y}$ denotes the inverse of the kinetic operator (i.e. a Green’s function). One then expands the exponent on the rhs of (15),

[TABLE]

applying Wick’s theorem to compute all the $\varphi$ contractions arising from the various powers $i^{n}\left(\int(\varphi s+g\varphi^{3}+\cdots)\right)^{n}/n!$ coming from the expansion of the exponential. (We henceforth set $\hbar=1$ for simplicity.) The lowest-order contribution comes from the term quadratic in $s$, namely $\frac{i^{2}}{2}\langle\int\mathrm{d}x\,\varphi(x)s(x)\int\mathrm{d}y\,\varphi(y)s(y)\rangle=\frac{i^{3}}{2}\iint\mathrm{d}x\mathrm{d}y\,{\cal K}^{-1}_{x,y}s(x)s(y)$. Factoring one power of $i$ this contributes $\frac{i^{2}}{2}\iint\mathrm{d}x\mathrm{d}y\,{\cal K}^{-1}_{x,y}s(x)s(y)=-\frac{1}{2}\iint\mathrm{d}x\mathrm{d}y\,{\cal K}^{-1}_{x,y}s(x)s(y)$, to the effective action $S_{\rm eff}$. This indeed coincides with the (correctly normalized) one-quantum exchange energy denoted $+\frac{1}{2}I$ above.

Summarizing: the quantum, Feynman-like computation of the PN-expanded gravitational action deals with a sum of Wick contractions from the powers of the interaction terms $\varphi s+g\varphi^{3}+\cdots$ in the original field + particle action. This calculation involves many imaginary units $i$. Because of a certain quantum tradition, these calculations have been done in $\bm{p}$-space, rather than in $\bm{x}$-space, using, e.g., elementary field contractions $\langle\varphi\varphi\rangle=i/(-p^{2})$ if the kinetic term is $\Box$. (We use the mostly plus signature.) In this approach one has to take care of correctly multiplying each diagram by the needed symmetry factor (which can be somewhat tricky when considering high-order contractions). In doing the explicit calculations at the $n$th PN order, there appear diagrams having up to $n$-loops, corresponding to integrating over $n$ independent loop momenta variables. [Note that though the Fourier-space integrals to compute are in one-to-one correspondence (modulo an overall Fourier transform) with the $\bm{x}$-space ones which enter the other approach, the computations are somewhat different, and the number of integrations to perform over intermediate points in the $\bm{x}$-space approach is generally not equal to the number of topological loops in the diagram.]

Let us discuss the equivalence between the two approaches in further detail, and also emphasize why it is useful to define the Green’s function ${\cal G}(x,y)$ associated with the kinetic operator $\int[\frac{1}{2}\varphi{\cal K}\varphi]$ as being minus the inverse of the kinetic term, say

[TABLE]

This was the convention of Damour:1995kt , and it leads, when coupling the field to a source $s(x)$, (i.e. $\int[\frac{1}{2}\varphi{\cal K}\varphi+s\varphi]$) to a leading-order effective action equal to $+\frac{1}{2}\int s(x){\cal G}(x,y)s(y)$. Actually, the usefulness of the minus sign in the Green’s function definition (18) is hidden in the usual “quantum” definition (16) of the elementary contraction of the field $\varphi$. Indeed, the rhs of Eq. (16) is really $-(\frac{i}{\hbar}{\cal K})^{-1}$, i.e. minus the inverse of the operator appearing in the exponent of the (functional) integral that one is dealing with. In other words, the imaginary units $i$ that crowd up the EFT computations are irrelevant. The essential point is that we have two different ways of approximating an integral of the type

[TABLE]

where $\epsilon$ is a formal small parameter, and where the functional measure is normalized so that $Z[s=0]=1$. As the perturbative calculation of $S_{\rm eff}=\epsilon\ln Z[s]$ is a purely algebraic matter, one can replace the quantum “small parameter” $\frac{\hbar}{i}$ by any formally small parameter $\epsilon$. One can even simplify the writing by assuming that the small parameter is absorbed in the definition of the quadratic form $\varphi{\cal K}\varphi$, and of the interaction terms. Doing so, the classical approximation to the integral (19) is to use the saddle-point approximation

[TABLE]

where $\varphi_{*}$ is the saddle point, i.e. the solution of

[TABLE]

In this approach, one solves the saddle point condition (II) by a perturbative series away from the unperturbed solution $\varphi=0$, namely [with ${\cal K}=-{\cal G}^{-1}$ according to the definition (18)]

[TABLE]

where the needed integrations over intermediate spacetime points are left implicit. This leads to an expansion of the effective action in powers of the source $s$:

[TABLE]

In the other, Feynman-like approach one approximates (at the exponential accuracy) the integral (19) by expanding the integrand away from the Gaussian term

[TABLE]

using the elementary contraction

[TABLE]

From the above reasoning, it is guaranteed that this will give the same result, (23), for the logarithm of $Z[s]$. But this reasoning shows that all the $i$’s are a useless complication (which can easily lead to sign errors when there are many of them), as we are computing a real effective action (when using the time-symmetric Green function appropriate to describing the conservative dynamics).

III Explicit expressions of the relevant four-loop, 4PN effective-action contributions

We focus, in this paper, on the few effective-action contributions that Ref. Foffa:2016rgu emphasized as being potentially problematic. As explained in Foffa:2016rgu these terms are purely “static” and follow from the simplified particle + field action

[TABLE]

where the (static) point-mass action is

[TABLE]

and where the field action Kol:2010si ; Foffa:2016rgu is

[TABLE]

Here, we followed the notation of Foffa:2016rgu , apart from the fact that we use $\lambda=1/\Lambda=\sqrt{32\pi G\ell_{0}^{d-3}}$. The gravitational field degrees of freedom are described by $\phi$ and $\sigma_{ij}$, with $\gamma_{ij}=\delta_{ij}+\lambda\sigma_{ij}$. In addition, $c_{d}\equiv 2\frac{d-1}{d-2}$, $\gamma=\det\gamma_{ij}$, $\sigma=\sigma_{ii}$, and $f_{,i}\equiv\partial_{i}f$. Note that, in this approximation, only $\phi$ is directly coupled to the particles. The tensor field $\sigma_{ij}$ is only excited through the cubic vertex following from the kinetic term of $\phi$:

[TABLE]

For the four-loop terms we are interested in, only the linear coupling of $\phi$ to the particles,

[TABLE]

matters. Here the Lagrangian density of the source is

[TABLE]

where $\delta_{a}\equiv\delta({\bm{x}}-{\bm{x}}_{a})$.

We can then describe the algebraic structure of the relevant particle + field Lagrangian as

[TABLE]

where

[TABLE]

includes the kinetic terms, the linear coupling to matter, and the cubic vertex between $\sigma_{ij}$ and $\phi$ coming from the $\phi$ kinetic term (III), namely

[TABLE]

[Note that, following Eq. (18), we have expressed the kinetic operators of $\phi$ and $\sigma$ in terms of the corresponding Green’s functions ${\cal G}_{\phi},{\cal G}_{\sigma}$.] The remaining, higher-order terms in the relevant 4PN action have the algebraic structure

[TABLE]

They respectively correspond to the $O(\sigma^{2})\partial_{i}\phi\partial_{j}\phi$ terms in the $\phi$ kinetic term (III), and to the sum of the last line in the field action (III), and of the terms coming from the kinetic terms of $\sigma_{ij}$ when considering the terms of order $\lambda\sigma_{..}$ in the expansion of

[TABLE]

using $\sqrt{\gamma}=1+\frac{1}{2}\lambda\sigma+O(\lambda^{2})$, $\gamma^{ij}=\delta_{ij}-\lambda\sigma_{ij}+O(\lambda^{2})$. Hence,

[TABLE]

As for the terms $b\,\sigma^{2}\phi^{2}$, they are explicitly given by

[TABLE]

with

[TABLE]

The saddle-point conditions (or field equations of motion) for $\phi$ and $\sigma$ have the structure

[TABLE]

As the solution of these field equations of motion is only needed for being replaced in the Lagrangian ${\cal L}(\phi,\sigma,s)$, it is well-known that it is enough to solve the equations of motion coming from ${\cal L}_{0}$, i.e. to take $\epsilon=0$ in the above field equations. Indeed, as $\delta{\cal L}/\delta\,{\rm field}=0$, the corrections to the field solution coming from $\epsilon{\cal L}_{1}$ contribute only at order $\epsilon^{2}$ to the Fokker action. [It is essentially this basic fact that, upon the suggestion one of us (TD), was used to simplify the recent 4PN harmonic-coordinates computation of the Fokker action Bernard:2016wrg .] To lowest-order in a non-linearity expansion in the source [i.e. in an expansion in powers of the two masses $m_{1},m_{2}$, see (31)], we immediately see that the solutions of the above field equations are

[TABLE]

From the above reasoning, we deduce the first result that the contribution of the action correction $\epsilon{\cal L}_{1}$ to the effective (Fokker) action is simply obtained by replacing in $\epsilon{\cal L}_{1}$ the fields $\phi$ and $\sigma$ by their lowest-order solution (because this is enough to get $\epsilon{\cal L}_{1}$ to order $s^{6}$), namely

[TABLE]

This result takes care of two of the contentious action contributions highlighted by Foffa:2016rgu , namely $L_{33}$, linked to $b\,\sigma^{2}\phi^{2}$, and $L_{50}$, linked to $c\,\sigma^{3}$, and allows one to compute them straightforwardly (including all numerical factors). The remaining contentious action contribution, $L_{49}$, is easily seen (from its diagram in Fig. 1 of Foffa:2016rgu ; see also below) to arise from the exchange of two cubic vertices, (34). Therefore, in a Fokker-type calculation, this term arise from solving the field equations of motion Eqs. (40), (41) to fourth order in the $\phi$-$\sigma$ coupling (34), that we had left in the zeroth order action ${\cal L}_{0}$, (33).

It is fairly easy to solve Eqs. (40), (41) (without the $\epsilon{\cal L}_{1}$ terms) to order $O(a^{4})$. First, let us note that we are talking here about a purely algebraic calculation that could be done by iterating polynomial expressions. The aim of our calculation is to get the correct numerical coefficient in front of the $O(a^{4})$ Fokker action contribution. This can be formally done by solving Eqs. (40), (41) as if $\phi$ and $\sigma$ were ordinary numbers. As Eq. (40) (without the $\epsilon{\cal L}_{1}$ term) is linear in $\phi$ we can solve $\phi$ in terms of $\sigma$ and replace the answer in the second equation. Denoting

[TABLE]

the solution of Eq. (40) reads $\phi(\sigma)=(1-2x)^{-1}{\cal G}_{\phi}s$, and its insertion in Eq. (41) (without the $\epsilon{\cal L}_{1}$ term) reads

[TABLE]

This is easily solved by iteration in powers of the source:

[TABLE]

Inserting this solution in ${\cal L}_{0}$ then easily leads to an expansion in even powers of $s$:

[TABLE]

Here, we are interested in the third term of order $s^{6}$, i.e. involving six masses. The aim of the above algebraic calculation was to safely derive the numerical factor in front of this contribution (which is linked to $L_{49}$). It is easy to understand which diagram this term is connected with by rewriting it as (denoting the linear-in-source solution as $\phi_{*}^{(1)}\equiv{\cal G}_{\phi}s$)

[TABLE]

where the nested brackets on each side (starting with $[2a[\cdots]\cdots]$) denote the third-order (in $s$) solution of the $\phi$ equation, i.e. the second term, $\phi_{*}^{(3)}$, in

[TABLE]

where

[TABLE]

denotes the lowest-order (quadratic in $s$) solution for $\sigma$, obtained by inserting $\phi_{*}^{(1)}$ in the effective source ($a\,\phi^{2}$) of $\sigma$ [see Eq. (42)]. The diagrammatic representation of this $O(a^{4}s^{6})$ contribution to the effective action is displayed in Fig. 1.

A useful way of reexpressing the $O(a^{4}s^{6})$ contribution to ${\cal L}_{0}^{\rm eff}[s]$ is to write it as

[TABLE]

where ${\cal L}_{\rm cubic}=a\sigma\phi^{2}$ is the cubic $\sigma$-$\phi$ coupling, Eq. (34), and where all the fields in $\delta{\cal L}_{\rm cubic}/{\delta\phi}=2a\,\sigma\phi$ on the rhs can be replaced by their lowest-order solutions.

IV Explicit computation of the contentious four-loop, 4PN effective-action contributions

We have given in the preceding section all the material needed to write down, in ${\bm{x}}$-space, all the integrals $L_{26}$ to $L_{50}$ in Fig. 1 of Foffa:2016rgu , i.e. all the $O(s^{6})$ diagrams where the $\phi$ field couples only linearly to the particles. [The other $O(s^{6})$ diagrams $L_{1}$ to $L_{25}$ in Fig. 1 of Foffa:2016rgu all involve some $\phi^{n}\cdot s$ coupling with $n\geq 2$.]

Among the integrals $L_{26}$ to $L_{50}$, we are only interested in reevaluating the three integrals ${L}_{33}$, ${L}_{49}$ and ${L}_{50}$, which contain the transcendental coefficient $\zeta(2)=\pi^{2}/6$, and whose evaluations in Ref. Foffa:2016rgu gave the problematic values (6). The method of computation used in Ref. Foffa:2016rgu was the Feynman-like one sketched above: in ${\bm{p}}$-space, with purely imaginary propagators $i{\cal K}^{-1}$, and with the use of integration by parts identities to reduce the multi-loop ${\bm{p}}$-space integrals to a subset of master integrals [one of them, ${\cal M}_{3,6}$ could only be evaluated numerically, though with such a high accuracy that they could recognize the presence of $\zeta(2)=\pi^{2}/6$ in it].

In the following three subsections we shall reevaluate the four-loop integrals ${L}_{33}$, ${L}_{49}$ and ${L}_{50}$, in ${\bm{x}}$-space, using ${\bm{x}}$-space integration by parts, and using as master integrals only the ones that have been used in our previous PN (and ADM) work, namely the original Riesz integration formula Riesz1949 [which was crucially used in the first complete computation of the (harmonic-coordinates) 2PN action (containing up to two-loop diagrams) Damour1982 ; Damour1983 ], together with the “generalized Riesz formula” (first derived in Jaranowski:1997ky for the computation of the 3PN Hamiltonian, and which was also sufficient for the computation of the local ADM 4PN Hamiltonian computation Jaranowski:2015lha ). To streamline the presentation of our computations, we will relegate most of the needed, general integration formulas to Appendix A.

IV.1 $L_{33}$

In ${\bm{x}}$-space, $L_{33}$ arises (together with its cousins $L_{26},L_{27},L_{28},L_{29},L_{30},L_{31},L_{32}$, and $L_{34}$ in Fig. 1 of Ref. Foffa:2016rgu ) from an integral of the form

[TABLE]

in which $\phi$ and $\sigma$ must be replaced by their lowest-order solutions, denoted $\phi_{*}^{(1)}$ and $\sigma_{*}^{(2)}$ above, so that $L_{\sigma\sigma\phi\phi}$ is of sixth order in the masses. The explicit expression of the integrand $\sigma\sigma\partial\phi\,\partial\phi$ is obtained from inserting Eq. (3.14) in Eq. (3.13), and reads

[TABLE]

When decomposing $\phi_{*}^{(1)}$ and $\sigma_{*}^{(2)}$ according to their mass content, i.e.

[TABLE]

and

[TABLE]

one recovers all the diagrams $L_{26}$ to $L_{27}$ (modulo some vanishing self-gravity ones, and the $1\leftrightarrow 2$ images of the previous ones).

But we are only interested in $L_{33}$, given by the spacetime diagram Fig. 2. (The thin lines represent the $\phi$ propagators, while the wavy lines represent the $\sigma$ propagators.) The $d$-dimensional projection of the diagram of $L_{33}$ is the two-point, four-loop diagram Fig. 3. (Here, the empty circles represent the two point-mass sources, i.e. the spatial projections of the thick, external worldlines in the corresponding spacetime diagram.) We see on its representation that this diagram (modulo the convention ${\cal L}_{33}^{\rm FMSS}=\frac{1}{2}L_{33}^{\rm here}$) is obtained from the general integral $L_{\sigma\sigma\phi\phi}$ by replacing each $\sigma$ by $m_{1}m_{2}\sigma_{12}$, and one $\phi$ by $m_{1}\phi_{1}$ and the other by $m_{2}\phi_{2}$, so that

[TABLE]

the factor 2 taking into account the two possibilities $\phi_{1}\phi_{2}$ vs $\phi_{2}\phi_{1}$.

We explained in Sec. III above the definitions of $\phi_{*}^{(1)}$ and $\sigma_{*}^{(2)}$ in terms of sources and propagators. In practical terms, the consideration of the Euler-Lagrange equations defined by the action (III) yields

[TABLE]

so that (using standard $d$-dimensional formulas recalled in Appendix A)

[TABLE]

where $r_{a}\equiv|{\bm{x}}-{\bm{x}}_{a}|$.

Writing the field equation for $\sigma_{ij}$ following from the action (III) yields (after a simple manipulation)

[TABLE]

so that

[TABLE]

In particular, we see that the mixed contribution $m_{1}m_{2}\sigma_{12}$ to $\sigma_{*ij}^{(2)}$ can be expressed (in ${\bm{x}}$-space) in terms of partial derivatives (with respect to ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$) of the $d$-dimensional potential $g_{d}$ defined by

[TABLE]

An explicit expression for $g_{d}({\bm{x}},{\bm{x}}_{1},{\bm{x}}_{2})$ was derived in Appendix C of Blanchet:2003gy . Let us only recall now that, when $\varepsilon=d-3\to 0$, one has the formal result

[TABLE]

so that one recovers the well-known fact (originally due to Fock Fock1955 ) that, in three dimensions,

[TABLE]

is a solution of $\Delta{\bar{g}}_{3}=r_{1}^{-1}r_{2}^{-1}$. (It will be convenient in the following to include the factor $\frac{1}{2}$ in the argument of the logarithm.)

In three dimensions, the explicit expression of $\sigma_{12}$ is

[TABLE]

where $\partial_{i}^{a}\equiv\partial/\partial x^{i}_{a}$ ($a=1,2$), while the functions $\phi_{a}$ ($a=1,2$) read

[TABLE]

It is then easily seen that, in $d=3$, the integral $L_{33}$ is convergent both in the ultraviolet (UV), i.e. near the point masses, and in the infrared (IR), i.e. at spatial infinity.

In three dimensions the integrand of $L_{33}$ [see Eq. (56) together with Eqs. (52)–(IV.1) and (64)–(65)] can be explicitly written as

[TABLE]

so that, for evaluating the integral $L_{33}$, it is enough to use the generalized Riesz formula $\hat{I}[a,b,c]$ with $c=-2$ (see Appendix A). Our final result is

[TABLE]

IV.2 $L_{49}$

In ${\bm{x}}$-space, $L_{49}$ arises (together with its cousins $L_{35},L_{36},L_{37},L_{38},L_{39},L_{41},L_{42},L_{44},L_{47}$, and $L_{48}$ in Fig. 1 of Ref. Foffa:2016rgu ) from the $O(a^{4}s^{6})$ contribution to ${\cal L}_{0}^{\rm eff}[s]$ given by Eq. (51). The corresponding spacetime diagram is displayed on Fig. 4, while its (two-point, four-loop) spatial projection is shown on Fig. 5.

Using the fact that the explicit expression of the $\sigma$-$\phi^{2}$ cubic vertex is given by Eq. (34), it is easily seen (after an integration by parts) that the latter contribution can be written as

[TABLE]

where

[TABLE]

in which $\phi$ and $\sigma$ must be replaced by their lowest-order solutions, denoted $\phi_{*}^{(1)}$ and $\sigma_{*}^{(2)}$ above. From the form of the diagram of $L_{49}$ one sees that one must keep in $\omega$ only the two pieces generated by $m_{1}m_{2}\sigma_{12}$ and bilinear in the two pieces of $\phi_{*}^{(1)}=m_{1}\phi_{1}+m_{2}\phi_{2}$. Defining

[TABLE]

where, for clariy, we put the mass labels $12$ of $\sigma$ as superscripts, we end up with

[TABLE]

where the extra factor $2$ takes into account the two orderings $\omega^{12}_{1}\omega^{12}_{2}$ vs $\omega^{12}_{2}\omega^{12}_{1}$.

The integral $L_{49}$ is IR convergent, but it has a mildly singular UV behavior because of the presence of two derivatives of $\phi_{a}$ in $\omega^{12}_{a}$ (when expanding its definition (70)). One must treat these derivatives in a distribution-theory way. After evaluating all differentiations present in $\omega_{1}^{12}$ one gets (in $d=3$)

[TABLE]

where

[TABLE]

It is not difficult to find the function $\chi_{1}$ such that $\Delta\chi_{1}=\omega_{1\,\text{fun}}^{12}$ in the sense of functions. It reads

[TABLE]

Computation of $\Delta\chi_{1}$ in the sense of distributions gives extra distibutional terms

[TABLE]

Hence, in the sense of distributions,

[TABLE]

Taking this result into account as an inverse Laplacian of $\omega_{1}$ we take

[TABLE]

Making use of Eqs. (72)–(73) and (IV.2), the integrand of $L_{49}$ can symbolically be written as

[TABLE]

where $a_{k}$, $b_{k}$, $c_{k}$, $a^{\prime}_{k}$, $b^{\prime}_{k}$, $c^{\prime}_{k}$ are integers and the coefficients $d_{k}$ and $d^{\prime}_{k}$ are rational numbers. The integral of the first part of (IV.2) is evaluated by means of the (ordinary) Riesz formula while the integral of the second part is computed by using Hadamard partie finie procedure. Our final result is

[TABLE]

IV.3 $L_{50}$

In ${\bm{x}}$-space, $L_{50}$ arises (together with its cousins $L_{40}$, $L_{43}$, $L_{45}$, and $L_{46}$ in Fig. 1 of Ref. Foffa:2016rgu ) from the effective action contribution denoted $c\sigma^{3}$ above, and defined in Eq. (III). The spacetime diagram of $L_{50}$ is displayed in Fig. 6, while its (non-planar, two-point, four-loop) spatial projection is shown in Fig. 7.

Again the term $L_{50}$ we are interested in is, as seen on its diagram, selected from this cubic expression in $\sigma$ by replacing each occurrence of $\sigma$ by its mixed piece $m_{1}m_{2}\sigma_{12}$, i.e., symbolically

[TABLE]

without any extra symmetry factor.

The integral $L_{50}$ is both IR and UV convergent. In three dimensions its integrand can be symbolically written as

[TABLE]

where $a_{k}$, $b_{k}$, $c_{k}$ are integers and the coefficients $d_{k}$ are rational numbers. For evaluation of the integral $L_{50}$ it is thus enough to use the generalized Riesz formula $\hat{I}[a,b,c]$ with $c=-3$. Our final result is

[TABLE]

IV.4 Total result, and comparison with Ref. Foffa:2016rgu

The crucial result of our new computations is that the transcendental coefficients $\sim\pi^{2}$ cancell in the sum of the three contributions $L_{33}$, $L_{49}$, and $L_{50}$:

[TABLE]

This cancellation comes about because, while our results for $L_{33}$, and $L_{49}$ agree with the corresponding results of Ref. Foffa:2016rgu recalled in Eq. (6) above, our result for $L_{50}$ differs from the corresponding result of Ref. Foffa:2016rgu by a factor $-3$:

[TABLE]

It would be interesting to understand the origin of such a missing factor $-3$ in Ref. Foffa:2016rgu . It might be caused by the presence of many $i$’s (including the ones linked to the Fourier transform of spatial derivatives $\partial_{j}\to ip_{j}$) in the quantum, ${\bm{p}}$-space calculation of $S_{\rm eff}$, together with an incorrect account of the pesky symmetry factors that enter any Wick-contraction calculation.

Anyway, we trust our result for $L_{50}$ because its normalization is very straightforwardly obtained in our ${\bm{x}}$-space computation. It would be also important to know if the error in $L_{50}^{\rm FMSS}$ has affected other integrals in Ref. Foffa:2016rgu . (Because of the cancellation of all the pole parts $\sim\frac{1}{\varepsilon}$ in the genuine $G^{5}$ contribution the cousins $L_{40}$, $L_{43}$, $L_{45}$, $L_{46}$ of $L_{50}$ cannot be uniformly affected by the same factor $-3$.)

Another reason for trusting our results is that they now reconcile the finding announced in Foffa:2016rgu that all the currently known $\pi^{2}$-dependent coefficients at order $G^{5}$ in the harmonic-ccordinates version of L_{{\text{4PN}}}^{\rm int}\big{|}^{O(G^{5})}_{{\rm linear}\ {\rm in}\ \dot{\bm{v}}_{a}} come from the double-zero reduction of the quadratic-in-acceleration terms in the original $O(G^{3})$ action, see Eq. (5). The correctness of the $O(m_{1}^{3}m_{2}^{3})$ sector of the harmonic-coordinates action of Bernard:2015njp was strongly expected in view of its agreement with the corresponding sector of the ADM action. [In terms of the $\mu$-reduced Hamiltonian, this corresponds to $O(\nu^{2})$ terms that had been unambiguously derived already in Ref. Jaranowski:2013lca .]

V Analytic computation of the master integral ${\cal M}_{3,6}$

The master integral denoted ${\cal M}_{3,6}$ in Ref. Foffa:2016rgu is the $d$-dimensional ${\bm{p}}$-space, four-loop integral depicted in Fig. 8, and defined by

[TABLE]

where ${\bar{\mathrm{d}}}k_{i}\equiv\mathrm{d}^{d}k_{i}/(2\pi)^{d}$ are normalized Fourier integrals over the loop momenta, and where the denominator is

[TABLE]

Modulo some normalization factors, this is the Fourier transform of the following $d$-dimensional ${\bm{x}}$-space integral

[TABLE]

where

[TABLE]

The diagrammatic representation (in ${\bm{x}}$-space) of the (scalar, massless) two-point, four-loop integral (87) is displayed in Fig. 9. Note that both the ${\bm{p}}$-space and ${\bm{x}}$-space representations of this Feynman integral have the shape of a four-spoked wheel.

The basic reason why the four-loop integral (87) can be analytically computed near $d=3$ by means of the generalized Riesz formula is seen in Eq. (4.11): near $d=3$, $g_{d}\equiv\Delta^{-1}(u_{1}u_{2})$ contains the (Fock) function $\ln s$, where

[TABLE]

Therefore, the integral (87) will contain (near $d=3$) a sum of terms of the type $\int\mathrm{d}^{3}x\,r_{1}^{-1}r_{2}^{-1}(\ln s)^{2}$ and $\int\mathrm{d}^{3}x\,r_{1}^{-1}r_{2}^{-1}\ln s$, which can be obtained by differentiating the generalized Riesz formula with respect to the exponent of $s$. However, there are some tricky details when implementing such a computation of (87), as will be now explained.

First, one must cope with the IR-divergence of (85), or equivalently (87). This IR-divergence is rooted in the IR-divergence of $g_{d}$ itself, which shows up in the $1/\varepsilon$ contribution (where we recall $\varepsilon\equiv d-3$) in Eq.  (4.11). Let us define

[TABLE]

and let us consider the new integral

[TABLE]

where

[TABLE]

The latter definition is such that ${\bar{g}}_{d}$ has a (point-wize) finite limit in 3 dimensions, namely

[TABLE]

We have

[TABLE]

where we defined

[TABLE]

and

[TABLE]

From Eq. (94), we see that we can reduce the computation of $I^{(d)}_{u_{1}u_{2}g^{2}}$ to that of the three integrals: $I^{(d)}_{u_{1}u_{2}{\bar{g}}^{2}}$, $I^{(d)}_{u_{1}u_{2}g}$ and $I^{(d)}_{u_{1}u_{2}}$. The last integral is trivially given by the standard $d$-dimensional777Because of the $\varepsilon$-singular factors $C_{0}\sim\frac{1}{\varepsilon}$ and $C_{0}^{2}\sim\frac{1}{\varepsilon^{2}}$ one needs to use the values of $I^{(d)}_{u_{1}u_{2}g}$ and $I^{(d)}_{u_{1}u_{2}}$ in $d$ dimensions. Riesz integral. After division by $\Omega_{d}$, where $\Omega_{d}=2\pi^{d/2}/\Gamma(d/2)$ denotes the surface of the unit sphere in $d$-dimensional Euclidean space, one finds

[TABLE]

The intermediate integral in Eq. (94), namely $I^{(d)}_{u_{1}u_{2}g}$, is a much simpler integral than $I^{(d)}_{u_{1}u_{2}{\bar{g}}^{2}}$ because it is a massless two-loop, two-point (scalar) Green function. It is depicted in Fig. 10. Such Green’s functions have been computed in the Feynman-integral literature. More precisely, the Fourier transform of $I^{(d)}_{u_{1}u_{2}g}$ (modulo some different normalization factors, including an overall sign) has been computed by Chetyrkin, Kataev, and Tkachov using Gegenbauer-polynomial, ${\bm{x}}$-space techniques888We note in passing that similar techniques have been used to compute $g_{d}$ itself in $d$ dimensions Blanchet:2003gy , and the generalized Riesz formula in 3 dimensions Jaranowski:1997ky . Chetyrkin:1980pr (see also Suzuki:2014hda ). It is trivial to compute the inverse Fourier transform of the result of Refs. Chetyrkin:1980pr ; Suzuki:2014hda (given in Appendix A), so as to compute the exact analytical expression of $I^{(d)}_{u_{1}u_{2}g}$, namely

[TABLE]

where the numerical factor (after the convenient factoring of $\Omega_{d}$, and some simplification) is found to be

[TABLE]

The $\varepsilon$-expansion of the latter numerical factor is

[TABLE]

Having the analytical expressions of $I^{(d)}_{u_{1}u_{2}g}$ and $I^{(d)}_{u_{1}u_{2}}$, the formula (94) reduces the computation of $I^{(d)}_{u_{1}u_{2}g^{2}}$ to that of $I^{(d)}_{u_{1}u_{2}{\bar{g}}^{2}}$. Though ${\bar{g}}_{d}$ has a finite limit when $d=3+\varepsilon\to 3$, and the coefficient of $I^{(d)}_{u_{1}u_{2}{\bar{g}}^{2}}$ is finite as $\varepsilon\to 0$ [so that it is enough to control $I^{(d)}_{u_{1}u_{2}{\bar{g}}^{2}}$ to $O(\varepsilon^{1})$ to get $I^{(d)}_{u_{1}u_{2}{g}^{2}}$ to $O(\varepsilon^{1})$], there are subtleties linked to the non uniformity of $\lim_{\varepsilon\to 0}{\bar{g}}_{d}$. Indeed, one must treat separately the contributions to the spatial integral $I^{(d)}_{u_{1}u_{2}{\bar{g}}^{2}}$ coming from some (large but) finite ball, say $|{\bm{x}}|<R$, and the contribution from spatial infinity, i.e. for $|{\bm{x}}|>R$. (Henceforth, we take the origin of space at the midpoint between ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$, because this significantly simplifies the asymptotic analysis at spatial infinity.) More precisely, let us write [where the factor $(1-\varepsilon)^{2}/\Omega_{d}$ is added for convenience]

[TABLE]

where

[TABLE]

and

[TABLE]

The first ($<$) integral has a limit as $\varepsilon\to 0$ which is simply given by

[TABLE]

To compute the $\varepsilon\to 0$ limit of the second ($>$) integral, we need an approximation to $(1-\varepsilon){\bar{g}}_{d}$ that is valid near spatial infinity, and in $d$ (rather than $3$) dimensions. For orientation, we recall that in $d=3$, the explicit knowledge of $\lim_{\varepsilon\to 0}{\bar{g}}_{d}\equiv{\bar{g}}_{3}=\ln\frac{s}{2}$ allows one to compute (when $r\equiv|{\bm{x}}|\to\infty$)

[TABLE]

The $d$-dimensional analog of this asymptotic expansion is obtained by combining the term-by-term inverse Laplacian of the asymptotic expansion of the source of $g_{d}$, namely

[TABLE]

with the general multipolar-expansion formula for the ($d$-dimensional) Poisson integral of an extended (but fast-decreasing at spatial infinity) source $s({\bm{x}})$:

[TABLE]

Actually, as the relevant source, $r_{1}^{2-d}\,r_{2}^{2-d}$, is not fast-decreasing when $d\approx 3$, one needs to adequately combine the two informations.999This way of combining two expansions to get the proper behavior of $d$-dimensional inverse Laplacians near spatial infinity was devised by Gerhard Schäfer and one of us (PJ) and it was never used so far in a published work. It is an IR analogue of the $d$-dimensional UV local analysis introduced in Ref. Damour:2001bu and completed (by the use of an explicit expression for the homogeneous contributions) in Appendix C4 of Ref. Jaranowski:2015lha . This leads to

[TABLE]

which is equivalent to

[TABLE]

Inserting the latter asymptotic expansion [together with the $(2-d)$th power of $r_{1}r_{2}=r^{2}(1+O(1/r^{2}))$, and $\mathrm{d}^{d}x=\Omega_{d}r^{d-1}\mathrm{d}r$] within the definition of $\frac{(1-\varepsilon)^{2}}{\Omega_{d}}I^{(d)\,>}_{u_{1}u_{2}{\bar{g}}^{2}}$ allows one to estimate the latter integral by means of a computable radial integral which yields

[TABLE]

where we introduced the function (of one variable)

[TABLE]

The appearance of the term $-f(R+\frac{r_{12}}{2})$ is exactly what is needed to define the Hadamard-regularization of the usual $3$-dimensional integral (104). Indeed, one checks that the difference

[TABLE]

has a finite limit as $R\to\infty$.

The next step is to recognize that the limit as $R\to\infty$ of (112) can be alternatively defined by an analytic continuation as $c\leadsto 0$ of the integral (over the full 3-dimensional space) of

[TABLE]

A subtle point here is that one obtains such a simple result [with the one-scale counterterm $f(R+\frac{r_{12}}{2})$] only when the exponents of $r_{1}$ and $r_{2}$ are both equal to $-1$. (Indeed, this guarantees that asymptotically $\frac{\mathrm{d}^{3}x}{4\pi}\,r_{1}^{-1}r_{2}^{-1}=\mathrm{d}r=\mathrm{d}S$, with $S\equiv r+\frac{r_{12}}{2}$.)

Let us then consider the following version of the generalized Riesz formula (with a normalization which is convenient for our present purpose)

[TABLE]

The restriction of the generalized Riesz formula to the special case $a=b=-1$ (keeping $c$ away from zero) then yields the following very simple result101010The simplicity of this result allows us to expand in powers of $c$ (i.e. to compute and integral involving integer powers of $\ln s$ by elementary means). The expansion in more general cases where $(a,b)$ deviate from $(-1,-1)$ (or other integer pairs) by $O(c)$ can also be analytically performed, though via more sophisticated techniques Moch:2001zr ; Huber:2007dx .

[TABLE]

The latter result can be easily derived from scratch by using elliptic coordinates. Indeed, in elliptic coordinates ($\xi\equiv\frac{r_{2}+r_{1}}{r_{12}}$, $\eta\equiv\frac{r_{2}-r_{1}}{r_{12}}$) one has $\mathrm{d}^{3}x/(r_{1}r_{2})=\frac{r_{12}}{2}\,\mathrm{d}\xi\,\mathrm{d}\eta\,\mathrm{d}\phi$. One then deduces that

[TABLE]

or, equivalently, in view of the previous reasonings, that

[TABLE]

Actually, the latter result can also be more directly derived simply by evaluating the $r<R$-truncated generalized Riesz integral in elliptic coordinates, which yields [for large $R$, modulo $O(1/R)$]

[TABLE]

Differentiating this result twice with respect to $c$ then yields (V).

Finally, putting together our results we can analytically compute the first three terms of the $\varepsilon$ expansion of the ${\bm{x}}$-space integral $I^{(d)}_{u_{1}u_{2}{g}^{2}}$, namely

[TABLE]

where the numerical factor (after the convenient factoring of $\Omega_{d}$) is found to be

[TABLE]

The Fourier transform (with respect to ${\bm{x}}_{1}-{\bm{x}}_{2}$) of this ${\bm{x}}$-space integral, and the addition of the various needed conventional, normalization coefficients then yields the first three terms of the $\varepsilon$ expansion of the master integral ${\cal M}_{3,6}$, namely

[TABLE]

where (with $\gamma$ denoting Euler’s constant)

[TABLE]

Our reasoning has analytically proven the latter expansion (which agrees with the result of Foffa:2016rgu ), and has, actually, reduced it to the evaluation of more elementary integrals: notably the two-loop integral $I^{(d)}_{u_{1}u_{2}g}$, and the integrals involving $\ln^{2}(s/2)$ discussed above, which were, actually, reduced to trivial integrals when using elliptic coordinates (and these trivial integrals did not involve any irrational coefficients).

Separately from the technical issue of analytically evaluating such integrals, let us note again that the evaluation of the contentious contributions to the four-loop effective action discussed in the previous sections involved only IR convergent integrals, while the master integral ${\cal M}_{3,6}$ is IR divergent (as shows up in its singular behavior as $\varepsilon\to 0$). This indicates that choosing ${\cal M}_{3,6}$ as one of the basis of elementary master integrals is probably not an optimal choice.

VI Conclusions

We have shown that remarkable cancellations take place within the four-loop, 4PN, $O(G^{5})$ “static”111111In the sense of being independent both on velocities and their time derivatives. contribution to the original, higher-time-derivative, harmonic-coordinates effective action of a gravitationally interacting binary point-mass system. Namely, the subset of diagrams $\propto G^{5}m_{1}^{3}m_{2}^{3}/(c^{8}r_{12}^{5})$ (denoted $L_{33}$, $L_{49}$, $L_{50}$ in Ref. Foffa:2016rgu ) that individually involve transcendental coefficients $\in\mathbb{Q}[\pi^{2}]=\mathbb{Q}[\zeta(2)]$ cancell against each other to leave a final, rational coefficient $+\frac{40}{3}$. On the one hand, this finding corrects a recent claim of Ref. Foffa:2016rgu , which found a final coefficient for the same terms equal to $\frac{1112}{9}-\frac{32}{3}\pi^{2}$. On the other hand, it confirms a previous lower-order finding of Foffa:2011ub , namely the fact that the corresponding highest-power-of-$G$, static terms at the previous PN level [three-loop, 3PN, $O(G^{4})$ level] did not involve any $\pi^{2}$ dependence, by contrast with the two-loop, 3PN, $O(G^{3}v^{2})$ terms. We leave to future work a deeper understanding of the rational-coefficient nature of such, highest-$G$-order, static terms at each PN order. As pointed out by Foffa et al. Foffa:2011ub ; Foffa:2016rgu , the same terms (at 3PN and 4PN) happen to be finite at $d=3$ (in dimensional regularization). At 4PN, this finiteness comes after the cancellation of poles $\propto 1/(d-3)$ present in individual diagrams. The latter cancellations can be seen rather easily, at 4PN, from the explicit ${\bm{x}}$-space expressions that we have given above for all the static 4PN diagrams (and not only $L_{33}$, $L_{49}$, $L_{50}$).

The cancellations discussed above are specific to the harmonic-gauge computation of the effective action. E.g. the situation is different in ADM gauge, where there are static, three-loop, 3PN, $O(G^{4})$ terms involving $\pi^{2}$, as well as static, four-loop, 4PN, $O(G^{5})$ terms involving $\pi^{2}$. It remains, however, true that the effective action for the gravitational interaction of point masses exhibit a remarkably small level of transcendentality. At one and two loops (at 1PN and 2PN), the action involves only rational coefficients. The 3PN, three-loop level introduces $\mathbb{Q}[\zeta(2)]$ coefficients, and this transcendentality level does not increase when going to the 3PN, four-loop level. Very-high-PN-order, analytical gravitational self-force studies of the EOB Hamiltonian Bini:2014nfa ; Bini:2015bla ; Kavanagh:2015lva have shown (for a subset of the diagrams) that the transcendentality level increases only quite slowly as the loop number (equal to the PN level) increases: the $\mathbb{Q}[\zeta(4)]=\mathbb{Q}[\pi^{4}]$ level is reached at six loops, and $\zeta(3)$ first appears at the seven-loop order. [Here, we are (roughly) subtracting the effects linked to non-local-in-time interactions which introduce Euler’s constant $\gamma$ and logarithms.] We leave to future work a better understanding of such facts.

Separately from the interest of finding special structures hidden in the gravitational effective action, our work provides a confirmation of the correctness of the 4PN-level $O(m_{1}^{3}m_{2}^{3})$ sector of the harmonic-coordinates action of Bernard:2015njp . This confirmation is independent of that following from its previously checked agreement with the corresponding sector of the 4PN, ADM action of Damour:2014jta ; Jaranowski:2015lha . [In terms of the $\mu$-reduced Hamiltonian, this corresponds to $O(\nu^{2})$ terms that had been first derived in Ref. Jaranowski:2013lca .] Having such independent confirmations is always useful. It would be useful that a full, independent 4PN, EFT-based computation of the 4PN effective action be performed. However, in view of the complications (and sign dangers) brought by working with purely imaginary propagators, and corresponding $i$-decorated vertices, we would advocate (as explained at the end of Sec. II above) to work with real propagators ${\cal G}=-{\cal K}^{-1}$, and corresponding $i$-free vertices (when viewed in ${\bm{x}}$-space).

Let us finally comment on the technicalities of the explicit, 4PN computation. We have shown in Sec. V above that the four-loop master integral ${\cal M}_{3.6}$ selected as basis element in Foffa:2016rgu , and that could only be numerically computed in the latter reference, could be analytically computed by means of what has been the standard tool in ADM computations since the 3PN level, namely the generalized Riesz formula Jaranowski:1997ky . It is remarkable that a tool set up for the three-loop level can analytically deal with a four-loop integral that resisted the state-of-the-art technologies in multi-loop computations. We think that this is due to two main facts (besides the special structure of the gravitational vertices): (i) our use of ${\bm{x}}$-space integration121212We note in this respect that the first five four-loop master integrals in Fig. 3 of Foffa:2016rgu are trivial to compute in ${\bm{x}}$-space. Indeed, repeated lines between two points just mean a power of $r_{12}$, $r_{1}$ or $r_{2}$, without any needed integration. The only needed integrations in ${\bm{x}}$-space go with the number of intermediate vertices; for instance, ${\cal M}_{1,1}$, ${\cal M}_{1,2}$, ${\cal M}_{1,3}$ and ${\cal M}_{1,4}$ only involve one intermediate-point integration, so that they are immediately derived from the normal Riesz formula, while ${\cal M}_{3,6}$ involves three intermediate-points integrations. and (ii) the fact that the repeated differentiation of the generalized Riesz formula with respect to the power of $s\equiv r_{1}+r_{2}+r_{12}$ allows one to compute integrals that can show up at an arbitrary high loop order. Indeed, the $n$th derivative with respect to the power of $s$ generates (3-dimensional) integrals of the type, say

[TABLE]

When $n_{1}=n_{2}=1$ and $n=2$ this corresponds to the four-loop ${\cal M}_{3.6}$ diagram. When taking higher values of $n$, $I_{n_{1},n_{2},n}$ describes higher-loop master integrals.

Acknowledgements.

T.D. thank Pierre Vanhove for informative discussions, and useful references, on Feynman integrals. We thank Stefano Foffa for clarifying the precise meaning of the notation for the kinetic terms in the bulk action for $\sigma_{ij}$ and $\phi$, and notably $(\overrightarrow{\nabla}\sigma)^{2}$. The work of P.J. was supported in part by the Polish NCN Grant No. UMO-2014/14/M/ST9/00707.

Appendix A Some useful formulas

A.1 $d$-dimensional results

The area of the $(d-1)$-dimensional unit sphere in $\mathbb{R}^{d}$ reads

[TABLE]

It is convenient to introduce the constant

[TABLE]

such that

[TABLE]

Then the $d$-dimensional Newtonian potential $u_{a}\equiv r_{a}^{2-d}$ fulfills the equation

[TABLE]

The (ordinary) Riesz formula in $d$ dimensions reads

[TABLE]

A three-dimensional generalization of the Riesz formula (128) for integrands of the form $r_{1}^{a}r_{2}^{b}(r_{1}+r_{2}+r_{12})^{c}$ was derived in Ref. Jaranowski:1997ky . It reads

[TABLE]

The function $I_{1/2}$ is defined as follows:

[TABLE]

where $B$ is the Euler beta function and $B_{1/2}$ is the incomplete beta function which can be expressed in terms of the Gauss hypergeometric function ${}_{2}F_{1}$:

[TABLE]

The $d$-dimensional Fourier transform of a power reads:

[TABLE]

The result of Chetyrkin:1980pr (and Suzuki:2014hda ) for the ${\bm{p}}$-space version of the two-loop diagram of Fig. 10 reads

[TABLE]

where

[TABLE]

A.2 Distributions in $d=3$ dimensions

In Sec. IV we have to compute different distributional derivatives. We collect here formulae which can be used for this goal. Let us start from identities involving Dirac delta distrubution and its derivatives:

[TABLE]

Because usually the function $f$ for which the above identities are used is singular at $\bm{x}=\bm{x}_{a}$, the symbol $f_{\text{reg}}(\bm{x}_{a})$ means the regularized “partie finie” value of the function $f$ at $\bm{x}=\bm{x}_{a}$ (for its definition and properties see, e.g., Appendix A4 of Ref. Jaranowski:2015lha ).

We have also employed distributional derivatives to calculate first and second partial derivatives of homogeneous functions $1/r_{a}$, $1/r_{a}^{2}$, and $1/r_{a}^{3}$ (for derivation and properties see, e.g., Appendix A5 of Ref. Jaranowski:2015lha ). The first partial derivatives read

[TABLE]

The second partial derivatives are

[TABLE]

Tracing the above formulas yields

[TABLE]

As an application of the above formulas let us note a useful expression which shows how to compute the Laplacian of the product of a (singular at $\bm{x}=\bm{x}_{a}$) function $f$ and $1/r_{a}^{3}$:

[TABLE]

where $\Delta\left(f/r_{a}^{3})\right)|_{\textrm{ord}}$ means the Laplacian computed using standard (i.e. non-distributional) rules of differentiations.