Emergence of Calabi–Yau manifolds in high-precision black-hole scattering

Integrand generation and integral family

We use the WQFT formalism^11,57,63 that quantizes the worldline deflections \(z_i^\mu (\tau )\) and graviton field h_μν(x) arising in the background field expansions \(x_i^\mu (\tau )=b_i^\mu +v_i^\mu \tau +z_i^\mu (\tau )\) and \(g_\mu \nu =\eta _\mu \nu +\sqrt32\rm\pi Gh_\mu \nu \), respectively (now setting c = 1). In the gravitational sector, we use a nonlinearly extended de Donder gauge that simplifies the three-graviton and four-graviton vertices (see Supplementary Information). The worldline actions (equation (1)) are improved by making use of the proper time gauge \(\dotx_i^2=1\) for the ith black hole:

This ensures a linear coupling to the graviton h_μν. At the present four-loop (G⁵ or 5PM) order, we require up to six-graviton vertices that derive from the Einstein–Hilbert action plus gauge-fixing term—taken in D = 4 − 2ϵ dimensions. We also require the single-graviton emission plus (0,…, 5)-deflection vertices derived from equation (8). We provide the explicit vertices and graviton gauge-fixing function in a Zenodo repository submission⁶⁴ accompanying this article; an analytic expression for the n-deflection worldline vertex was given in refs. ^11,63.

The full 5PM integrand is generated using the Berends–Giele-type recursion relation discussed in ref. ⁶⁵ and sorted into five self-force (SF) sectors according to their scaling with the masses m₁ and m₂:

The powers of the masses follow from the number of times a worldline is touched in a given graph. Here we compute the sub-leading self-force (1SF) contributions \(\Delta p_\rm1SF^(5)\mu \) and \(\Delta p_\overline1\rmS\rmF^(5)\mu \), as well as reproducing the 0SF contributions \(\Delta p_\rm0SF^(5)\mu \) and \(\Delta p_\overline0\rmS\rmF^(5)\mu \) that follow from the geodesic motion in a Schwarzschild background. The resulting integrand is reduced to scalar integrals by means of tensor reduction and ‘planarized’ using partial fraction (eikonal) identities as described in ref. ⁵³. In summary, all integrals are mapped to the 5PM-1SF planar family \(\int _\ell :=\int d^D\ell /(2\pi )^D\), \(\bar\delta (x)\,:=\,2\pi \delta (x)\)

in which σ and n denote causal i0⁺ prescriptions and integer powers of propagators, respectively. The four worldline propagators D_i(σ_i) appearing are (i = 2, 3, 4)

and the 14 gravitons propagators D_IJ with I = (0, 1, i, q) are

There are at most three bulk graviton propagators D_1i that may go on-shell at 5PM order.

IBP reduction

IBP identities^66,67,68 are used to reduce to master integrals. We use a future release of KIRA 3.0 (refs. ^69,70) adapted to our needs that uses the FireFly^71,72 library for reconstructing rational functions through finite-field sampling. We have 45 top-level sectors in the 5PM-1SF family that have been described in ref. ⁵³. The integrals encountered in the planar family (equations (10a)–(10c)) have propagator powers in the range n_i/IJ ∈ [−9, 8]. The strategies applied to reduce the runtime of the IBP reduction are comparable with the conservative case⁵³. The final set of needed IBP replacement rules to master integrals generated comprises about 30 GB of data and can be made available on request.

Differential equations and function space

The method of differential equations^73,74,75 is used, in which the matrices in equation (3) depend on the parameters x and the dimensional regulator ϵ = (4 − D)/2. We take the physical limit ϵ → 0 to compute our observables. Therefore, we need the solutions of the integrals expanded in ϵ. To systematically compute this expansion, we transform equation (3) into canonical form⁷⁵ such that the ϵ dependence is factored out of the differential equation matrix:

with J(x, ϵ) = \(\widehat\bfT(x,\epsilon )\)I(x, ϵ) and \(\epsilon \widehatA(x)=(\widehatT(x,\epsilon )\widehatM(x,\epsilon )+d\widehatT(x,\epsilon )/dx)T(x,\epsilon )^-1\). The solution is then a path-ordered matrix exponential:

in which j encodes the boundary values of our integrals at x = 1.

We take a bottom-up approach to determine the required transformation \(\widehatT(x,\epsilon )\). For this, we sort our integrals into groups sharing the same set of propagators. These so-called sectors are ordered from lower (fewer propagators) to higher (more propagators), resulting in the block diagonal matrix in Fig. 3. We begin by ϵ-factorizing the lower sectors and then move on to the higher sectors. First, we transform the diagonal blocks, which are identified with the maximal cuts⁷⁶, into ϵ-form and then proceed to the off-diagonal contributions. The IBP reductions were fully completed for all integrals before the diagonal blocks were identified. Particularly for handling sectors coupled to the CY3 diagonal sector, it is important to choose a good initial basis of integrals such that the relevant couplings are as simple as possible. As we proceed to canonicalize, we adapt and improve our choice of initial basis accordingly.

The simplest diagonal blocks to canonicalize are those containing only multiple polylogarithms^77,78,79, depicted in lilac in Fig. 3. The algorithm CANONICA⁸⁰ finds the necessary transformation to ϵ-form for these blocks by making a suitable ansatz. It is noteworthy that the complexity of this transformation, as well as the runtime, depends highly on the choice of initial integrals. In general, we pick our initial basis integrals so that the regulator ϵ does not appear non-trivially in the denominators of the differential equation (3).

Diagonal blocks containing a K3 surface, depicted in purple in Fig. 3, are handled using INITIAL⁸¹. The INITIAL algorithm requires a pure seed integral to construct an ϵ-factorized differential equation through an ansatz tailored to the specific seed integral. Note that, in this context, a pure integral is given as a linear combination of iterated integrals (equation (4)) having no non-trivial pre-factors in x. For non-pure integrals, these pre-factors are non-trivial functions and are known as leading singularities. For each K3 sector, we find an appropriate seed integral by analysing the diagonal block at ϵ = 0. We decouple each diagonal K3 block by switching to a derivative basis: \(\bfI_\rmK3=(I_1,I_2,I_3)\to (I_1,I_1^\prime ,I_1^\prime\prime )\) or \(\bfI_\rmK3=(I_1,I_2,I_3,I_4)\to (I_1,I_1^\prime ,I_1^\prime\prime ,I_4)\), depending on the size of the block, in which I_i are the master integrals of this sector—I₄ is chosen so that it decouples from I₁ and its derivatives as much as possible. The choice of I₁ ensures that its third-order differential (Picard–Fuchs) equation (\(\widehat\theta =x\frac\rmd\rmdx\)),

has the explicit solution \( _\epsilon =0\propto \varpi _\rmK3=\left(\frac2\pi \right)^2K^2(1-x^2)\), that is, it is proportional to a K3 period. Unlike the polylogarithmic case, this third-order differential equation is not factorizable into first-order equations. Using the normalized integral I₁/ϖ_K3 as the seed, INITIAL may then construct an ϵ-form for the diagonal parts of our K3 sectors (similar to the 4PM case).

To canonicalize the single diagonal CY3 block, depicted in red in Fig. 3, we follow the discussion in ref. ⁵². Similar to K3, we first pick a suitable starting integral and make a basis change \(\bfI_\rmCY3=(I_1,I_2,I_3,I_4)\to (I_1,I_1^\prime ,I_1^\prime\prime ,I_1^\prime\prime\prime )\). The starting integral I₁ is chosen so that, when ϵ = 0, it satisfies the Picard–Fuchs equation (5), which is a hypergeometric system. This implies that the periods of our CY3 geometry are given in terms of hypergeometric functions, for example, \(x_4F_3\left[\frac12,\frac12,\frac12,\frac12;1,1,1;x^4\right]\). Moreover, this also gives rise to intriguing arithmetic properties discussed in ref. ⁵⁹. From this equation and its four fundamental solutions ϖ₀,…, ϖ₃, which we collect in the Wronskian matrix W(x) = (∂^jϖ_i) with 0 ≤ i, j ≤ 3, we construct in several steps the rotation matrix into ϵ-form. This approach was invented in ref. ⁸² and further developed in ref. ⁸³, and, in the K3 case, it is equivalent to the INITIAL algorithm. The process involves the following three steps:

1. 1.

   We split the Wronskian matrix into a semi-simple and unipotent part \(\widehatW=\widehatW^\rmss\times \widehatW^\rmu\). Naively, we can understand this splitting as a decomposition of the maximal cuts of the CY3 sector into their leading singularities and pure parts. The unipotent part is named after the unipotent differential equation it fulfils:

   $$(\rmd-A^\rmu(x))\widehatW^\rmu(x)=0,$$

   (14)

   in which A^u(x) is nilpotent. The matrix A^u(x) can generally be written in terms of invariants, known as Y-invariants, of the CY variety and was explicitly given for our CY3 in ref. ⁵².

2. 2.

   We rotate our integrals with the inverse of W^ss, which strips them of their leading singularities. (The analogous operation for K3 is normalizing I₁ with the K3 period ϖ_K3, leaving only the pure part.) To parametrize all degrees of freedom in W^ss of a CY3, we need the holomorphic solution ϖ₀, an extra function

   $$\alpha _1=\frac\varpi _^2x(\varpi _\varpi _1^\prime -\varpi _^\prime \varpi _1),$$

   (15)

   called the structure series, and their derivatives. The appearance of α₁ is new compared with the K3 case and shows the increased complexity in structure of a CY3. In ref. ⁸⁴, the structure series α₁ (and more generally the Y-invariants) are generally defined and used to construct a normal form of a CY differential equation; more specifically, for our CY3, they were derived in ref. ⁵². For ϵ-factorizing our differential equations, it is important that this normal form of the CY differential equation is in a factorized form with respect to its derivatives. To eliminate all redundancies in this step, we must use Griffiths transversality—an essential property of CY geometries that yields quadratic relations between their periods—to simplify the form of \(\widehatW^\rmss\).

3. 3.

   After completing these steps and appropriately rescaling in ϵ to arrange the weights of the integrals, the diagonal block of our CY3 looks like:

   $$\beginarrayl\frac\rmd\rmdx\widetilde\bfI_\rmCY3\,=\,\mathop\sum \limits_i=-2^1\epsilon ^i\widehat\bfM_\rmCY3^i(x)\widetilde\bfI_\rmCY3,\,\rmwith\\ \widetilde\bfI_\rmCY3\,=\,T_\epsilon -\rmscalings(\widehatW^\rmss)^-1\bfI_\rmCY3.\endarray$$

   (16)

   We find an ϵ-form by acting with a suitable set of transformation matrices on \(\widetilde\bfI_\rmCY3\), working order by order in ϵ starting from ϵ⁻². The process requires us to introduce four new functions G_i(x) (i = 1,…, 4), which obey a first-order differential equation containing ϖ₀, α₁, their derivatives and G_j(x) functions with j < i. For example,

   $$G_1^\prime (x)=-\frac96x(x^4+1)\varpi _(x)^2(x-1)^2(x+1)^2(x^2+1^2)\alpha _1(x).$$

   (17)

   Because of this structure, the functions G_i(x) are all expressible as iterated integrals of CY periods and associated functions. These functions were previously introduced in terms of a different variable in ref. ⁵² and are listed for our conventions in the Supplementary Information of this article.

We now have the ϵ-form of the diagonal part of the CY3 sector and, thus, a canonical form of all diagonal blocks. We refer to this intermediate basis, in which all diagonal blocks are in ϵ-form, as \(\mathfrakJ\). The next stage in canonicalization involves tackling the off-diagonal blocks. To do so, we distinguish between off-diagonal blocks coupled to the CY3 sector, which require special care, and the rest.

We have developed our own algorithm to transform the off-diagonal entries of our differential equation that do not couple to the CY3 sector but can depend on K3 functions. This algorithm uses FINITEFLOW⁸⁵ and MultivariateApart⁸⁶ and provides suitable ansätze also including elliptic contributions for the required transformations. It is similar to algorithms used for polylogarithmic off-diagonals, such as those found in CANONICA or Libra⁸⁷.

For sectors polylogarithmic or K3 on their diagonal blocks, yet also coupled to the CY3 sector, a good initial basis of integrals is essential to minimize these couplings. One possibility to identify such candidates is to perform an integrand analysis, usually done in the Baikov representation^88,89 of the integrals. However, we instead found it simpler to use the diagonals themselves to derive constraints on the initial choice of integrals, expanding on the ideas of ref. ⁹⁰. Having now found canonical masters on the diagonals, that is, the maximal cuts, our strategy is to choose initial candidate integrals that are related as closely as possible to these canonical masters within their respective diagonal blocks. More precisely, we search for candidates that, on their diagonal blocks, are given by a linear combination of the canonical maximal cut integrals and overall functions of ϵ and x:

in which the c_k are constant numbers.

We expect that such a ‘good’ choice of candidate integrals only requires minimal corrections to form a canonical basis. For certain sectors that are polylogarithmic on their corresponding diagonal block, we need to enlarge these types of constraint by combining different polylogarithmic sectors and requiring equation (18) to hold beyond a single diagonal block. Thus, we obtain further conditions resulting from off-diagonal couplings between separate polylogarithmic blocks. In some instances, we can also relax the condition in equation (18) by considering the c_k(x) as functions of x and still find easy transformations. The use of IBPs makes this procedure efficient and allows us to find all transformations for the coupling to the CY3 manually, proceeding similarly as for the diagonal of the CY3 sector. We build successive transformations, removing iteratively all nonlinear-in-ϵ contributions, starting with the highest negative power of ϵ. By doing so, for some integrals, we need to introduce 16 new functions G₅(x),…, G₂₀(x), which again satisfy first-order differential equations listed in the Supplementary Information. More specifically, for the mixings between K3 and CY3 sectors, we introduce new functions whose first-order differential equations contain the periods of both the K3 and the CY3. For example,

This concludes the canonicalization process of the whole differential equation system.

Having converted our matrix into its canonical form, \(\widehatA(x)\) provides all integration kernels needed to express the master integrals as iterated integrals (equation (4)). We selected a set of linearly independent kernels by examining their small velocity expansion. Our observables consist of iterated integrals that include K3, CY3 and mixed integration kernels, functions from the rotation matrix \(\widehatT(x,\epsilon )\) and algebraic functions from the decomposition in terms of initial master integrals. We need at most four-times-iterated integrals; all multiple polylogarithms are constructed from the kernels \(\\frac1x,\frac11\pm x,\fracx1+x^2\\) and have a maximum transcendental weight of 3. Let us also note that the K3 and CY3 periods occurring above are related to those of the Legendre curve by a symmetric and a Hadamard product, respectively^52,59.

Boundary fixing

A complete solution to the differential equation (3) requires the determination of integration constants in the form of boundary integrals, that is, master integrals in the static limit x → 1 (v → 0), which are functions only of ϵ. As the integration and x → 1 limit do not commute, we use the method of regions^91,92,93 to isolate contributions with definite (ϵ-dependent) scalings in the velocity v and series-expand integrals at the level of the integrand. These so-called regions are associated with different velocity scalings of the bulk graviton momenta ℓ_i, which can be either potential (P) or radiative (R):

There are three propagators D₁₂, D₁₃, D₁₄ that may enter both regions; the rest are kinematically restricted to P (including the velocity-suppressed P: (v², v)) by the presence of energy-conserving delta functions δ(ℓ_i · v_j). We thus denote the four possible regions as (PPP), (PPR), (PRR) and (RRR). We needed to evaluate 14 + 14 (even + odd) boundary integrals in the (PPR), 5 + 5 in the (PRR) and 4 + 4 in the (RRR) sectors, as well as the 28 + 18 (PPP) boundary integrals that were already determined in the conservative case⁵³ (here there is no distinction between Feynman and retarded bulk propagators). We perform all integrals analytically and check them numerically using pySecDec^{94,95,96,97,98,99} and AMPred^{100,101,102,103}.

All new boundary integrals for us to evaluate, as compared with ref. ⁵³, contain radiative graviton propagators. Our main strategy is to perform them by means of Schwinger parametrization, in which we also benefit by explicitly integrating out loops involving only gravitons, leaving lower-loop integrals. In doing so, we evaluate at most two-loop integrals in all but two cases. These two cases are genuine three-loop integrals for which we did not manage the integration over Schwinger parameters and therefore moved over to the time domain to establish identities at the level of the series coefficients. In general, our integrations yield generalized hypergeometric functions _pF_q, which can be series-expanded in ϵ by numerically expanding them to high precision and reconstructing analytic expressions using an ansatz and an integer relation algorithm. Expressions for all boundary integrals, and our methodology for deriving them, are elaborated in the Supplementary Information.

Results

Our main results are full expressions—including dissipation—for the 1SF and \(\overline1\rmS\rmF\) parts of the momentum impulse, \(\Delta p_\rm1SF^(5)\mu \) and \(\Delta p_\overline1\rmS\rmF^(5)\mu \) (equation (9)). They are expanded on basis vectors:

also pulling out an overall factor of the impact parameter b. This decomposition constitutes a split into parts originating from integrals of even and odd parity in \(v_i^\mu \). The basis vectors are the impact parameter \(\widehatb^\mu =(b_2^\mu -b_1^\mu )/b\) and dual velocity vectors \(\checkv_1^\mu =(\gamma v_2^\mu -v_1^\mu )/(\gamma ^2-1)\) and \(\checkv_2^\mu =(\gamma v_1^\mu -v_2^\mu )/(\gamma ^2-1)\). The main coefficients of interest here are \(\c_b,\barc_b\\) and \(\c_v,\barc_v\\), because the set \(\c_v^\prime ,\barc_v^\prime \\) being determined by lower-PM results from using preservation of mass \(p_1^2=(p_1+\Delta p_1)^2\). The coefficients are further decomposed into those with an even or odd number of radiative gravitons in the boundary fixing:

with w ∈ b, v. The even part is defined from (PPP) + (PRR) and the odd part from (PPR) + (RRR). The two parts have distinctive parities under flipping the sign of the relative velocity v → −v. Even and odd sectors then give rise to integer and half-integer post-Newtonian (PN) orders, respectively.

The 1SF and \(\overline1\rmS\rmF\) parts of the impulse therefore each consist of four non-trivial coefficients, labelled by b or v, each having an even or odd number of radiative gravitons R. We expand each of these in sets of basis functions F(γ) with coefficient functions d(γ), being polynomial up to factors of (γ² − 1),

barred coefficients being expanded in the same way. Here z ∈ even, odd counts the parity of radiative gravitons. The second line with a logarithm produced from the cancellation of 1/ϵ poles between the different boundary regions is associated with tails. It is relevant for all coefficients except c_b,odd, \(\barc_b,\rmodd\) and c_v,even. The b-type basis functions are all multiple polylogarithms of maximum weight three and, thus, relatively simple. By contrast, the v-type basis functions are much more complex. They generally have the structure of equation (4) with kernels depending on the CY periods. All basis functions and polynomial coefficients are provided in the ancillary file in our repository submission⁶⁴.

The total loss of four-momentum at 5PM order, using equation (7), is given schematically by (see, for example, ref. ¹⁰⁴):

Our \(1\rmSF/\overline1\rmS\rmF\) result fixes all coefficients except r₄(γ). Similarly, we may derive the relative scattering angle^105,106 \(\theta =\arccos (\bfp_\rmin\cdot \bfp_\rmout/| \bfp_\rmin| | \bfp_\rmout| )\). Here p_in = p₁ = −p₂ is the incoming momentum in the centre-of-mass frame and

is the (relative) outgoing momentum. We expand the scattering angle in G:

Both the angle and \(P_\rmrad^(5)\mu \) expansion coefficients are separated into parts even and odd under v → −v and expanded on suitable basis functions. All of these results can be found in the observables.m ancillary file in our repository submission⁶⁴.

Checks

Our result for the impulse has been checked internally by the cancellation of 1/ϵ poles and obeying the mass condition \(p_1^2=(p_1+\Delta p_1)^2\)—verifying the \(\c_v^\prime ,\barc_v^\prime \\) coefficients (equations (21a) and (21b)). We have also checked that the (PPP) 1SF contribution is related to its \(\overline1\rmS\rmF\) counterpart by symmetry. Furthermore, we have performed several checks using the non-relativistic (v → 0) limit. First, PN results for the relative scattering angle at 5PM and to first order in self-force are known to 5.5PN order^104,105,107. These include conservative terms at integer PN orders, starting from 0PN order, which were already matched in ref. ⁵³, plus dissipative terms appearing at 2.5PN, 3.5PN, 4PN, 4.5PN, 5PN and 5.5PN orders, which we reproduce. Furthermore, and in the same works^104,105,107, dissipative PN results for the radiation of energy and recoil of the two-body system were reported at a relative 3PN order to their leading order. These allow us to perform non-trivial checks on r₁(γ), r₂(γ) and r₃(γ) of equation (24) to relative 3PN order.

As an example of our results, we print explicitly the series expansion in v of the G⁵ component of E_rad:

The terms in the square brackets up to v⁶ reproduce the known PN results of E_rad, with the remaining three lines providing further hitherto unknown terms. Naturally, this series expansion can be extended to any order in v using our present results. Explicit results relevant for the PN checks of the relative scattering angle and recoil are given in the Supplementary Information.