# Vanishing of quadratic Love numbers of Schwarzschild black holes

Simon Iteanu, Massimiliano Maria Riva, Luca Santoni, Nikola Savić, Filippo Vernizzi

► To cite this version:

Simon Iteanu, Massimiliano Maria Riva, Luca Santoni, Nikola Savić, Filippo Vernizzi. Vanishing of quadratic Love numbers of Schwarzschild black holes. Journal of High Energy Physics, 2025, 10.1007/JHEP02(2025)174. hal-04970788

HAL Id: hal-04970788

https://hal.science/hal-04970788v1

Submitted on 28 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Vanishing of quadratic Love numbers of Schwarzschild black holes

Simon Iteanu<sup>D, a,b</sup> Massimiliano Maria Riva<sup>D, c,d</sup> Luca Santoni<sup>D, b</sup> Nikola Savić<sup>D, e</sup>
and Filippo Vernizzi<sup>D, e</sup>

E-mail: simon.iteanu@ijclab.in2p3.fr, mr4533@columbia.edu,
santoni@apc.in2p3.fr, nikola.savic@ipht.fr, filippo.vernizzi@ipht.fr

**ABSTRACT:** The induced conservative tidal response of self-gravitating objects in general relativity is parametrized in terms of a set of coefficients, which are commonly referred to as Love numbers. For asymptotically-flat black holes in four spacetime dimensions, the Love numbers are famously zero in the static regime. In this work, we show that this result continues to hold upon inclusion of nonlinearities in the theory for Schwarzschild black holes. We first solve the quadratic Einstein equations in the static limit to all orders in the multipolar expansion, including both even and odd perturbations. We show that the second-order solutions take simple analytic expressions, generically expressible in the form of finite polynomials. We then define the quadratic Love numbers at the level of the point-particle effective field theory. By performing the matching with the full solution in general relativity, we show that quadratic Love number coefficients are zero to all orders in the derivative expansion, like the linear ones.

**KEYWORDS:** Black Holes, Classical Theories of Gravity, Effective Field Theories

ARXIV ePRINT: 2410.03542

<footnote><sup>a</sup>Ecole Normale Supérieure de Lyon, Laboratoire de Physique,
Lyon F-69342, France</footnote>

<footnote><sup>b</sup>Université Paris Cité, CNRS, Astroparticule et Cosmologie,
10 Rue Alice Domon et Léonie Duquet, F-75013 Paris, France</footnote>

<footnote><sup>c</sup>Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany</footnote>

<footnote><sup>d</sup>Department of Physics, Center for Theoretical Physics, Columbia University,
New York, 538 West 120th Street, NY 10027, U.S.A.</footnote>

<footnote><sup>e</sup>Université Paris-Saclay, CNRS, CEA, Institut de Physique Théorique,
91191 Gif-sur-Yvette, France</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Static quadratic perturbations</td><td>4</td></tr><tr><td>2.1</td><td>Metric parametrization and gauge choice</td><td>5</td></tr><tr><td>2.2</td><td>Linear perturbations</td><td>6</td></tr><tr><td>2.3</td><td>Quadratic perturbations</td><td>9</td></tr><tr><td>2.4</td><td>Selection rules</td><td>12</td></tr><tr><td>3</td><td>Solutions</td><td>13</td></tr><tr><td>3.1</td><td>Linear solutions</td><td>13</td></tr><tr><td>3.2</td><td>Quadratic solutions</td><td>15</td></tr><tr><td>3.3</td><td>Quadratic solutions induced by quadrupole × quadrupole</td><td>17</td></tr><tr><td>4</td><td>Worldline effective field theory</td><td>19</td></tr><tr><td>4.1</td><td>Point-particle action and nonlinear Love numbers</td><td>20</td></tr><tr><td>4.2</td><td>Background field method</td><td>22</td></tr><tr><td>5</td><td>Matching with the point particle</td><td>25</td></tr><tr><td>5.1</td><td>The external field</td><td>25</td></tr><tr><td>5.2</td><td>Quadratic response</td><td>28</td></tr><tr><td>5.3</td><td>Matching</td><td>31</td></tr><tr><td>6</td><td>Higher multipoles</td><td>32</td></tr><tr><td>6.1</td><td>General polynomial form of second-order static solutions</td><td>33</td></tr><tr><td>6.2</td><td>Matching to all orders in derivatives</td><td>37</td></tr><tr><td>7</td><td>Conclusions</td><td>39</td></tr><tr><td>A</td><td>Integration of tensor harmonics</td><td>40</td></tr><tr><td>B</td><td>Coefficients</td><td>42</td></tr><tr><td>C</td><td>From spherical to Cartesian coordinates</td><td>44</td></tr><tr><td>D</td><td>Symmetric trace-free tensor relations</td><td>45</td></tr><tr><td>E</td><td>Recovering the Schwarzschild metric</td><td>46</td></tr></table>

## 1 Introduction

Gravitational waves from coalescing binary systems represent a unique channel to study the fundamental properties of gravity and compact objects in the strong-field regime. As the number of merger events observed by the LIGO-Virgo-KAGRA network [1] increases and, thanks also to future facilities and envisioned new detectors, we move toward an era of precision physics with gravitational waves, ever more accurate waveform templates
will be necessary [2–7]. This requires in particular to have a precise understanding of the
conservative and dissipative dynamics of two-body systems (see e.g., [8–10]), which includes
tidal effects [11–13].

The tidal deformability of a compact object is described in terms of a set of coefficients, which capture the conservative and dissipative induced response when the object is acted upon by an external long-wavelength tidal field [14–16]. The coefficients associated with the conservative response are often called Love numbers and, together with the dissipative numbers, they carry relevant information about the object’s structure and interior dynamics. This makes the tidal response coefficients an important target for current and future gravitational-wave observations from merging binary systems, which can offer valuable insights into the equation of state of neutron stars [11, 17–19], the physics at the horizon of black holes [14–16, 20–28], as well as the existence of more exotic compact objects [29–34].

In the Post-Newtonian (PN) regime of the inspiral of two compact objects, the leading-order tidal effects are associated with tidal heating, and start affecting the gravitational waveform at 2.5PN order for spinning objects and 4PN order for non-rotating ones [8, 35, 36]. Conservative tidal effects are instead estimated to become relevant from 5PN order, in the limit of static tides (see, e.g., [37, 38]). Although these are the most prominent tidal effects that are expected at the level of the gravitational waveform, there are in general very good reasons to study *subleading* corrections to the physics of compact sources beyond the point-particle approximation. First of all, the calculation of subleading effects, in addition to putting leading-order results on solid grounds, may be necessary to probe and get insight into interesting properties of the system that do not enter at leading order. Moreover, subleading effects might be useful to recalibrate effective models, improve waveform and break degeneracy among observables. They are especially relevant whenever leading-order effects happen to be vanishing, making the next-to-leading-order corrections the actual dominant contributions. One notable example is given by black holes in General Relativity (GR), for which the static Love numbers — associated with the leading-order, conservative, finite-size terms in an expansion in the number of time derivatives — are absent for all multipoles [14–16, 20–22, 24–28]. This fact makes the study of subleading corrections to the conservative response of black holes particularly interesting at both theoretical and observational level.

One main motivation, at the theoretical level, to compute subleading tidal effects of
black holes has to do with symmetries. The vanishing of the static Love numbers in GR
has for long time been recognized as an outstanding naturalness puzzle in gravity from
an effective-field-theory perspective [39, 40]. Recently, symmetry explanations have been
proposed as possible resolutions to the puzzle [23, 24, 41–45]. It is therefore particularly
interesting to understand what is the fate of these symmetries, whether they are an ‘accident’
of the static regime of the linearized Einstein equations for the perturbations, or whether
they instead persist beyond leading order.

At subleading order, corrections to the tidal Love numbers of black holes are of two types: they can either result from finite-frequency effects, or from gravitational nonlinearities.<fnref n="1" /> The former, which are often referred to as *dynamical* Love numbers, capture the conservative

<footnote n="1">¹In the language of effective field theories, the former correspond to higher-order time-derivative operators, while the latter correspond to operators with higher number of fields. See section 4 for details.</footnote>
response induced by the presence of a time-dependent tidal field; in the case of black holes, they have been studied in e.g. [27, 46–49]. The latter, which will be the focus of the present work, stem instead from the field nonlinearities in the Einstein equations. Similarly to electromagnetism, where nonlinear polarization theory can be used to describe nonlinear optical effects, such as the nonlinear polarization of an optical medium [50], one can formulate an analogous question in gravity and ask what is the *nonlinear* tidal deformability of a compact object. In this work we make progress in this direction by explicitly computing the nonlinear Love numbers of Schwarzschild black holes at second-order in perturbation theory.

Nonlinear corrections to the Love numbers of black holes have been previously studied in [21, 47, 51–53].<fnref n="2" /> In particular, refs. [21, 51] provide evidence that the nonlinear Love numbers of Schwarzschild black holes vanish to all orders in perturbation theory, but their analysis is restricted to axisymmetric perturbations only. On the other hand, ref. [47] examines the quadratic tidal deformation of a non-rotating compact body using a post-Newtonian definition for the tidally induced multipole moments, focusing on the parity-even sector of perturbations. Moreover, in ref. [52] some of us showed that the electric-type quadrupolar Love numbers vanish at quadratic order by explicitly performing the matching with the point-particle Effective Field Theory (EFT), focusing on the parity-even sector of perturbations.

In this work, following up on and extending the methods proposed in [52], we present the first comprehensive study of quadratic corrections to the tidal deformability of Schwarzschild black holes in GR, integrating and expanding upon the results in the literature in multiple directions. First, as in [52], we adopt a definition for the quadratic Love numbers based on the worldline EFT [55, 56] (see also [10, 38, 39, 57–59] for some reviews). The EFT provides a robust definition of the (nonlinear) Love numbers as coupling constants of higher-dimensional operators, which allows to systematically address potential ambiguities resulting from gauge invariance in the theory and nonlinearities. In addition, our analysis now includes both even and odd perturbations and is not restricted to axisymmetric configurations. Finally, we complete the program initiated in [52] by computing the quadratic Love number couplings for all multipoles, beyond the leading quadrupole. We will show explicitly that quadratic Love numbers of asymptotically-flat Schwarzschild black holes in four-dimensional GR vanish exactly for generic multipoles.

In particular, in section 2, we derive the perturbation equations for the Schwarzschild metric up to second order in the Regge-Wheeler gauge, in the static limit, for both even and odd components. The quadratic terms in these equations describe how the coupling of linear perturbations sources the metric, respecting selection rules dictated by the symmetries of the setup, which we derive and discuss in detail. These equations are then used to derive the second-order metric in section 3, where we focus on the monopole, quadrupole, and hexadecapole induced by the quadratic coupling of two quadrupoles. Section 4 reviews the definition of quadratic Love numbers within the worldline EFT formalism and the background field method. In section 5, we use this to study the quadratic response of the Schwarzschild metric induced by an external field, which is then compared to the second-order solutions obtained in full GR. Section 6 is instead devoted to the analysis of second-order static perturbations of arbitrary multipole: first, we check that, thanks to the special structure

<footnote n="2">²See also [54] for a scalar-field analysis.</footnote>
of the source, the quadratic solutions take a simple polynomial form; we then perform the matching to the worldline EFT to all orders in the spatial derivatives, and conclude that all quadratic Love number couplings of Schwarzschild black holes vanish in GR.

More technical aspects of the calculations are finally collected in the appendices A–E. In particular, in appendix A we review the integrals between three tensor harmonics used in the main text, appendix B collects the explicit numerical values of some of the coefficients in the main text, appendix C explains how to express the tidal field in Cartesian coordinates, appendix D derives some useful relations between trace-free tensors and, finally, appendix E explains how to compute the Schwarzschild metric perturbatively within our approach.

**Notations and conventions.** We will work in units such that $c = \hbar = 1$ and use the mostly plus signature for the metric, $(-, +, +, +)$. The reduced Planck mass is defined by $M_{\text{Pl}}^{-2} \equiv 8\pi G$. For convenience, we will sometimes introduce the parameter $\kappa \equiv \sqrt{32\pi G} = 2M_{\text{Pl}}^{-1}$. For a black hole of mass $M$, we denote the Schwarzschild radius by $r_s = 2GM$. Our convention for the time Fourier transform is $\Psi(t, r, \theta, \phi) = \int \frac{d\omega}{2\pi} e^{-i\omega t} \tilde{\Psi}(\omega, r, \theta, \phi)$. Moreover, we will decompose quantities in spherical harmonics using the convention $\tilde{\Psi}(\omega, r, \theta, \phi) = \sum_{\ell,m} \tilde{\Psi}^{\ell m}(\omega, r) Y_{\ell}^{m}(\theta, \phi)$, where $Y_{\ell}^{m}(\theta, \phi)$ are normalized as $\int d\Omega Y_{\ell}^{m*}(\theta, \phi) Y_{\ell'}^{m'}(\theta, \phi) = \delta_{\ell\ell'}\delta^{mm'}$. To simplify the notation, we will often omit arguments and the tilde to denote the Fourier transform, relying on the context to discriminate between the different meanings. Greek indices run over spacetime coordinates, latin indices run over spatial coordinates and capital latin indices run over the angular variables $\theta, \phi$. We define the shorthand $\int_k \equiv \int \frac{d^4k}{(2\pi)^4}$, and the notation $\delta^{-(d)}(x) \equiv (2\pi)^d\delta^{(d)}(x)$. The symbol $(\cdots)$ denotes symmetrization over the enclosed indices e.g., $A_{(\mu}B_{\nu)} = \frac{1}{2}(A_{\mu}B_{\nu} + A_{\nu}B_{\mu})$, while we use $\langle \cdots \rangle$ to denote the traceless symmetrization. Depending on the context, we will use different metrics to raise and lower spacetime indices. The full metric $g_{\mu\nu}$ is used for the complete GR calculations in section 2, and when introducing the EFT worldline approach in section 4.1. In the EFT calculation, from section 4.2 to section 5.3, we will use the flat Minkowski metric $\eta_{\mu\nu}$ for this purpose. Additionally, we distinguish between the 2d Levi-Civita tensor $\epsilon$, whose indices are raised and lowered with the metric on the unit sphere, and the spacetime Levi-Civita symbol $\epsilon$, whose indices are raised and lowered with $\eta_{\mu\nu}$.

## 2 Static quadratic perturbations

Black hole perturbation theory has a long history, dating back to the original work by Regge and Wheeler on the study of linearized perturbations of Schwarzschild spacetime in the late 1950s [60]. More recently, the study of black hole perturbations has been extended to second order, see e.g. [61–70].

In this section we study the deviations from the Schwarzschild metric up to second order in the metric perturbations in the static limit. We first expand the Einstein-Hilbert (EH) action,

$$S_{\text{EH}} = \int d^4x \sqrt{-g} \frac{M_{\text{Pl}}^2}{2} R, \qquad (2.1)$$

where $R$ is the Ricci scalar, up to quadratic and cubic order in the metric perturbations. We then vary the action with respect to the metric perturbations, first at linear and then at second order, and we derive the second-order equations for the metric degrees of freedom.
## 2.1 Metric parametrization and gauge choice

We decompose the metric $g_{\mu\nu}$ as

$$g_{\mu\nu} = g_{\mu\nu}^{\text{Sch}} + \delta g_{\mu\nu}, \qquad (2.2)$$

where $g_{\mu\nu}^{\text{Sch}}$ is the static and spherically symmetric background solution. In spherical polar coordinates, $(t, r, \theta, \phi)$, this is given by the Schwarzschild metric,

$$g_{\mu\nu}^{\text{Sch}} = \begin{pmatrix} -f(r) & 0 & 0 & 0 \\ * & f^{-1}(r) & 0 & 0 \\ * & * & r^2 & 0 \\ * & * & * & r^2 \sin^2 \theta \end{pmatrix}, \qquad f(r) \equiv 1 - \frac{r_s}{r}, \qquad (2.3)$$

where the asterisks denote symmetric components and $r_s \equiv 2GM$, with M the mass of the black hole. In eq. (2.2), $\delta g_{\mu\nu}$ describes the deviation from the Schwarzschild metric.

To write $\delta g_{\mu\nu}$ in four spacetime dimensions, we follow the standard Regge-Wheeler (RW) parameterization [60]. It is convenient to classify the metric perturbations according to their behavior under parity transformation, $(\theta, \phi) \rightarrow (\pi - \theta, \phi + \pi)$. We thus distinguish between even, or polar, perturbations, and odd, or axial, perturbations. Moreover, we use the index + to label parity-even quantities and the index - for parity-odd quantities.

With this decomposition and notation, the most general parametrization of $\delta g_{\mu\nu}$ in four spacetime dimensions takes the form $\delta g_{\mu\nu} = \delta g_{\mu\nu}^{+} + \delta g_{\mu\nu}^{-}$, with

$$\delta g_{\mu\nu}^{+} = \begin{pmatrix} f(r)H_0 & H_1 & \nabla_A \mathcal{H}_0 \\ * & f^{-1}(r)H_2 & \nabla_A \mathcal{H}_1 \\ * & * & r^2 K \gamma_{AB} + r^2 \left(\nabla_A \nabla_B - \frac{1}{2} \gamma_{AB} \right) G \end{pmatrix}, \qquad (2.4)$$

$$\delta g_{\mu\nu}^{-} = \begin{pmatrix} 0 & 0 & -\epsilon_A^C \nabla_C h_0 \\ * & 0 & -\epsilon_A^C \nabla_C h_1 \\ * & * & -\frac{1}{2}(\epsilon_A^C \nabla_C \nabla_B + \epsilon_B^C \nabla_C \nabla_A)h_2 \end{pmatrix}. \qquad (2.5)$$

Here $\nabla_A$ and $\epsilon_{AB}$ are the covariant derivative and Levi-Civita tensor on a sphere, respectively, as reviewed in appendix A. Two-dimensional tensors labeled by capital latin letters, A, B, C, ..., are raised and lowered using the metric on the unit sphere $\gamma_{AB}$, with line element $\gamma_{AB}dx^A dx^B = d\theta^2 + \sin^2\theta d\phi^2$.

The quantities $H_a$, $\mathcal{H}_a$, K, G and $h_a$ above are functions of the coordinates $(t, r, \theta, \phi)$. However, since here we consider only the static limit, they are independent of time. Moreover, given the symmetries of the background solution $g_{\mu\nu}^{\text{Sch}}$, it is convenient to decompose these quantities in spherical harmonics. For instance, for $H_0$ we write

$$H_0(r, \theta, \phi) = \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} H_0^{\ell m}(r) Y_{\ell}^{m}(\theta, \phi), \qquad (2.6)$$

where $\Delta_S^2 Y_{\ell}^m \equiv \nabla^A \nabla_A Y_{\ell}^m = -\ell(\ell+1)Y_{\ell}^m$ and $i\partial_{\phi} Y_{\ell}^m = -mY_{\ell}^m$. While the even quantities $H_a$, $\mathcal{H}_a$, K, G appearing in $\delta g_{\mu\nu}^{+}$ do not transform under parity, $h_a$ flip sign under a parity transformation, ensuring that $g_{\mu\nu}$ transforms under parity as a 2-tensor. Therefore, from
the behavior of $Y_{\ell}^{m}(\theta, \phi)$ under parity, $Y_{\ell}^{m}(\pi - \theta, \phi + \pi) = (-1)^{\ell}Y_{\ell}^{m}(\theta, \phi)$, an even harmonic component, e.g. $H_{0}^{\ell m}(r)$, picks up a factor of $(-1)^{\ell}$ under parity transformation, while an odd component, e.g. $h_{0}^{\ell m}(r)$, picks up a factor of $(-1)^{\ell+1}$.

For later convenience, we rewrite eqs. (2.4) and (2.5) upon using the decomposition (2.6) and expressing the covariant derivatives in terms of partial derivatives. One obtains

$$ \delta g_{\mu\nu}^{+} = \begin{pmatrix} f H_0^{\ell m} & H_1^{\ell m} & H_0^{\ell m} \partial_{\theta} & H_0^{\ell m} \partial_{\phi} \\ * & f^{-1} H_2^{\ell m} & H_1^{\ell m} \partial_{\theta} & H_1^{\ell m} \partial_{\phi} \\ * & * & r^2 (K^{\ell m} + G^{\ell m} W) & r^2 G^{\ell m} V \\ * & * & * & r^2 \sin^2 \theta (K^{\ell m} + G^{\ell m} W) \end{pmatrix} Y_{\ell}^{m}, \quad (2.7) $$

$$ \delta g_{\mu\nu}^{-} = \begin{pmatrix} 0 & 0 & -\frac{1}{\sin\theta} h_0^{\ell m} \partial_{\phi} \sin\theta h_0^{\ell m} \partial_{\theta} \\ * & 0 & -\frac{1}{\sin\theta} h_1^{\ell m} \partial_{\phi} \sin\theta h_1^{\ell m} \partial_{\theta} \\ ** & * & -\frac{1}{\sin\theta} h_2^{\ell m} V \sin\theta h_2^{\ell m} W \\ ** & * & \sin\theta h_2^{\ell m} V \end{pmatrix} Y_{\ell}^{m}, \quad (2.8) $$

where the symbols $V$ and $W$ stand for differential operators defined as $V \equiv \partial_{\theta}\partial_{\phi} - \frac{\cos\theta}{\sin\theta}\partial_{\phi}$ and $W \equiv \frac{1}{2}(\partial_{\theta}^{2} - \frac{\cos\theta}{\sin\theta}\partial_{\theta} - \frac{1}{\sin^{2}\theta}\partial_{\phi}^{2})$, respectively.

To simplify the derivation of the perturbation equations and their solutions, it is useful to fix a gauge. A convenient choice, given the symmetries of the problem, is the RW gauge [60, 65, 66], which sets

$$ \mathcal{H}_0 = \mathcal{H}_1 = G = h_2 = 0. \quad (2.9) $$

This condition completely fixes the gauge freedom and can be directly applied to the action without losing any constraints [71].<fnref n="3" /> We then expand the action up to cubic order in the remaining metric variables and study the equations obtained from the variation of the action perturbatively.<fnref n="4" /> Specifically, for the metric variable $H_0$, we write $H_0 = {}^{(1)}H_0 + {}^{(2)}H_0 + \dots$, where ${}^{(1)}H_0$ is the solution to the linear equations, obtained by varying the quadratic action. ${}^{(2)}H_0$ is the solution to the second-order equations, derived by varying the cubic action and by using the linear solution in the terms that are quadratic in the metric perturbations. Note that the expansion parameter is the amplitude of the perturbations. Unlike the EFT calculation in section 5, no expansion is performed in $r_s/r$. In the following sections, we first study the linear order, and then the quadratic order.

## 2.2 Linear perturbations

Let us first briefly review the derivation of the linearized equations, obtained from the variation of the quadratic action with respect to the metric variables. Content and results

<footnote n="3">³A residual (large) gauge freedom remains when considering the lowest ($\ell = 0$ and $\ell = 1$) multipoles of the metric perturbations, as discussed in more details at the end of section 2.2 for linear perturbations and in section 3.3 for quadratic perturbations.</footnote>

<footnote n="4">⁴Note that, instead of going through the action, one could have chosen to work directly at the level of the Einstein equations, expanded up to second order in the fields. In the gauge (2.9), the two ways of proceeding are clearly equivalent and lead to the same result. We choose here to start from the action just for practical convenience [22, 72, 73]; having to deal with less equations, this facilitates the analysis at second order.</footnote>
of this subsection are well known; however, it will give us the opportunity to set up the notation and introduce the ingredients that will be necessary in section 2.3 for the analysis of the second-order perturbations.

There are four metric variables in the even sector, i.e. $H_0$, $H_1$, $H_2$ and $K$, and two in the odd sector, i.e. $h_0$ and $h_1$. However, not all of them are independent degrees of freedom. Indeed, in vacuum, the gravitational field is expected to describe two propagating degrees of freedom: one even and one odd. Since at linear order the even and odd components do not mix [60], we can treat them separately. Let us first focus on the even metric quantities. The four equations obtained by varying the action with respect to even quantities must contain three constraints and the dynamics of a single degree of freedom, as we are going to show next.

Instead of $H_2$, it will be convenient to use the variable $\psi$ obtained by substituting in the action [72, 74, 75]

$$H_2 = f(r)\psi + r\partial_r K + \frac{1}{f(r)}\left(1 - \frac{3r_s}{2r} + \frac{1}{2}\Delta_{S^2}\right) K, \quad (2.10)$$

where $\Delta_{S^2} \equiv \partial_{\theta}^2 + \frac{\cos\theta}{\sin\theta}\partial_{\theta} + \frac{1}{\sin^2\theta}\partial_{\phi}^2$. This field redefinition has the advantage of simplifying the linear field equations, making them of first order in $r$ rather than second order.

The linear equations are obtained by varying the quadratic action with respect to $H_0$, $\psi$, $K$ and $H_1$ and expanding in spherical harmonics. We first focus on $\ell \ge 2$; we discuss the monopole and dipole cases separately, below. Defining for convenience the shorthand notation

$$\lambda \equiv \ell(\ell + 1), \quad (2.11)$$

and using $\Delta_{S^2}Y_{\ell}^m = -\ell(\ell + 1)Y_{\ell}^m$ to replace the angular derivatives, we obtain

$$\partial_r {}^{(1)}\psi^{\ell m} + \frac{(\lambda+2)r+2r_s}{2f(r)r^2} {}^{(1)}\psi^{\ell m} - \frac{\lambda[(\lambda-2)r+3r_s]}{4f(r)^3r^2} {}^{(1)}K^{\ell m} = 0, \quad (2.12)$$

$$\partial_r {}^{(1)}K^{\ell m} - 2f(r)\partial_r {}^{(1)}H_0^{\ell m} - \frac{\lambda+1}{f(r)r} {}^{(1)}K^{\ell m} - \frac{\lambda}{r_s} {}^{(1)}H_0^{\ell m} + \frac{2f(r)}{r_s} {}^{(1)}\psi^{\ell m} = 0, \quad (2.13)$$

$$\partial_r {}^{(1)}\psi^{\ell m} - \frac{\lambda[(\lambda-2)r+3r_s]}{2f(r)^2rr_s} {}^{(1)}H_0^{\ell m} + \frac{\lambda[(\lambda-2)r(r-2r_s)-3r_s^2]}{2f(r)^3r^2r_s} {}^{(1)}K^{\ell m} + \frac{(\lambda+1)r+r_s}{f(r)r^2} {}^{(1)}\psi^{\ell m} = 0, \quad (2.14)$$

$$\lambda r^2 f(r) {}^{(1)}H_1^{\ell m} = 0. \quad (2.15)$$

As expected from spherical symmetry, there is no $m$ dependence in these equations.

The last equation simply implies ${}^{(1)}H_1 = 0$. Combining the first and third equation allows to remove $\partial_r {}^{(1)}\psi^{\ell m}$ and obtain a linear algebraic equation relating three variables,

$${}^{(1)}\psi^{\ell m} - \frac{(\lambda - 2)r + 3r_s}{f(r)r_s} {}^{(1)}H_0^{\ell m} - \frac{r(\lambda - 2)(3r_s - 2r) + 3r_s^2}{2f(r)^2rr_s} {}^{(1)}K^{\ell m} = 0. \quad (2.16)$$

We choose to describe the independent even degree of freedom in terms of the quantity $H_0$. Thus, this equation and its first derivative with respect to $r$, together with eqs. (2.12)
and (2.13), form a system of first-order differential equations which can be solved algebraically for $^{(1)}$ψ and $^{(1)}$K in terms of $^{(1)}$H$_{0}$ and $\partial_{r}$$^{(1)}$H$_{0}$, i.e.,

$$^{(1)}\psi^{lm} = -\frac{3(r_s - 2r)r_s^2 - r^2[rf(r)\lambda + r_s]}{2(\lambda - 2)r^3 f(r)^3} {}^{(1)}H_0^{lm} - \frac{[r(2r - 3r_s)(\lambda - 2) - 3r_s^2]}{2(\lambda - 2)rf(r)^2} \partial_r {}^{(1)}H_0^{lm}, \quad (2.17)$$

$$^{(1)}K^{lm} = -\frac{r_s^2 + rr_s(\lambda - 4) - r^2(\lambda - 2)}{(\lambda - 2)r^2 f(r)} {}^{(1)}H_0^{lm} + \frac{r_s}{\lambda - 2} \partial_r {}^{(1)}H_0^{lm}. \quad (2.18)$$

Therefore, replacing $^{(1)}$ψ and $^{(1)}$K (and their first derivatives) in eq. (2.16) by means of these expressions, we obtain a second-order equation for $^{(1)}$H$_{0}$,

$$\partial_r^2 {}^{(1)}H_0^{lm}(r) + \frac{2r - r_s}{r^2 f(r)} \partial_r {}^{(1)}H_0^{lm}(r) - \frac{\lambda r^2 f(r) + r_s^2}{r^4 f(r)^2} {}^{(1)}H_0^{lm}(r) = 0. \quad (2.19)$$

This equation describes the parity-even sector for linear perturbations; the other even metric quantities can be derived from the solution of this equation.

Similarly, the two equations obtained by varying the action with respect to odd quantities represent one constraint and one degree of freedom. Indeed, variation of the quadratic action with respect to $h_0$ and $h_1$ gives, respectively,

$$\partial_r^2 {}^{(1)}h_0^{lm}(r) - \frac{\lambda r - 2r_s}{r^3 f(r)} {}^{(1)}h_0^{lm}(r) = 0, \quad (2.20)$$

$$\lambda(\lambda - 2)rf(r) {}^{(1)}h_1^{lm}(r) = 0. \quad (2.21)$$

The second equation implies $^{(1)}h_1 = 0$, while the first equation describes the parity-odd linear perturbation. Note that $h_0$ and $h_1$ appear in the metric (2.5) with at least one angular derivative acting on them; therefore, there is no $l = 0$ odd-type mode and $\lambda \neq 0$.

In conclusion, in the static limit the parity-even and parity-odd sectors of the linear perturbations of a Schwarzschild black hole decouple and are described, respectively, by eq. (2.19) for $^{(1)}$H$_{0}$ and by eq. (2.20) for $^{(1)}$h$_{0}$. The two even variables $^{(1)}$H$_{2}$ (or $^{(1)}$ψ) and $^{(1)}$K are given in terms of $^{(1)}$H$_{0}$ from the constraint equations, while $^{(1)}$H$_{1}$ and $^{(1)}$h$_{1}$ vanish.

These considerations apply for $l \ge 2$. The cases $l = 0$ and $l = 1$ (i.e., $\lambda = 0$ and $\lambda = 2$, respectively) are different because some of the metric components in eqs. (2.7) and (2.8) are absent, leaving a residual gauge freedom, as discussed in [60] for the static case and in [76] for the general case. Specifically, for the monopole, the components $H_0$, $H_1$, and $G$ are missing and there are no odd perturbations, leaving us with only four even perturbations. Moreover, eqs. (2.12) and (2.14) become the same, while eq. (2.15) vanishes identically, reducing the system to only two equations instead of four. Using time and radial coordinate transformations, two metric components can be fixed, for instance $^{(1)}$K = $^{(1)}$H$_{1}$ = 0 [76]. The remaining components can be removed via a redefinition of the black hole mass, in accordance with Birkhoff's theorem.

For the dipole, the perturbations $G$ and $h_2$ are absent. Among the remaining even-type perturbations, three can be fixed by constraints, while the other three can be entirely gauged away by shifting the origin of the coordinate system. Of the two remaining odd perturbations, one can be eliminated through a coordinate transformation, and the other corresponds to a rotation of the black hole. For a non-rotating black hole, this can be set to zero. Therefore, in the following, we set all linear monopole and dipole perturbations to zero.
## 2.3 Quadratic perturbations

The second-order equations can be obtained by varying the cubic action with respect to the metric variables. The linear part of these equations is identical to the left-hand side of eqs. (2.12)–(2.15), (2.20) and (2.21), upon replacement of first-order quantities $^{(1)}$H$_{0}$, $^{(1)}$ψ, etc., by second-order quantities $^{(2)}$H$_{0}$, $^{(2)}$ψ, etc. Additionally, at second order there is a quadratic part, a source term, which describes how quadratic couplings of linear perturbations source the linear part of the equations. These quadratic terms include both even and odd quantities, i.e., they mix the two parity sectors. Therefore, at quadratic order, the even and odd components are no longer independent.

The procedure and steps to derive the equation describing the evolution of $^{(2)}$H$_{0}$ are analogous to the ones used to derive eq. (2.19) described above; the only difference is the presence of quadratic terms in the equation, involving products of linear quantities. Once expanded in spherical harmonics, these equations display tensor harmonics up to rank two. In particular, these are the even and odd vector harmonics, respectively defined as

$$Z_{\ell}^{m}{}_{A} \equiv \nabla_{A} Y_{\ell}^{m}, \qquad (2.22)$$

$$X_{\ell}^{m}{}_{A} \equiv \epsilon_{A}{}^{B} \nabla_{B} Y_{\ell}^{m}, \qquad (2.23)$$

and the even and odd 2-tensor harmonics, respectively defined as

$$Z_{\ell}^{m}{}_{AB} \equiv \left( \nabla_{A} \nabla_{B} Y_{\ell}^{m} - \frac{1}{2} \gamma_{AB} \Delta_{S^{2}} Y_{\ell}^{m} \right), \qquad (2.24)$$

$$X_{\ell}^{m}{}_{AB} \equiv \epsilon_{A}{}^{C} \left( \nabla_{C} \nabla_{B} Y_{\ell}^{m} - \frac{1}{2} \gamma_{CB} \Delta_{S^{2}} Y_{\ell}^{m} \right). \qquad (2.25)$$

These quantities form a basis, respectively of vectors and tensors, on the sphere.

To obtain an equation for the spherical harmonic coefficients, we multiply the resulting equations by $Y_{\ell}^{m}$ and integrate over the angles. Generally, the quadratic terms contain the following integrals,<fnref n="5" />

$$^{0}I_{\ell\ell_1\ell_2}^{mm_1m_2} = \int d\Omega Y_{\ell}^{m*}(\theta, \phi) Y_{\ell_1}^{m_1}(\theta, \phi) Y_{\ell_2}^{m_2}(\theta, \phi), \qquad (2.26)$$

$$^{1}I_{\ell\ell_1\ell_2}^{mm_1m_2} = \int d\Omega Y_{\ell}^{m*}(\theta, \phi) Z_{\ell_1 A}^{m_1}(\theta, \phi) Z_{\ell_2}^{m_2 A}(\theta, \phi), \qquad (2.27)$$

$$^{2}I_{\ell\ell_1\ell_2}^{mm_1m_2} = \int d\Omega Y_{\ell}^{m*}(\theta, \phi) Z_{\ell_1 AB}^{m_1}(\theta, \phi) Z_{\ell_2}^{m_2 AB}(\theta, \phi), \qquad (2.28)$$

$$^{1}J_{\ell\ell_1\ell_2}^{mm_1m_2} = \int d\Omega Y_{\ell}^{m*}(\theta, \phi) X_{\ell_1 A}^{m_1}(\theta, \phi) Z_{\ell_2}^{m_2 A}(\theta, \phi), \qquad (2.29)$$

$$^{2}J_{\ell\ell_1\ell_2}^{mm_1m_2} = \int d\Omega Y_{\ell}^{m*}(\theta, \phi) X_{\ell_1 AB}^{m_1}(\theta, \phi) Z_{\ell_2}^{m_2 AB}(\theta, \phi), \qquad (2.30)$$

where $d\Omega = \sin\theta d\theta d\phi$. We will explain how to solve these integrals in appendix A. There we show that

$$^{s}I_{\ell\ell_1\ell_2}^{mm_1m_2} = \frac{1}{2} [1 + (-1)^{\ell+\ell_1+\ell_2}] s I_{\ell\ell_1\ell_2}^{mm_1m_2} = \frac{1}{2} (s I_{\ell\ell_1\ell_2}^{mm_1m_2} + s I_{\ell\ell_2\ell_1}^{mm_2m_1}), \qquad (2.31)$$

$$^{s}J_{\ell\ell_1\ell_2}^{mm_1m_2} = -\frac{i}{2} [1 - (-1)^{\ell+\ell_1+\ell_2}] s I_{\ell\ell_1\ell_2}^{mm_1m_2} = -\frac{i}{2} (s I_{\ell\ell_1\ell_2}^{mm_1m_2} - s I_{\ell\ell_2\ell_1}^{mm_2m_1}), \qquad (2.32)$$

<footnote n="5">⁵In general, there are also integrals quadratic in odd tensors, such as $\int d\Omega Y_{\ell}^{m*}(\theta, \phi) X_{\ell_1 A}^{m_1}(\theta, \phi) X_{\ell_2}^{m_2 A}(\theta, \phi)$, but these can be reduced to integrals at most linear in odd tensors upon contractions of Levi-Civita symbols.</footnote>
where we defined

$$ {}^{s}\mathcal{I}_{\ell\ell_1\ell_2}^{m_1m_2} = (2 - \delta_{0s})2^{-s}(-1)^s \sqrt{\frac{(1+2\ell_1)(1+2\ell_2)}{4\pi(1+2\ell)}} \frac{(\ell_1+|s|)!(\ell_2+|s|)!}{(\ell_1-|s|)!(\ell_2-|s|)!} C_{\ell_1\ell_2\ell}^{m_1m_2\ell} C_{\ell_1\ell_2\ell}^{-ss0}, \quad (2.33) $$

and where $C_{\ell_1\ell_2\ell}^{m_1m_2\ell}$ are the Clebsch-Gordan coefficients, often denoted as $\langle \ell_1 m_1 \ell_2 m_2 | \ell m \rangle$. We note that transformation of spherical harmonics under parity implies that ${}^sI_{\ell\ell_1\ell_2}^{m_1m_2}$ vanishes for $\ell + \ell_1 + \ell_2$ odd, while ${}^sJ_{\ell\ell_1\ell_2}^{m_1m_2}$ vanishes for $\ell + \ell_1 + \ell_2$ even. In addition, from their definitions it is manifest that ${}^sI_{\ell\ell_1\ell_2}^{m_1m_2}$ is symmetric while ${}^sJ_{\ell\ell_1\ell_2}^{m_1m_2}$ is antisymmetric under the exchange of $\ell_1 m_1$ and $\ell_2 m_2$.

In the end, from the same manipulations that we used in the linear case, with the inclusion of the quadratic terms, we find the differential equation

$$ \partial_r^2 {}^{(2)}H_0^{\ell m}(r) + \frac{2r-r_s}{r^2 f(r)} \partial_r {}^{(2)}H_0^{\ell m}(r) - \frac{\lambda r^2 f(r) + r_s^2}{r^4 f(r)^2} {}^{(2)}H_0^{\ell m}(r) = \sum_{\ell_1, m_1, \ell_2, m_2} S_{H_0 \ell\ell_1\ell_2}^{m_1m_2}(r). \quad (2.34) $$

As expected, the left-hand side is the same as in the linear equation (2.19), with ${}^{(1)}H_0$ replaced by ${}^{(2)}H_0$. The source term on the right-hand side is given by a sum over angular and magnetic numbers $\ell_1, m_1, \ell_2$ and $m_2$, and the selection rules due to the Clebsch-Gordan coefficients imply that $|\ell_1 - \ell_2| \le \ell \le \ell_1 + \ell_2$, $|m_1| \le \ell_1$, $|m_2| \le \ell_2$ and $m = m_1 + m_2$. The source term is quadratic in both the even and odd linear harmonic variables ${}^{(1)}H_0^{\ell m}$, ${}^{(1)}\psi^{\ell m}$, ${}^{(1)}K^{\ell m}$ and ${}^{(1)}h_0^{\ell m}$, and their first and second radial derivatives, multiplied by one of the integrals in eqs. (2.26)–(2.30). The explicit expression of the source is cumbersome, but it can be simplified by using eqs. (2.17) and (2.18) to express ${}^{(1)}\psi$ and ${}^{(1)}K$ in terms of ${}^{(1)}H_0$ and $\partial_r {}^{(1)}H_0$. Additionally, one can use the linear equations (2.19) and (2.20) to eliminate the second-order radial derivatives and write the source term solely in terms of ${}^{(1)}H_0$, ${}^{(1)}h_0$, $\partial_r {}^{(1)}H_0$, and $\partial_r {}^{(1)}h_0$. By doing so, one obtains a lengthy but manageable expression,

$$ S_{H_0 \ell\ell_1\ell_2}^{m_1m_2}(r) = \frac{r_s^2}{\lambda(\lambda-2)f(r)^2 r^4} \left\{ \begin{aligned} & 0 I_{\ell\ell_1\ell_2}^{m_1m_2} \left[ \left( a_2 \frac{r^2}{r_s^2} + a_1 \frac{r}{r_s} + a_0 \right) {}^{(1)}H_0^{\ell_1 m_1} {}^{(1)}H_0^{\ell_2 m_2} \right. \\ & \qquad + f(r) \frac{r^2}{r_s} b_0 \partial_r {}^{(1)}H_0^{\ell_1 m_1} {}^{(1)}H_0^{\ell_2 m_2} - \frac{1}{2} \lambda(\lambda-2) f(r)^2 \frac{r^4}{r_s^2} \partial_r {}^{(1)}H_0^{\ell_1 m_1} \partial_r {}^{(1)}H_0^{\ell_2 m_2} \\ & \qquad + \frac{r_s^2}{f(r)^2 r^4} \left( c_4 \frac{r^4}{r_s^4} + c_3 \frac{r^3}{r_s^3} + c_2 \frac{r^2}{r_s^2} + c_1 \frac{r}{r_s} + c_0 \right) {}^{(1)}h_0^{\ell_1 m_1} {}^{(1)}h_0^{\ell_2 m_2} \\ & \qquad + \frac{r_s}{f(r)r^2} \left( d_2 \frac{r^2}{r_s^2} + d_1 \frac{r}{r_s} + d_0 \right) \partial_r {}^{(1)}h_0^{\ell_1 m_1} {}^{(1)}h_0^{\ell_2 m_2} \\ & \qquad + \left( e_2 \frac{r^2}{r_s^2} + e_1 \frac{r}{r_s} + e_0 \right) \partial_r {}^{(1)}h_0^{\ell_1 m_1} \partial_r {}^{(1)}h_0^{\ell_2 m_2} + (\ell_1 m_1 \leftrightarrow \ell_2 m_2) \end{aligned} \right. \\ & \qquad + {}^2 I_{\ell\ell_1\ell_2}^{m_1m_2} \left[ \frac{r_s^2}{f(r)^2 r^4} \left( f_3 \frac{r^3}{r_s^3} + f_2 \frac{r^2}{r_s^2} + f_1 \frac{r}{r_s} - \frac{19}{2} \right) {}^{(1)}h_0^{\ell_1 m_1} {}^{(1)}h_0^{\ell_2 m_2} \right. \\ & \qquad + \frac{r_s}{f(r)r^2} \left( g_2 \frac{r^2}{r_s^2} + g_1 \frac{r}{r_s} - 11 \right) \partial_r {}^{(1)}h_0^{\ell_1 m_1} {}^{(1)}h_0^{\ell_2 m_2} \\ & \qquad + \left. \left( 4 \frac{r}{r_s} - 2 \right) \partial_r {}^{(1)}h_0^{\ell_1 m_1} \partial_r {}^{(1)}h_0^{\ell_2 m_2} + (\ell_1 m_1 \leftrightarrow \ell_2 m_2) \right] \right\}, \quad (2.35) $$
where the coefficients $a_2, a_1$, etc., are functions of $\lambda, \lambda_1$ and $\lambda_2$, i.e., $a_2 = a_2(\lambda, \lambda_1, \lambda_2)$, etc., with $\lambda_1$ and $\lambda_2$ defined as $\lambda_i \equiv \ell_i(\ell_i + 1)$ in analogy with (2.11), with an extra subscript to distinguish the two different linear modes that combine to generate the second-order one. Explicit expressions for all the coefficients can be found in appendix B.

The source for $^{(2)}H_0$ contains both even and odd linear modes; however, only quadratic terms of the type even × even or odd × odd are present. As already noted, the symbols $sI_{\ell_1\ell_2}^{mm_1m_2}$ are symmetric under the exchange $\ell_1m_1 \leftrightarrow \ell_2m_2$ and vanish for $\ell + \ell_1 + \ell_2 = \text{odd}$. Therefore the source for $^{(2)}H_0$ receives contribution only when $\ell + \ell_1 + \ell_2$ is even and has to be symmetric under the exchange $\ell_1m_1 \leftrightarrow \ell_2m_2$. The other even variables of the second-order metric, $^{(2)}\psi$ and $^{(2)}K$, can be obtained from equations analogous to eqs. (2.17) and (2.18) extended at second order, which we do not write explicitly here. We have verified that, like $^{(2)}H_0$, they are sourced by even × even or odd × odd quadratic terms for $\ell + \ell_1 + \ell_2 = \text{even}$.

The remaining even variable, $H_1$, no longer vanishes at second order but it is determined by quadratic terms of the linear perturbations. Indeed, varying the cubic action with respect to $H_1$, multiplying the resulting equation by $Y_\ell^{m*}$ and integrating over the angles, yields

$$ ^{(2)}H_1^{\ell m} = \sum_{\ell_1, m_1, \ell_2, m_2} S_{H_1 \ell_1 \ell_2}^{m_1 m_2}(r), \qquad (2.36) $$

where the sources on the right-hand side are given by

$$ S_{H_1 \ell_1 \ell_2}^{m_1 m_2}(r) = \frac{1}{2\lambda f(r) r^2} J_{\ell_1 \ell_2}^{m_1 m_2} \left[ r_s^{(1)} H_0^{\ell_1 m_1(1)} h_0^{\ell_2 m_2} \right. \\ \left. \qquad + r^2 f(r) \left( \partial_r^{(1)} H_0^{\ell_1 m_1(1)} h_0^{\ell_2 m_2} - (1) H_0^{\ell_1 m_1} \partial_r^{(1)} h_0^{\ell_2 m_2} \right) - (\ell_1 m_1 \leftrightarrow \ell_2 m_2) \right]. \qquad (2.37) $$

The source above contains only quadratic terms of the type even × odd. The symbols $sJ_{\ell_1\ell_2}^{m_1m_2}$ are antisymmetric under the exchange $\ell_1m_1 \leftrightarrow \ell_2m_2$ and vanish whenever $\ell + \ell_1 + \ell_2$ even or $m_1 = m_2 = 0$ or $\ell_1m_1 = \ell_2m_2$ (see appendix A). Hence, $^{(2)}H_1$ does not receive any contribution in such cases. The only contribution can come from modes for which $\ell + \ell_1 + \ell_2$ is odd.

We can now turn to the odd quantities $h_0$ and $h_1$ at second order. The equations governing their behavior can be found by proceeding similarly to what we described above for the even sector. For $^{(2)}h_0$ we find

$$ \partial_r^{2(2)} h_0^{\ell m}(r) - \frac{\lambda r - 2r_s}{r^3 f(r)} \, ^{(2)}h_0^{\ell m}(r) = \sum_{\ell_1, m_1, \ell_2, m_2} S_{h_0 \ell_1 \ell_2}^{m_1 m_2}(r), \qquad (2.38) $$

where the left-hand side has the same form as in the linear equation (2.20), and the source

<footnote>$^6$The denominator in the prefactor of the curly brackets vanishes for $\lambda = 0$ and $\lambda = 2$, corresponding to $\ell = 0$ and $\ell = 1$, respectively. However, as the curly brackets also vanish in these cases, the source remains finite, as we will show in section 3.3. Similar apparent divergences occur in the other sources discussed in this section.</footnote>
<table><tr><td>ℓ + ℓ<sub>1</sub> + ℓ<sub>2</sub> even</td><td>(1) H<sub>0</sub><sup>ℓ<sub>2</sub>m<sub>2</sub></sup></td><td>(1) h<sub>0</sub><sup>ℓ<sub>2</sub>m<sub>2</sub></sup></td></tr><tr><td>(1) H<sub>0</sub><sup>ℓ<sub>1</sub>m<sub>1</sub></sup></td><td>(2) H<sub>0</sub><sup>ℓm</sup></td><td>(2) h<sub>0</sub><sup>ℓm</sup></td></tr><tr><td>(1) h<sub>0</sub><sup>ℓ<sub>1</sub>m<sub>1</sub></sup></td><td>(2) h<sub>0</sub><sup>ℓm</sup></td><td>(2) H<sub>0</sub><sup>ℓm</sup></td></tr></table>

<table><caption>**Table 1.** Selection rules for quadratic sources, for the case $ℓ + ℓ_1 + ℓ_2 = \text{even}$, left table, and $ℓ + ℓ_1 + ℓ_2 = \text{odd}$, right table. Here by $H_0$ we generically denote the trio $H_0, H_2$ and $K$, which are related by constraints.</caption><tr><td>ℓ + ℓ<sub>1</sub> + ℓ<sub>2</sub> odd</td><td>(1) H<sub>0</sub><sup>ℓ<sub>2</sub>m<sub>2</sub></sup></td><td>(1) h<sub>0</sub><sup>ℓ<sub>2</sub>m<sub>2</sub></sup></td></tr><tr><td>(1) H<sub>0</sub><sup>ℓ<sub>1</sub>m<sub>1</sub></sup></td><td>(2) h<sub>1</sub><sup>ℓm</sup></td><td>(2) H<sub>1</sub><sup>ℓm</sup></td></tr><tr><td>(1) h<sub>0</sub><sup>ℓ<sub>1</sub>m<sub>1</sub></sup></td><td>(2) H<sub>1</sub><sup>ℓm</sup></td><td>(2) h<sub>1</sub><sup>ℓm</sup></td></tr></table>

term is given by

$$ S_{h_0 \ell_1 \ell_2}^{m_1 m_2}(r) = \frac{1}{2\lambda f(r) r^2} \left\{ \begin{aligned} & 0 I_{\ell_1 \ell_2}^{m_1 m_2} \left[ \frac{r_s^2}{f(r) r^2} \left( \tilde{a}_2 \frac{r^2}{r_s^2} + \tilde{a}_1 \frac{r}{r_s} + \tilde{a}_0 \right) \right. \\ & \qquad \left. + r_s \left( \tilde{b}_1 \frac{r}{r_s} + \tilde{b}_0 \right) \partial_r \left( H_0^{\ell_1 m_1} H_0^{\ell_2 m_2} + (\ell_1 m_1 \leftrightarrow \ell_2 m_2) \right) \right] \\ & - 2 I_{\ell_1 \ell_2}^{m_1 m_2} \left[ H_0^{\ell_1 m_1} H_0^{\ell_2 m_2} + (\ell_1 m_1 \leftrightarrow \ell_2 m_2) \right] \end{aligned} \right\}, \quad (2.39) $$

where the explicit expressions of the coefficients $\tilde{a}_2, \tilde{a}_1, \tilde{a}_0, \tilde{b}_1$ and $\tilde{b}_0$ are given in appendix B. The odd quantity $^{(2)}h_0^{\ell m}$ is sourced by quadratic terms of the type even × odd, proportional to the $^s I_{\ell_1 \ell_2}^{m_1 m_2}$ symbols. Thus, by a similar logic as above, we conclude that $S_{h_0 \ell_1 \ell_2}^{m_1 m_2}$ must vanish for $ℓ + ℓ_1 + ℓ_2 = \text{odd}$ and must be symmetric under the exchange of $ℓ_1 m_1$ and $ℓ_2 m_2$.

The remaining second-order odd metric quantity, $^{(2)}h_1$, is determined by even × even and odd × odd quadratic terms, which are nonzero only for $ℓ + ℓ_1 + ℓ_2 = \text{odd}$ and antisymmetric under $ℓ_1 m_1 \leftrightarrow ℓ_2 m_2$. Indeed, we find

$$ ^{(2)}h_1^{\ell m} = \sum_{\ell_1, m_1, \ell_2, m_2} S_{h_1 \ell_1 \ell_2}^{m_1 m_2}(r), \quad (2.40) $$

with

$$ S_{h_1 \ell_1 \ell_2}^{m_1 m_2}(r) = \frac{1}{2\lambda(\lambda - 2)f(r)r} I_{\ell_1 \ell_2}^{m_1 m_2} \left[ r^3 f(r) {}^{(1)}H_0^{\ell_1 m_1} \partial_r {}^{(1)}H_0^{\ell_2 m_2} \right. \\ \left. \qquad + \lambda_2 {}^{(1)}h_0^{\ell_1 m_1} (2 - r\partial_r) {}^{(1)}h_0^{\ell_2 m_2} - (\ell_1 m_1 \leftrightarrow \ell_2 m_2) \right]. \quad (2.41) $$

## 2.4 Selection rules

Let us conclude this subsection with a brief summary of the selection rules encountered above. They arise from the symmetries of the system, specifically, the *staticity* (stationarity and time-reversal invariance) and *spherical symmetry* of the background, along with the invariance of the action under parity and time-reversal. As previously discussed, since $H_2$ and $K$ can be expressed in terms of $H_0$ and its derivatives using constraint equations, only four metric variables remain in the RW action, leading to 20 possible cubic combinations of metric variables (without counting derivatives) in the action. However, because $^{(1)}H_1 = 0 = {}^{(1)}h_1$,
cubic interactions involving at least two of these fields do not contribute to the quadratic equations of motion, excluding 10 cubic combinations. Additionally, due to the static background and the time-independence of the perturbations, it is impossible to construct cubic interactions with an odd number of $\delta g_{it}$ perturbations, such as $H_1 H_0^2$ or $h_0^3$. This leaves us with only five combinations, which, upon variation of the action, are the ones that yield the terms in the sources (2.35), (2.37), (2.39) and (2.41).

The spherical symmetry of the background allows for the decomposition of the fields into spherical harmonics, with angular momentum conservation in the cubic operators enforced by the Clebsch-Gordan coefficients. In general, the selection rules imposed by the Clebsch-Gordan coefficients dictate that $|\ell_1 - \ell_2| \le \ell \le \ell_1 + \ell_2$ and $m = m_1 + m_2$. If we define $p$ as the parity of a quadratic mode — where $p = +1$ denotes an even mode and $p = -1$ denotes an odd mode — then we have seen that $p = (-1)^{\ell+\ell_1+\ell_2} \times p_1 \times p_2$ [63, 64, 66, 67, 70], where $p_1$ and $p_2$ are the parities of the first and second linear modes generating the quadratic mode, respectively. Indeed, we can consider two cases:

* $\ell + \ell_1 + \ell_2 = \text{even}$: by parity and the properties of the integrals involving three tensor harmonics, this case arises from cubic couplings in the action with an even number of odd tensor harmonics, i.e. $H_0^3$ and $H_0 h_0^2$. In this case, quadratic terms with $p_1 = p_2 = \pm 1$ can source $p = +1$, corresponding to $^{(2)}H_0$ (or $^{(2)}H_2$ and $^{(2)}K$, using the constraints). Conversely, quadratic terms with $p_1 = -p_2 = \pm 1$ can source $p = -1$, corresponding to $^{(2)}h_0$. See table 1, left-hand side.

* $\ell + \ell_1 + \ell_2 = \text{odd}$: this case arises from the cubic couplings involving an odd number of odd tensor harmonics, i.e. $H_0^2 h_1$, $h_0^2 h_1$, and $H_0 H_1 h_0$. In this case, quadratic terms with $p_1 = p_2 = \pm 1$ can source $p = -1$, corresponding to $^{(2)}h_1$ (we recall that $^{(1)}H_1 = 0$, so that $^{(2)}h_0$ is not sourced at second order in this case), while quadratic terms with $p_1 = -p_2 = \pm 1$ can source $p = +1$, corresponding to $^{(2)}H_1$. See table 1, right-hand side. Additionally, as discussed above, the selection rules of the Clebsch-Gordan coefficients in this case imply that when $m_1 = m_2 = 0$ and $\ell_1 m_1 = \ell_2 m_2$, the quadratic sources vanish.

# 3 Solutions

In this section we solve the equations derived previously, upon imposing suitable boundary conditions at the horizon. We start by reviewing the derivation of the linear solutions.

## 3.1 Linear solutions

We solve here eqs. (2.19) and (2.20), respectively for $^{(1)}H_0^{\ell m}(r)$ and $^{(1)}h_0^{\ell m}(r)$, and find the regular solutions at the horizon $r = r_s$. The other even variables, $^{(1)}\psi$ and $^{(1)}K$, can be extracted using eqs. (2.17) and (2.18), once $^{(1)}H_0$ is known, while we remind that $^{(1)}H_1 = 0 = h_1$. Let us start with $^{(1)}H_0$. Whenever the $m$-dependence is trivial, we will suppress it from our notation in this section.

Equation (2.19) is a second-order ordinary differential equation with three regular singular points, i.e. $r = 0, r_s, \infty$. It is a special case of the hypergeometric equation, and it can be
recast in the form of an associated Legendre equation [77, 78]:

$$ (1 - z^2) \partial_z^2 {}^{(1)}H_0^{\ell m} - 2z \partial_z {}^{(1)}H_0^{\ell m} + \left[ \ell(\ell + 1) - \frac{4}{1 - z^2} \right] {}^{(1)}H_0^{\ell m} = 0, \quad z \equiv \frac{2r}{r_s} - 1. \quad (3.1) $$

The two independent solutions are thus<fnref n="7" />

$$ H_{0,1}^{\ell}(r) = P_{\ell}^{2} (2r/r_s - 1), \quad (3.2) $$

$$ H_{0,2}^{\ell}(r) = Q_{\ell}^{2} (2r/r_s - 1), \quad (3.3) $$

where $P_\ell^2$ and $Q_\ell^2$ are the associated Legendre functions of the first and second kind respectively, of degree $\ell$ and order 2, with near-horizon behaviors

$$ H_{0,1}^{\ell}(r) = -(\ell - 1)\ell(\ell + 1)(\ell + 2) \frac{r - r_s}{2r_s} + \mathcal{O}((r - r_s)^2), \quad (3.4) $$

$$ H_{0,2}^{\ell}(r) = -\frac{r_s}{2(r - r_s)} + \mathcal{O}((r - r_s)^0). \quad (3.5) $$

The second solution is divergent at the horizon, so we discard it, while the first one is regular. In particular, since $\ell$ takes integer values, $H_{0,1}^\ell$ is a finite polynomial in $r$, while $H_{0,2}^\ell$ contains instead $\log(\frac{r-r_s}{r})$ terms. Imposing regularity at the horizon leads us to discard the solution containing logarithms.<fnref n="8" /> Including a constant amplitude factor, we will write the physical linear solution for the even variable ${}^{(1)}H_0^{\ell m}$ as

$$ {}^{(1)}H_0^{\ell m}(r) = -\mathcal{E}_{\ell m}^{+} \frac{\Gamma(\ell+1)^2}{\ell(\ell-1)\Gamma(2\ell+1)} P_{\ell}^2(2r/r_s-1) = \mathcal{E}_{\ell m}^{+} \frac{\Gamma(\ell+1)^2}{\ell(\ell-1)\Gamma(2\ell+1)} r(r-r_s) \partial_r^2 P_{\ell}(2r/r_s-1), \quad (3.6) $$

where the normalization is chosen to ensure that $H_0^{\ell m}(r) = \mathcal{E}_{\ell m}^{+}(r/r_s)^\ell + \mathcal{O}(r^{\ell-1})$ for large $r$, with $\mathcal{E}_{\ell m}^{+}$ being the amplitude of the external even tidal field.

Let us now turn to ${}^{(1)}h_0$. Through the field redefinition

$$ {}^{(1)}h_0^{\ell m}(r(x)) = x^{-1}(1-x)u(x), \quad x \equiv \frac{r_s}{r}, \quad (3.7) $$

eq. (2.20) can be recast in the standard hypergeometric form

$$ x(1-x)u''(x) + [c - (a + b + 1)x] u'(x) - ab u(x) = 0, \quad (3.8) $$

with parameters

$$ a = -\ell - 1, \quad b = \ell, \quad c = -2, \quad (3.9) $$

$$ H_{0,1}^{\ell}(r) = \frac{\Gamma(\ell+2)\Gamma(\ell+3)}{\Gamma(\ell)} (-1)^{\ell+1} \frac{r}{r-r_s} {}_2F_1 \left(-\ell, \ell+1, 3; \frac{r}{r_s} \right), \\ H_{0,3}^{\ell}(r) = \left(\frac{r_s}{r}\right)^{\ell+1} \frac{r}{r-r_s} {}_2F_1 \left(\ell-1, \ell+1, 2\ell+2; \frac{r_s}{r} \right). $$

<footnote n="7">⁷In terms of hypergeometric functions, the two independent solutions can be written as [52]</footnote>

<footnote n="8">⁸One can express the Zerilli variable [76] in terms of $H_0$ and $\partial_r H_0$, and check that the logarithmic part of the homogeneous solution gives a divergent Zerilli variable [22].</footnote>
which satisfy $a + b - c = 1$. For this combination of parameters, the two independent solutions are<fnref n="9" />

$$h_{0,4}^{\ell}(r) = -\frac{r_s(r-r_s)}{r} \left(-\frac{r_s}{r}\right)^{\ell} {}_2F_1\left(\ell, \ell+3; 2\ell+2; \frac{r_s}{r}\right), \quad (3.10)$$

$$h_{0,5}^{\ell}(r) = -\frac{r^2(r-r_s)}{r_s^2} {}_2F_1\left(2-\ell, \ell+3; 4; \frac{r}{r_s}\right). \quad (3.11)$$

Since the first argument of ${}_2F_1(2-\ell, \ell+3; 4; r_s/r)$ in $h_{0,5}$ is a non-positive integer, while its third argument is a positive integer, it follows from the properties of the hypergeometric functions that $h_{0,5}$ is a finite polynomial with only positive powers of $r/r_s$.<fnref n="10" /> On the other hand, ${}_2F_1(\ell, \ell+3; 2\ell+2; r_s/r)$ in $h_{0,4}$ contains an inverse power of $r-r_s$ and a $\log(r-r_s)$, going as ${}_2F_1(\ell, \ell+3; 2\ell+2; r_s/r) = \frac{r_s}{r-r_s} \frac{\Gamma(2\ell+2)}{\Gamma(\ell)\Gamma(\ell+3)} + \frac{\Gamma(2\ell+2)}{\Gamma(\ell-1)\Gamma(\ell+2)} \log\left(\frac{r}{r_s}-1\right) + \mathcal{O}((r-r_s)^0)$ for $r \to r_s$. Note that, despite the divergent behavior of ${}_2F_1(\ell, \ell+3; 2\ell+2; r_s/r)$ near the horizon, $h_{0,4}$ remains finite because of the $r-r_s$ prefactor in (3.10). Following [22], we will demand that the Regge-Wheeler variable $\Psi_{RW}$, which is related to the odd field $h_0$ via $\Psi_{RW} \propto (2h_0 - r\partial_r h_0)$ [22, 76], does not diverge at the horizon. This excludes $h_{0,4}$, selecting $h_{0,5}$ as the correct physical odd solution, which we will rewrite for convenience as

$$^{(1)}h_0^{\ell m}(r) = \mathcal{E}_{\ell m}^{-} \frac{(4)_{\ell-2}\Gamma(\ell-1)}{(2-\ell)_{\ell-2}(\ell+3)_{\ell-2}} \frac{r^2(r-r_s)}{r_s^2} {}_2F_1\left(2-\ell, \ell+3; 4; \frac{r}{r_s}\right), \quad (3.13)$$

where $\mathcal{E}_{\ell m}^{-}$ is a dimensionless constant, corresponding to the amplitude of the linear odd field, and $(a)_n \equiv a(a+1)\cdots(a+n-1)$ is the Pochhammer symbol. The normalization factor is chosen in such a way that $^{(1)}h_0^{\ell}(r) = r_s \mathcal{E}_{\ell m}^{-}(r/r_s)^{\ell+1} + \mathcal{O}(r^\ell)$ for large $r$.

Equations (3.6) and (3.13) are the final expressions for the linear solutions that we will use in the next section to evaluate the quadratic sources in the second-order equations, to compute the metric perturbations at second order.

## 3.2 Quadratic solutions

We have now all the ingredients to solve eqs. (2.34) and (2.38) for the metric at second order, induced by the coupling of two linear modes. The right-hand sides of (2.34) and (2.38) are computed on the linear solutions (3.6) and (3.13), and are therefore known functions of $r$. From the standard properties of linear ordinary differential equations, the most general solutions to the inhomogeneous equations (2.34) and (2.38) are obtained from a linear superposition of homogeneous and particular solutions. The particular solutions can be found

$$h_{0,5}^{\ell}(r) = -\frac{r^2(r-r_s)}{r_s^2} \sum_{j=0}^{\ell-2} (-1)^j \frac{6(\ell-2)!}{j!(\ell-j-2)!} \frac{\Gamma(\ell+j+3)}{\Gamma(\ell+3)\Gamma(j+4)} \left(\frac{r}{r_s}\right)^j. \quad (3.12)$$

<footnote n="9">⁹See, e.g., ref. [79]. This corresponds to case 27 in the table in section 2.2.2 of that reference.</footnote>

<footnote n="10">¹⁰To see this, it might be convenient to rewrite $h_0$ in (3.11) (or, equivalently, (3.13) up to an overall constant prefactor) as</footnote>
using standard Green's function methods:

$$^{(2)}H_0^{lm}(r) = \sum_{\ell_1,m_1,\ell_2,m_2} \int_{r_s}^{\infty} dr' \mathcal{G}_{H_0}^{\ell}(r,r') \mathcal{S}_{H_0}^{m_1 m_2}{}^{\ell_1 \ell_2}(r'), \quad (3.14)$$

$$^{(2)}h_0^{lm}(r) = \sum_{\ell_1,m_1,\ell_2,m_2} \int_{r_s}^{\infty} dr' \mathcal{G}_{h_0}^{\ell}(r,r') \mathcal{S}_{h_0}^{m_1 m_2}{}^{\ell_1 \ell_2}(r'), \quad (3.15)$$

where $\mathcal{G}_{H_0}^{\ell}(r, r')$ and $\mathcal{G}_{h_0}^{\ell}(r, r')$ are the Green's functions associated with (2.34) and (2.38), satisfying

$$\partial_r^2 \mathcal{G}_{H_0}^{\ell}(r, r') + \frac{2r - r_s}{r^2 f(r)} \partial_r \mathcal{G}_{H_0}^{\ell}(r, r') - \frac{\lambda r^2 f(r) + r_s^2}{r^4 f(r)^2} \mathcal{G}_{H_0}^{\ell}(r, r') = \delta(r - r'), \quad (3.16)$$

$$\partial_r^2 \mathcal{G}_{h_0}^{\ell}(r, r') - \frac{\lambda r - 2r_s}{r^3 f(r)} \mathcal{G}_{h_0}^{\ell}(r, r') = \delta(r - r'), \quad (3.17)$$

respectively. Following the standard procedure, it is easy to check that the most general solutions for $\mathcal{G}_{H_0}^{\ell}(r, r')$ and $\mathcal{G}_{h_0}^{\ell}(r, r')$ that respect the correct boundary condition at the horizon, and are continuous across $r = r'$, are given by [52]:<fnref n="11" />

$$\mathcal{G}_{H_0}^{\ell}(r, r') = \frac{r'(r' - r_s)}{r_s W_{H_0}^{\ell}} \left[ H_{0,2}^{\ell}(r) H_{0,1}^{\ell}(r') \Theta(r - r') + H_{0,1}^{\ell}(r) H_{0,2}^{\ell}(r') \Theta(r' - r) \right], \quad (3.18)$$

$$\mathcal{G}_{h_0}^{\ell}(r, r') = \frac{1}{r_s W_{h_0}^{\ell}} \left[ h_{0,4}^{\ell}(r) h_{0,5}^{\ell}(r') \Theta(r - r') + h_{0,5}^{\ell}(r) h_{0,4}^{\ell}(r') \Theta(r' - r) \right], \quad (3.19)$$

where $\Theta$ is the Heaviside step function, while $H_{0,1}^{\ell}$, $H_{0,2}^{\ell}$, $h_{0,4}^{\ell}$ and $h_{0,5}^{\ell}$ are the homogeneous solutions (3.2)–(3.3) and (3.10)–(3.11). The quantities $W_{H_0}^{\ell}$ and $W_{h_0}^{\ell}$ in the expressions above are $r$-independent (but, $\ell$-dependent) dimensionless constants associated to the Wronskians of the differential equations, and can be expressed as

$$W_{H_0}^{\ell} = -\frac{1}{2}(\ell - 1)\ell(\ell + 1)(\ell + 2), \quad W_{h_0}^{\ell} = -\frac{3\Gamma(2\ell + 3)}{(\ell + 1)\Gamma(\ell + 2)\Gamma(\ell + 3)}. \quad (3.20)$$

The homogeneous solutions at second order to add to the particular ones are already known: they are given by eqs. (3.2)–(3.3) and (3.10)–(3.11), since the left-hand sides of (2.34) and (2.38) are identical to the left-hand sides in the linear equations (2.19) and (2.20). The corresponding arbitrary integration constants are chosen as follows [52, 54]: at second order, we will set to zero the constant associated to the solution that has growing profile at large $r$, since it is degenerate with the linear tidal field solution (in other words, it can always be reabsorbed through a redefinition of the linear tidal field's amplitude); we will choose the remaining constant, associated to the amplitude of the decaying solution at infinity, in such a way to remove any possible divergence at the horizon in the particular solutions (3.14) and (3.15). Since the Green's functions in (3.14) and (3.15) have already been chosen in (3.18) and (3.19) in such a way to respect regularity, we will also set the second integration constant to zero.

<footnote n="11">¹¹Recall that, in section 3.1, we argued that the static solutions associated with regular behavior at the horizon are $H_{0,1}^{\ell}(r)$ in (3.2) and $h_{0,5}^{\ell}(r)$ in (3.11).</footnote>
All in all, the inhomogeneous solutions for $^{(2)}$H$_0^{\ell m}$ and $^{(2)}$h$_0^{\ell m}$ that satisfy regularity conditions at the horizon can be written in full generality as [52]:

$$^{(2)}H_0^{\ell m}(r) = \frac{1}{r_s W_{H_0}^{\ell}} \sum_{\ell_1, m_1, \ell_2, m_2} \left[ H_{0,2}^{\ell}(r) \int_{r_s}^{r} dr' r'(r' - r_s) H_{0,1}^{\ell}(r') S_{H_0 \ell \ell_1 \ell_2}^{m_1 m_2}(r') - H_{0,1}^{\ell}(r) \int^{r} dr' r'(r' - r_s) H_{0,2}^{\ell}(r') S_{H_0 \ell \ell_1 \ell_2}^{m_1 m_2}(r') \right], \quad (3.21)$$

$$^{(2)}h_0^{\ell m}(r) = \frac{1}{r_s W_{h_0}^{\ell}} \sum_{\ell_1, m_1, \ell_2, m_2} \left[ h_{0,4}^{\ell}(r) \int_{r_s}^{r} dr' h_{0,5}^{\ell}(r') S_{h_0 \ell \ell_1 \ell_2}^{m_1 m_2}(r') - h_{0,5}^{\ell}(r) \int^{r} dr' h_{0,4}^{\ell}(r') S_{h_0 \ell \ell_1 \ell_2}^{m_1 m_2}(r') \right], \quad (3.22)$$

where the integrals on the second line of each equation are indefinite integrals.<fnref n="12" />

### 3.3 Quadratic solutions induced by quadrupole × quadrupole

It is instructive to provide an explicit example of the previous expressions. We will focus here on the case of two linear quadrupolar fields (i.e., $\ell_1 = 2$ and $\ell_2 = 2$), generating a second-order mode $\ell$ [52]. The $\ell = 2$ linear even and odd solutions can be read off from the general expressions (3.6) and (3.13):

$$^{(1)}H_0^{2m}(r) = \mathcal{E}_{2m}^{+} \frac{r}{r_s} \left( \frac{r}{r_s} - 1 \right), \quad (3.23)$$

$$^{(1)}h_0^{2m}(r) = \mathcal{E}_{2m}^{-} \frac{r^2}{r_s} \left( \frac{r}{r_s} - 1 \right), \quad (3.24)$$

where in the amplitudes of the even and odd linear modes, respectively $\mathcal{E}^+$ and $\mathcal{E}^-$, we kept the dependence on $m$ generic.

Based on the angular momentum and parity selection rules, the induced quadratic mode can only take the values $\ell = 0, \ell = 2$, or $\ell = 4$.<fnref n="13" /> In the following, we will discuss each case separately, starting with the quadrupole and addressing the somewhat special case of the monopole at the end. The analysis of the second-order solutions induced by generic $\ell_1$ and $\ell_2$ modes will be considered in section 6.

#### 3.3.1 Quadrupole

Let us start by analyzing the second-order $\ell = 2$ mode induced by two linear quadrupoles. In the even sector, from the general formula (3.21), we obtain for $^{(2)}H_0^{2m}$:

$$^{(2)}H_0^{2m} = \sum_{m_1, m_2} \mathcal{I}_{222}^{m_1 m_2} \frac{1}{4} \left[ \mathcal{E}_{2m_1}^{-} \mathcal{E}_{2m_2}^{-} \frac{3}{2} \frac{r}{r_s} \left( 12 \frac{r^3}{r_s^3} - 14 \frac{r^2}{r_s^2} + \frac{r}{r_s} - 1 \right) - \mathcal{E}_{2m_1}^{+} \mathcal{E}_{2m_2}^{+} \frac{r^2}{r_s^2} \left( 2 \frac{r^2}{r_s^2} + \frac{r}{r_s} - 3 \right) \right], \quad (3.25)$$

<footnote n="12">¹²Alternatively, they can also equivalently rewritten as integrals between $r$ and some other arbitrary radius $r_0$. The latter is however completely immaterial, because it can be removed in general via a redefinition of the integration constant of the growing homogeneous solutions $H_{0,1}^{\ell}(r)$ and $h_{0,5}^{\ell}(r)$ [52].</footnote>

<footnote n="13">¹³Quadratic perturbations with $\ell = 1$ and $\ell = 3$ vanish. This occurs because in these cases $\ell + \ell_1 + \ell_2$ is odd and, according to the selection rules discussed in table 1, only $^{(2)}H_1$ and $^{(2)}h_1$ can be generated. However, these also vanish because the radial profiles of their sources are antisymmetric under $\ell_1 m_1 \leftrightarrow \ell_2 m_2$ and do not explicitly depend on the magnetic quantum numbers.</footnote>
where the sums run over $|m_1| \le 2$ and $|m_2| \le 2$. By solving the constraint equations at second order one can also easily find the other even scalars in the metric,

$$^{(2)}H_2^{2m} = \sum_{m_1,m_2} {}^0I_{222}^{mm_1m_2} f(r) \frac{1}{4} \left[ \mathcal{E}_{2m_1}^{-} \mathcal{E}_{2m_2}^{-} \frac{3}{2} \frac{r}{r-r_s} \left( 56 \frac{r^3}{r_s^3} - 54 \frac{r^2}{r_s^2} + \frac{r}{r_s} - 1 \right) - \mathcal{E}_{2m_1}^{+} \mathcal{E}_{2m_2}^{+} \frac{r^3}{r_s^3} \left( 4 \frac{r}{r_s} + 1 \right) \right], \quad (3.26)$$

$$^{(2)}K^{2m} = \sum_{m_1,m_2} {}^0I_{222}^{mm_1m_2} \frac{1}{16} \left[ \mathcal{E}_{2m_1}^{-} \mathcal{E}_{2m_2}^{-} 3 \left( 42 \frac{r^4}{r_s^4} - 24 \frac{r^3}{r_s^3} + 2 \frac{r^2}{r_s^2} - 1 \right) - \mathcal{E}_{2m_1}^{+} \mathcal{E}_{2m_2}^{+} \left( 2 \frac{r^4}{r_s^4} + 8 \frac{r^2}{r_s^2} - 7 \right) \right]. \quad (3.27)$$

Since $\ell + \ell_1 + \ell_2$ is even, $^{(2)}H_1^{2m} = 0$, for all $m$.

In the odd sector, from (3.22) and eq. (2.39) for the source, we find instead

$$^{(2)}h_0^{2m} = - \sum_{m_1,m_2} {}^0I_{222}^{mm_1m_2} \mathcal{E}_{2m_1}^{-} \mathcal{E}_{2m_2}^{+} \frac{1}{2} \frac{r^3}{r_s^2} \left( \frac{r^2}{r_s^2} + \frac{r}{r_s} - 2 \right). \quad (3.28)$$

In addition, the other odd variable in the metric vanishes, $^{(2)}h_1^{2m} = 0$.

### 3.3.2 Hexadecapole

For the induced $\ell = 4$ mode, the derivation follows the same logic as for the quadrupole. The second-order solutions are

$$^{(2)}H_0^{4m} = \sum_{m_1,m_2} {}^0I_{422}^{mm_1m_2} \frac{1}{36} \left[ \begin{aligned} & \mathcal{E}_{2m_1}^{+} \mathcal{E}_{2m_2}^{+} \frac{r^2}{r_s^2} \left( 73 \frac{r^2}{r_s^2} - 107 \frac{r}{r_s} + 34 \right) \\ & + \mathcal{E}_{2m_1}^{-} \mathcal{E}_{2m_2}^{-} \frac{r}{r_s} \left( -139 \frac{r^3}{r_s^3} + 119 \frac{r^2}{r_s^2} + 17 \frac{r}{r_s} - 3 \right) \end{aligned} \right], \quad (3.29)$$

$$^{(2)}H_2^{4m} = \sum_{m_1,m_2} {}^0I_{422}^{mm_1m_2} f(r) \frac{1}{36} \left[ \begin{aligned} & \mathcal{E}_{2m_1}^{+} \mathcal{E}_{2m_2}^{+} \frac{r^3}{r_s^3} \left( 97 \frac{r}{r_s} - 58 \right) \\ & - \mathcal{E}_{2m_1}^{-} \mathcal{E}_{2m_2}^{-} \frac{r^2}{r_s(r-r_s)} \left( 7 \frac{r^3}{r_s^3} + \frac{r^2}{r_s^2} - 17 \frac{r}{r_s} + 3 \right) \end{aligned} \right], \quad (3.30)$$

$$^{(2)}K^{4m} = \sum_{m_1,m_2} {}^0I_{422}^{mm_1m_2} \frac{1}{108} \left[ \begin{aligned} & \mathcal{E}_{2m_1}^{+} \mathcal{E}_{2m_2}^{+} \left( 291 \frac{r^4}{r_s^4} - 266 \frac{r^3}{r_s^3} - 12 \frac{r^2}{r_s^2} + 21 \right) \\ & - \mathcal{E}_{2m_1}^{-} \mathcal{E}_{2m_2}^{-} \left( 21 \frac{r^4}{r_s^4} + 52 \frac{r^3}{r_s^3} - 30 \frac{r^2}{r_s^2} + 1 \right) \end{aligned} \right], \quad (3.31)$$

$$^{(2)}h_0^{4m} = \sum_{m_1,m_2} {}^0I_{422}^{mm_1m_2} \mathcal{E}_{2m_1}^{-} \mathcal{E}_{2m_2}^{+} \frac{1}{5} \frac{r^3}{r_s^2} \left( 8 \frac{r^2}{r_s^2} - 13 \frac{r}{r_s} + 5 \right). \quad (3.32)$$

By parity, $^{(2)}H_1^{4m} = 0 = {}^{(2)}h_1^{4m}$.

### 3.3.3 Monopole

As discussed at the end of section 2.2, the case of the monopole is special because the metric before gauge fixing contains only four (even) components: $H_0$, $H_1$, $H_2$ and $K$. At the linear level, a monopolar perturbation can be absorbed into the definition of the Schwarzschild mass, consistently with Birkhoff's theorem. However, this is no longer the case at second
order because the two linear quadrupoles, $\ell_1 = 2$ and $\ell_2 = 2$, break the spherical symmetry. Therefore, the generated quadratic monopole in RW gauge is a physical perturbation.

Despite its apparent singular behavior, the source in RW gauge in eq. (2.35) remains finite for $\ell = 0$ once the linear solutions are inserted and the limit $\lambda \to 0$ is appropriately taken. Therefore, the quadratic equation (2.34) can be used to compute the monopole of $^{(2)}H_0$. The result is

$$ ^{(2)}H_0^{00} = \sum_m {}^0I_{022}^{0-mm} \left[ \begin{aligned} & \mathcal{E}_{2-m}^{+} \mathcal{E}_{2m}^{+} \frac{r^2}{2r_s(r-r_s)} \left( -\frac{r^3}{r_s^3} + \frac{r^2}{r_s^2} - 2\frac{r}{r_s} - 3 \right) \\ & + \mathcal{E}_{2-m}^{-} \mathcal{E}_{2m}^{-} \frac{3r_s}{4r} \left( 6\frac{r^5}{r_s^5} + 2\frac{r^4}{r_s^4} - 6\frac{r^3}{r_s^3} - \frac{r}{r_s} + 1 \right) \end{aligned} \right]. \quad (3.33) $$

Moreover, the second-order constraint equations that relate $H_2$ and $K$ to $H_0$ can also be extended for $\lambda = 0$, yielding

$$ ^{(2)}H_2^{00} = \sum_m {}^0I_{022}^{0-mm} \frac{1}{f(r)} \frac{1}{2} \left[ \begin{aligned} & \mathcal{E}_{2-m}^{+} \mathcal{E}_{2m}^{+} \frac{r^2}{r_s(r-r_s)} \left( 7\frac{r^3}{r_s^3} - 23\frac{r^2}{r_s^2} + 26\frac{r}{r_s} - 11 \right) \\ & - \mathcal{E}_{2-m}^{-} \mathcal{E}_{2m}^{-} \frac{3r_s}{2r} \left( 98\frac{r^5}{r_s^5} - 138\frac{r^4}{r_s^4} + 42\frac{r^3}{r_s^3} + \frac{r}{r_s} - 1 \right) \end{aligned} \right], \quad (3.34) $$

$$ ^{(2)}K^{00} = \sum_m {}^0I_{022}^{0-mm} \frac{1}{2} \left[ \begin{aligned} & \mathcal{E}_{2-m}^{+} \mathcal{E}_{2m}^{+} \left( 2\frac{r^4}{r_s^4} - 4\frac{r^3}{r_s^3} + 2\frac{r^2}{r-s^2} - 1 \right) \\ & + \mathcal{E}_{2-m}^{-} \mathcal{E}_{2m}^{-} \frac{3r_s}{2r} \left( -20\frac{r^5}{r_s^5} + 32\frac{r^4}{r_s^4} - 12\frac{r^3}{r_s^3} - \frac{r}{r_s} + 2 \right) \end{aligned} \right], \quad (3.35) $$

while $^{(2)}H_1^{00} = 0$.

It is important to note that this procedure has implicitly fixed a gauge. Notice that when $\ell = 0$ the metric components corresponding to $H_0$, $H_1$ and $G$ vanish trivially, see eq. (2.7). Therefore, (2.9) does not fix any gauge: one can still perform a gauge transformation with the freedom of choosing two arbitrary functions of $r$. In the way we took the limit, such functions have been implicitly chosen to yield precisely the solutions written above. To understand this better, let us follow a different logic. Instead of taking the limit $\ell \to 0$ in the general-$\ell$ results, let us set $\ell = 0$ from the very beginning. In particular, let us choose for instance the gauge such that $^{(2)}K^{00} = 0$ and $^{(2)}H_1^{00} = 0$. One shall then follow the standard procedure, and find the solutions for $^{(2)}H_0^{00}$ and $^{(2)}H_2^{00}$. One can verify explicitly that these solutions are indeed related to (3.33)–(3.35) via a second-order gauge transformation, up to a homogeneous solution which can be reabsorbed into a redefinition of the Schwarzschild radius in the background metric.

# 4 Worldline effective field theory

In this section we leave the framework of full GR to enter that of the EFT. First introduced in [55, 56], the worldline EFT allows for a clear and unambiguous definition of the Love numbers as Wilson coefficients of higher-dimensional operators. We will first briefly review its main features (see refs. [10, 38, 39, 57–59] for more details), which we will extend to account for nonlinear response. We will then introduce the background field method [80–82]
and explain its advantages for our calculations. In this context, the metric perturbation is computed using Feynman diagrams, which we will provide at the end of this section, and solving the resulting (loop) integrals.

Only in the next section we will perform a similar calculation to the one discussed in the previous section, but within the EFT context. By matching these results with those obtained in full GR, we will be able to uniquely determine the values of the Love number couplings induced by the quadrupole-quadrupole interaction. The extension to other multipoles will be considered in section 6.

## 4.1 Point-particle action and nonlinear Love numbers

From far away, any object appears in first approximation as a point source. For non-rotating bodies, the action describing the dynamics of a point particle in GR is

$$S = S_{\text{EH}} + S_{\text{pp}}, \qquad (4.1)$$

where $S_{\text{EH}}$ is the usual Einstein-Hilbert term, while $S_{\text{pp}}$ is the point-particle action,

$$S_{\text{pp}} = -M \int d\tau = -M \int \sqrt{-g_{\mu\nu}} dx^{\mu} dx^{\nu}, \qquad (4.2)$$

with $M$ the mass of the source, and $\tau$ the point particle's proper time. In the second equality we used that $d\tau = \sqrt{-g_{\mu\nu}}dx^{\mu}dx^{\nu}$. In terms of the proper time, we shall define the four-velocity $u^{\mu} = dx^{\mu}/d\tau$.

Clearly, for a generic object eq. (4.1) cannot be exact. Finite-size effects and the tidal deformation of the object can be described by adding to (4.1) a series of higher-dimensional operators, organized in a derivative expansion and localized on the particle's worldline. At quadratic order, they take the form [56, 83]

$$S_{\text{T}} = \sum_{\ell=2}^{\infty} \int d\tau (Q_{E}^{\mu L} E_{\mu L} + Q_{B}^{\mu L} B_{\mu L}) + \dots, \qquad (4.3)$$

where ellipses stand for nonlinear couplings between $Q_{E,B}$ and powers of $E$ and $B$, and where we have introduced the multi-index notation $\mu_L \equiv \mu_1 \cdots \mu_\ell$ and the operators [52, 84]

$$E_{\mu_1 \cdots \mu_\ell} \equiv \nabla_{\langle \mu_1}^{\perp} \cdots \nabla_{\mu_{\ell-2}}^{\perp} E_{\mu_{\ell-1} \mu_\ell \rangle}, \qquad (4.4)$$

$$B_{\mu_1 \cdots \mu_\ell} \equiv \nabla_{\langle \mu_1}^{\perp} \cdots \nabla_{\mu_{\ell-2}}^{\perp} B_{\mu_{\ell-1} \mu_\ell \rangle}, \qquad (4.5)$$

with

$$E_{\mu\nu} \equiv C_{\mu\rho\nu\sigma} u^{\rho} u^{\sigma}, \qquad B_{\mu\nu} \equiv \frac{1}{2} \epsilon_{\gamma\alpha\beta\langle\mu} C_{\nu\rangle\delta}^{\alpha\beta} u^{\delta} u^{\gamma}, \qquad (4.6)$$

corresponding to the electric and magnetic parts of the Weyl tensor $C_{\mu\rho\nu\sigma}$, respectively. The symbol $\langle \cdots \rangle$ denotes the traceless symmetrization over the enclosed indices, while

$$\nabla_{\mu}^{\perp} \equiv (\delta_{\mu}^{\nu} + u^{\nu} u_{\mu}) \nabla_{\nu} \qquad (4.7)$$

is the covariant derivative projected on the hypersurface orthogonal to $u^{\mu}$. In the following, we will mostly work in the rest frame of the point particle, such that the only nontrivial
components of $E_{\mu\nu}$ and $B_{\mu\nu}$ become the spatial ones, i.e. $E_{ij} = C_{titj}$ and $B_{ij} \equiv \frac{1}{2}\epsilon_{tkl}\langle iC_{j} \rangle^{kl}$. In practice, we will replace $\mu_L$ in (4.3) with the spatial indices $i_L \equiv i_1 \cdots i_\ell$.

In the action (4.3), the quantities $Q_{E,B}^{\mu L}$ carry all the information about the response of the body, as a result of the interaction between the worldline and the external fields $E_{\mu\nu}$ and $B_{\mu\nu}$. Such effective couplings, which have been kept fully general so far, do not capture only the conservative tidal response of the body, but can be used also to describe dissipative effects. In the latter case, the quantities $Q_{E,B}^{\mu L}$ are interpreted as composite operators that depend on additional gapless degrees of freedom, which live on the worldline and parametrize the unknown microscopic physics that is responsible for dissipation, such as the internal dynamics of the object, or absorption across the horizon in the case of black holes [56, 85, 86]. In the spirit of the EFT, we will remain completely agnostic about the microscopic details of the system, which are integrated out, and we will parametrize $Q_{E,B}^{\mu L}$ in (nonlinear) response theory.

Concretely, one can express the one-point function of $Q_{E,B}^{\mu L}$ in the most general way as

$$ \langle Q_{E,B}^{i_L}(\tau) \rangle = \sum_{n=1}^{\infty} \sum_{k=0}^{n} \int d\tau_1 \cdots d\tau_n \langle {}^{(n)}\mathcal{R}_{E,B}^{i_L|j_L_1 \cdots j_L_k|j_{L_{k+1}} \cdots j_{L_n}} (\tau - \tau_1, \ldots, \tau - \tau_n) \\ \times E_{j_{L_1}}(\tau_1) \cdots E_{j_{L_k}}(\tau_k) B_{j_{L_{k+1}}} (\tau_{k+1}) \cdots B_{j_{L_n}}(\tau_n), \quad (4.8) $$

where ${}^{(n)}\mathcal{R}_{E,B}$ represents the $n^{th}$-order response function. As an example, when $n=1$, eq. (4.8) gives the linear response of the system and ${}^{(1)}\mathcal{R}_{E,B}$ reduces to the standard retarded Green's function [85, 86]. For higher values of $n$, eq. (4.8) captures instead the nonlinear response.<fnref n="14" /> This is similar to nonlinear polarization theory in electromagnetism, where the nonlinear polarization of an optical medium is parametrized as an expansion in powers of the external field (see, e.g., [50]). Note also that, as opposed to the linear case where the even and odd sectors are decoupled, at higher nonlinear orders an electric response (analogously, a magnetic response) can be induced by combinations of both external electric and magnetic fields [52, 84, 87]. Besides obvious symmetry properties under the exchange of its multi-indices, there are some nontrivial constraints on the form of the tensor ${}^{(n)}\mathcal{R}$. One comes from causality, which forces ${}^{(n)}\mathcal{R}$ to vanish whenever any of its arguments $\tau - \tau_j$, for $j=1, \dots, n$, turns negative. Moreover, the tensorial structure of ${}^{(n)}\mathcal{R}$ is further dictated by the building blocks in the theory: these are the metric, the Levi-Civita tensor and possibly the spin vector, in the case of rotating objects.

Since we are interested here in performing the matching with Schwarzschild black holes at leading order in the static limit, a certain number of simplifications occur. In the absence of rotation, the tensorial form of ${}^{(n)}\mathcal{R}$ boils down to simply a product of Kronecker deltas, involving only spatial indices in the rest frame of the point particle [52, 84]. In addition, in the limit of static tides, the system is non-dissipative and the time dependence in ${}^{(n)}\mathcal{R}$ reduces to a product of Dirac delta-functions [48, 56, 85, 86]. For instance, to leading order in the number of gradients, one has, schematically,

$$ {}^{(n)}\mathcal{R}^{ij|i_1 j_1 \cdots i_n j_n} = \sum_k \lambda_k \text{perm}_k \left[ \delta_{|j_n|}^{\langle i_j} \delta_{\langle i_1}^{j_i} \delta_{j_1}^{\langle i_2} \delta_{\langle i_3}^{j_3} \cdots \delta_{\langle i_n|}^{j_{n-1}} \right] \delta(\tau - \tau_1) \cdots \delta(\tau - \tau_n), \quad (4.9) $$

<footnote n="14">¹⁴The case $n=0$ captures the intrinsic multipoles of the object, which are irrelevant for our discussion below.</footnote>
where the sum runs over non-redundant permutation $k$ of the indices (see, e.g., discussion in ref. [84]; see also section 6.2 for higher orders in derivatives). Plugging the nonlinear response theory solution (4.8) into the action (4.3) including nonlinear terms, one then finds a series of effective interactions which, under the previous assumptions, take the form of local operators involving contractions of $E_{ij}$ and $B_{ij}$, and derivatives thereof.

For concreteness, up to cubic order in the number of fields and to leading order in derivatives, the finite-size action $S_T$ reads

$$S_T = \int d\tau \left[ \lambda_2^{(E)} E_{ij} E^{ij} + \lambda_2^{(B)} B_{ij} B^{ij} \right. \\ \left. \qquad + \lambda_{222}^{(E^3)} E^i_j E^j_k E^k_i + \lambda_{222}^{(B^3)} B^i_j B^j_k B^k_i + \lambda_{222}^{(E^2 B)} E^i_j E^j_k B^k_i + \lambda_{222}^{(EB^2)} E^i_j B^j_k B^k_i \right]. \qquad (4.10)$$

The coefficients $\lambda_2^{(E)}$ and $\lambda_2^{(B)}$ are the linear Love number couplings, which are famously zero for black holes in GR [14–16, 20–22, 24–27, 40, 46, 88]. In ref. [52], we verified that this result is robust against nonlinear corrections, and we further showed that $\lambda_{222}^{(E^3)} = 0$. On the other hand, for parity reasons, one immediately concludes that $\lambda_{222}^{(B^3)} = \lambda_{222}^{(E^2 B)} = 0$ in GR [87]. In section 5, we will extend our previous result [52] to the odd sector: we will explicitly perform the matching at the next-to-leading order in powers of $r_s$ including the operator $E^i_j B^j_k B^k_i$, and show that $\lambda_{222}^{(EB^2)} = 0$ for Schwarzschild black holes. We will finally generalize the matching to all orders in the number of spatial derivatives in section 6.

## 4.2 Background field method

Given the EFT presented in the previous section, the starting point of our computation is the following action

$$S[g, x(\tau)] = S_{\text{EH}}[g] + S_{\text{pp}}[g, x(\tau)] + S_{\text{T}}[g, x(\tau)], \qquad (4.11)$$

where the terms on the right-hand side are given in eqs. (2.1), (4.2) and (4.10), respectively. The goal is to compute the response induced by an external field. More precisely, we separate the full metric as

$$g_{\mu\nu} = \bar{g}_{\mu\nu} + \kappa h_{\mu\nu}, \qquad (4.12)$$

where $\kappa \equiv \sqrt{32\pi G}$. Here $\bar{g}_{\mu\nu}$ is an external “background” metric which includes the tidal field and, therefore, has to satisfy the full nonlinear Einstein equations in vacuum, i.e.,

$$\bar{R}_{\mu\nu} - \frac{1}{2}\bar{g}_{\mu\nu}\bar{R} = 0, \qquad (4.13)$$

with $\bar{R}_{\mu\nu}$ and $\bar{R}$ representing the Ricci tensor and Ricci scalar of $\bar{g}_{\mu\nu}$, respectively. The second term, $\kappa h_{\mu\nu}$, denotes the response.

Given this split, it is convenient to use the background field method [80–82]. Specifically, we insert $g_{\mu\nu}$ from eq. (4.12) into eq. (4.11) and rewrite the action into two pieces, $S[g, x(\tau)] = S_0[\bar{g}, x(\tau)] + S_h[\bar{g}, h, x(\tau)]$, where the first term, $S_0$, is independent of the response field $h_{\mu\nu}$.
Then, we treat the massive object as an external source of the gravitational field. This allows us to compute the effective action of the system, $W[J]$, by formally performing the following path integral,

$$e^{iW[J]} = e^{iS_0[\bar{g},x(\tau)]} \int \mathcal{D}[h] \exp \left\{ iS_h[\bar{g},h, x(\tau)] + iS_{\text{GF}}[\bar{g},h] + i \int d^4y \sqrt{-\bar{g}} \mathcal{J}^{\mu\nu}(y) h_{\mu\nu}(y) \right\}. \quad (4.14)$$

The additional action, $S_{\text{GF}}$, is the usual gauge-fixing term arising from a Faddeev-Popov procedure. We will work in the background de Donder gauge, i.e. we define

$$S_{\text{GF}} \equiv - \int d^4 x \sqrt{-\bar{g}} \bar{g}^{\mu\nu} \bar{\Gamma}_{\mu} \bar{\Gamma}_{\nu}, \qquad \bar{\Gamma}_{\mu} \equiv \bar{g}^{\rho\sigma} \bar{\nabla}_{\rho} h_{\sigma\nu} - \frac{\bar{g}^{\rho\sigma}}{2} \bar{\nabla}_{\nu} h_{\rho\sigma}, \quad (4.15)$$

where $\bar{\nabla}_{\rho} \bar{g}_{\mu\nu} = 0$. This choice ensures that the final result is invariant under diffeomorphisms of the external metric $\bar{g}_{\mu\nu}$.

We are interested in computing the one-point expectation value of the field $h_{\mu\nu}$, which in terms of the effective action $W[J]$ above is given by<fnref n="15" />

$$\langle h_{\mu\nu}(x) \rangle = \frac{1}{\sqrt{-\bar{g}}} \frac{\delta W[J]}{i \delta \mathcal{J}^{\mu\nu}(x)} \Big|_{\mathcal{J}=0}. \quad (4.16)$$

In practice, we will compute the above equation by considering all Feynman diagrams with one external $h_{\mu\nu}(x)$ field. Although we are performing a path integral, we are ultimately interested in the classical limit. Following [55], we enforce this limit from the outset by employing a saddle-point approximation. Diagrammatically, this means neglecting all diagrams that contain close loops of the field $h_{\mu\nu}$. This is also why we do not need to add any ghost fields in eq. (4.14).

One last step before proceeding is to write

$$\bar{g}_{\mu\nu} = \eta_{\mu\nu} + H_{\mu\nu}, \quad (4.17)$$

where $\eta_{\mu\nu}$ is the flat Minkowski metric and $H_{\mu\nu}$ is the external (static) tidal field, which is determined by solving the vacuum Einstein equations and requiring regularity at the origin. We will discuss the form of $H_{\mu\nu}$ in more in details in section 5.1. Equation (4.17) allows to perform a perturbative expansion of the action in both the tidal field amplitude and in the coupling constant $\kappa$. By performing this expansion, one schematically obtains, for the bulk action,

$$S_{\text{EH}}[g] + S_{\text{GF}}[\bar{g}, h] = -\frac{1}{2} \int d^4 x \left[ \partial_\mu h_{\alpha\beta} \partial^\mu h^{\alpha\beta} - \frac{1}{2} \partial_\mu h \partial^\mu h \right] \\ + \int d^4 x \left[ H(\partial h)^2 + h(\partial H) \partial h + \kappa h(\partial h)^2 \right] \\ + \int d^4 x \left[ H^2(\partial h)^2 + H h(\partial H) \partial h + \kappa H h(\partial h)^2 + \kappa^2 h^2(\partial h)^2 \right] + \dots, \quad (4.18)$$

where here and below indices are raised and lowered with $\eta_{\mu\nu}$. From the first term on the right-hand side we obtain the usual graviton propagator in de Donder gauge, while the

<footnote n="15">¹⁵Since $S_0$ is a constant term, we shall drop it in (4.14).</footnote>
higher-order terms in the second and third lines give the self-interaction vertices (which we have written here schematically without making index contractions explicit). Notice that, on the right-hand side, we have not explicitly written any terms containing only the external field $H_{\mu\nu}$, as they belong to $S_0[\bar{g}, x(\tau)]$ and are immaterial for the following calculations. Furthermore, since $H_{\mu\nu}$ satisfies the vacuum nonlinear Einstein equations, we do not have operators with a single $h_{\mu\nu}$. Additionally, for the choice of the gauge-fixing term (4.15), all interaction vertices containing at least one tidal field $H_{\mu\nu}$ receive a contribution from $S_{\text{GF}}$, while the self-interactions of $h_{\mu\nu}$ remain unaffected by $S_{\text{GF}}$.

Finally, we have to expand the matter part of the action, i.e. $S_{\text{pp}}$ and $S_{\text{T}}$. Rather than using the proper time, it is convenient to choose a frame in which the massive object is at rest at the origin, i.e. we choose the time $t$ such that

$$x^{\mu}(t) = (t, 0, 0, 0), \qquad v^{\mu}(t) \equiv \frac{dx^{\mu}}{dt} = (1, 0, 0, 0). \qquad (4.19)$$

The point-particle action then takes the form

$$S_{\text{pp}} = -M \int dt \frac{d\tau}{dt} = -M \int dt \sqrt{-g_{\mu\nu}v^{\mu}v^{\nu}}, \qquad (4.20)$$

and similarly for $S_{\text{T}}$.

In order to write the Feynmann rules, let us introduce the following conventions:

<figure><img src="image_2.png" /></figure>

Then, introducing the tensor

$$P_{\mu\nu\rho\sigma} = -\frac{1}{2}(\eta_{\mu\rho}\eta_{\nu\sigma} + \eta_{\mu\sigma}\eta_{\nu\rho} - \eta_{\mu\nu}\eta_{\rho\sigma}) \qquad (4.21)$$

to define the graviton propagator in de Donder gauge, the Feynman rules obtained from the bulk action (4.18) are

$$\mu\nu \xrightarrow{k} \rho\sigma = \frac{i}{k^2 - i0^+} P_{\mu\nu\rho\sigma}, \qquad (4.22)$$

$$k_1 \xrightarrow{k_2} \alpha_2 \beta_2 \qquad k_3 \xrightarrow{k_3} \alpha_3 \beta_3 = i\kappa\delta^{-(4)}(k_1 + k_2 + k_3)\mathcal{V}_{\alpha_1\beta_1\alpha_2\beta_2\alpha_3\beta_3}, \qquad (4.23)$$

$$k_2 \xrightarrow{k_3} \alpha_2 \beta_2 \qquad k_3 \xrightarrow{k_3} \alpha_3 \beta_3 = iH_{\alpha_1\beta_1}(-k_2 - k_3)V_3^{\alpha_1\beta_1}_{\alpha_2\beta_2\alpha_3\beta_3}, \qquad (4.24)$$

$$k_2 \xrightarrow{k_3} \alpha_2 \beta_3 \qquad k_3 \xrightarrow{k_3} \alpha_3 \beta_4 = i \int_q H_{\alpha_1 \beta_1}(q) H_{\alpha_2 \beta_2}(-k_2 - k_3 - q) V_4^{\alpha_1 \beta_1 \alpha_2 \beta_2}_{\alpha_3 \beta_3 \alpha_4 \beta_4}, \qquad (4.25)$$
where $\delta^{-(4)}(k) \equiv (2\pi)^4\delta^{(4)}(k)$. $V$ is the standard trilinear graviton vertex in de Donder gauge, while $V_3$ and $V_4$ are cubic and quartic vertices in the gauge (4.15). Note that we can ignore the $i0^+$ prescription in the propagator (4.22) because it is irrelevant for the following static computations. On the other hand, we will only need the following rule for the point-particle part of the action:<fnref n="16" />

$$ M \underset{\sim}{\stackrel{k}{\rightleftharpoons}}^{\mu\nu} = i\kappa \frac{M}{2} \int dt \, e^{-ik \cdot x(t)} v^{\mu} v^{\nu}. \qquad (4.26) $$

# 5 Matching with the point particle

In this section we will compute the linear and quadratic response to an external tidal field in the EFT approach outlined above and match it to the full GR calculation of section 3.2. First, we will derive the tidal field $H_{\mu\nu}$ at both first and second order, expressing it in a form that is suitable for the calculation of Feynman diagrams in the EFT approach.

## 5.1 The external field

As mentioned above, by construction $H_{\mu\nu}$ satisfies the vacuum Einstein equations and is regular at the location of the point particle, which we place for convenience at the origin of the coordinate system. We shall parameterize it as described in section 2.1, and split it into even and odd components, $H_{\mu\nu}(x) = H_{\mu\nu}^{+}(x) + H_{\mu\nu}^{-}(x)$. We can then apply the decomposition from eqs. (2.4) and (2.5). Specifically, we will denote the components of $H_{\mu\nu}^{\pm}$ with a bar and write

$$ H_{\mu\nu}^{+} = \begin{pmatrix} f(r)\bar{H}_0 & \bar{H}_1 & \nabla_A \bar{H}_0 \\ * & f^{-1}(r)\bar{H}_2 & \nabla_A \bar{H}_1 \\ * & * & r^2 \bar{K} \gamma_{AB} + r^2 (\nabla_A \nabla_B - \frac{1}{2} \gamma_{AB}) \bar{G} \end{pmatrix}, \qquad (5.1) $$

$$ H_{\mu\nu}^{-} = \begin{pmatrix} 0 & 0 & -\epsilon_A C \nabla_C \bar{h}_0 \\ * & 0 & -\epsilon_A C \nabla_C \bar{h}_1 \\ * & * & -\frac{1}{2}(\epsilon_A C \nabla_C \nabla_B + \epsilon_B C \nabla_C \nabla_A) \bar{h}_2 \end{pmatrix}. \qquad (5.2) $$

Moreover, since the gauge fixing in eq. (4.15) ensures that the final result is invariant under diffeomorphisms of $\bar{g}_{\mu\nu} = \eta_{\mu\nu} + H_{\mu\nu}$, we can choose $H_{\mu\nu}$ in any gauge we prefer. To be closer to the GR calculation, we will work in the RW gauge (see eq. (2.9)). However, since all calculations using Feynman diagrams are performed in Cartesian coordinates, it is useful to first rewrite (5.1) and (5.2) in a Lorentz-covariant form, which will facilitate their conversion to Cartesian coordinates when required. To this end, it is convenient to define the spatial radial vector

$$ r^{\mu} = \mathbb{P}_{\nu}^{\mu} x^{\nu} \equiv (\delta^{\mu}_{\nu} + v^{\mu} v_{\nu}) x^{\nu}, \qquad (5.3) $$

and the unit vector $\hat{r}^{\mu} = r^{\mu}/r$, $\hat{r}^{\mu}\hat{r}_{\mu} = 1$, where we remind the reader that in this section indices are raised and lowered using the flat metric. As derived in appendix C, eqs. (5.1)

<footnote n="16">¹⁶There are no direct couplings with $H_{\mu\nu}$ coming from eq. (4.20) because $H_{\mu\nu}$ vanishes on the worldline.</footnote>
and (5.2) can be written in Lorentz-covariant form as

$$H_{\mu\nu}^{+} = (\bar{H}_{0} + \bar{K})v_{\mu}v_{\nu} + 2\bar{H}_{1}v_{(\mu}\hat{r}_{\nu)} + (\bar{H}_{2} - \bar{K})\hat{r}_{\mu}\hat{r}_{\nu} + \bar{K}\eta_{\mu\nu}, \quad (5.4)$$

$$H_{\mu\nu}^{-} = -2v^{\alpha}\varepsilon_{\alpha\rho\sigma(\mu}v_{\nu)}\hat{r}^{\rho}\partial^{\sigma}\bar{h}_{0} + 2v^{\alpha}\varepsilon_{\alpha\rho\sigma(\mu}\hat{r}_{\nu)}\hat{r}^{\rho}\partial^{\sigma}\bar{h}_{1}, \quad (5.5)$$

where $\varepsilon_{\alpha\rho\sigma\mu}$ is the flat-space Levi-Civita symbol ($\varepsilon_{0123} = +1$).

Next, we insert these expressions into the vacuum Einstein equations, expand them perturbatively in the amplitude of the tidal field, and solve for $H_{\mu\nu}$ order by order, ensuring regularity at the origin. To take advantage of the spherical symmetry of the background when solving the Einstein equations, it is convenient to decompose into spherical harmonics, as in section 2.1. In order to work in Cartesian coordinates, we shall then express the spherical harmonics in terms of constant symmetric trace-free (STF) tensors $\mathcal{Y}_{\alpha_1...\alpha_\ell}^m$ (see e.g., [89, 90]), as

$$Y_{\ell}^{m}(\theta, \phi) = \mathcal{Y}_{\alpha_1...\alpha_\ell}^{m} \hat{r}^{\alpha_1} \dots \hat{r}^{\alpha_\ell}. \quad (5.6)$$

Here the STF tensors $\mathcal{Y}_{\alpha_1...\alpha_\ell}^m$ have covariant indices, but they are purely spatial, i.e., $\mathbb{P}_{\beta}^{\alpha_i} \mathcal{Y}_{\alpha_1...\alpha_i...\alpha_\ell}^m = \mathcal{Y}_{\alpha_1...\beta...\alpha_\ell}^m$, $\forall i \in \{1, \dots, \ell\}$. This allows us to rewrite $H_{\mu\nu}^{\pm}$ in Cartesian coordinates as

$$H_{\mu\nu}^{\pm} = \sum_{\ell} \mathcal{H}_{\mu\nu|\alpha_1 \dots \alpha_\ell}^{\pm} x^{\alpha_1} \dots x^{\alpha_\ell}, \quad (5.7)$$

where we introduced the tidal tensors $\mathcal{H}_{\mu\nu|\alpha_1...\alpha_\ell}^{\pm}$, which are purely spatial constant tensors.

To illustrate the procedure concretely, let us expand the tidal tensor perturbatively, as $\mathcal{H}_{\mu\nu|\alpha_1...\alpha_\ell}^{\pm} = (1)\mathcal{H}_{\mu\nu|\alpha_1...\alpha_k}^{\pm} + (2)\mathcal{H}_{\mu\nu|\alpha_1...\alpha_\ell}^{\pm} + \dots$, and derive the quadratic tidal tensor sourced by two linear quadrupoles. This is what is required to match the GR calculation performed in section 3. Solving the linear vacuum Einstein equations for $\ell = 2$ and requiring regularity at the origin, we find the even and odd linear tidal fields in RW gauge,

$$(1) H_{\mu\nu}^{\pm}(x) = (1) \mathcal{H}_{\mu\nu|\alpha\beta}^{\pm} x^{\alpha} x^{\beta}, \quad (5.8)$$

with

$$(1) \mathcal{H}_{\mu\nu|\alpha\beta}^{+} = \sum_{m} \frac{\mathcal{E}_{2m}^{+}}{r_s^2} (\eta_{\mu\nu} + 2v_{\mu}v_{\nu}) \mathcal{Y}_{\alpha\beta}^{m}, \quad (5.9a)$$

$$(1) \mathcal{H}_{\mu\nu|\alpha\beta}^{-} = -4 \sum_{m} \frac{\mathcal{E}_{2m}^{-}}{r_s^2} v^{\tau} v_{(\mu} \varepsilon_{\nu)} \alpha_{\tau}{}^{\sigma} \mathcal{Y}_{\sigma\beta}^{m}, \quad (5.9b)$$

where $\mathcal{E}_{\ell m}^{\pm}$ are the amplitudes of the even/odd tidal fields, chosen so to match the $r_s/r \ll 1$ limit of the full GR tidal field solution, see eqs. (3.23) and (3.24).

We now use these linear solutions to solve the vacuum Einstein equations at second order. As discussed in section 3.3, given a $\ell = 2$ mode at linear order, the quadratic solution is given by a superposition of $\ell = 0, 2, 4$ multipoles. Solving the vacuum Einstein equations at second order and imposing regularity, we find

$$(2) H_{\mu\nu}^{\pm} = (2) \mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{\pm} x^{\alpha_1} x^{\alpha_2} x^{\alpha_3} x^{\alpha_4}, \quad (5.10)$$

where the tidal tensor can be divided into monopole, quadrupole and hexadecapole parts, i.e.,

$$(2) \mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{\pm} = (2) \mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{\pm, \ell=0} + (2) \mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{\pm, \ell=2} + (2) \mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{\pm, \ell=4}, \quad (5.11)$$
Their explicit forms are given by

$$^{(2)}\mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{+,-0} = \frac{1}{2\sqrt{\pi}} \sum_{m_1,m_2} I_{022}^{0m_1m_2} \sum_{p=\pm} \frac{\mathcal{E}_{2m_1}^p \mathcal{E}_{2m_2}^p}{r_s^4} \left[ (c_{\eta}^{(0,p)}\eta_{\mu\nu} + c_v^{(0,p)}v_{\mu}v_{\nu}) \mathbb{P}_{\alpha_1\alpha_2} + c_r^{(0,p)}\mathbb{P}_{\alpha_1\mu}\mathbb{P}_{\alpha_2\nu} \right] \mathbb{P}_{\alpha_3\alpha_4}, \qquad (5.12a)$$

$$^{(2)}\mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{+,-2} = \sum_{m_1,m_2} I_{222}^{(m_1+m_2)m_1m_2} \sum_{p=\pm} \frac{\mathcal{E}_{2m_1}^p \mathcal{E}_{2m_2}^p}{r_s^4} \left[ (c_{\eta}^{(2,p)}\eta_{\mu\nu} + c_v^{(2,p)}v_{\mu}v_{\nu}) \mathbb{P}_{\alpha_1\alpha_2} + c_r^{(2,p)}\mathbb{P}_{\alpha_1\mu}\mathbb{P}_{\alpha_2\nu} \right] \mathcal{Y}_{\alpha_3\alpha_4}^{m_1+m_2}, \qquad (5.12b)$$

$$^{(2)}\mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{-,\ell=2} = -4 \sum_{m_1,m_2} I_{222}^{(m_1+m_2)m_1m_2} \frac{\mathcal{E}_{2m_1}^{+} \mathcal{E}_{2m_2}^{-}}{r_s^4} v^{\alpha} v_{(\mu\varepsilon_{\nu})\alpha\sigma\alpha_3} d^{(2)}\mathbb{P}_{\alpha_1\alpha_2} \mathcal{Y}_{\sigma\alpha_4}^{m_1+m_2}, \qquad (5.12c)$$

$$^{(2)}\mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{+,-4} = \sum_{m_1,m_2} I_{422}^{(m_1+m_2)m_1m_2} \sum_{p=\pm} \frac{\mathcal{E}_{2m_1}^p \mathcal{E}_{2m_2}^p}{r_s^4} [c_{\eta}^{(4,p)}\eta_{\mu\nu} + c_v^{(4,p)}v_{\mu}v_{\nu}] \mathcal{Y}_{\alpha_1\alpha_2\alpha_3\alpha_4}^{m_1+m_2}, \qquad (5.12d)$$

$$^{(2)}\mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{-,\ell=4} = -8 \sum_{m_1,m_2} I_{422}^{(m_1+m_2)m_1m_2} \frac{\mathcal{E}_{2m_1}^{+} \mathcal{E}_{2m_2}^{-}}{r_s^4} v^{\alpha} v_{(\mu\varepsilon_{\nu})\alpha\sigma\alpha_1} d^{(4)}\mathcal{Y}_{\sigma\alpha_2\alpha_3\alpha_4}^{m_1+m_2}. \qquad (5.12e)$$

The numerical coefficients $c^{(\ell,p)}$ and $d^{(\ell)}$ appearing in these expressions are given in appendix B.

In order to perform the computation in the next sections, we need the tidal field in Fourier space. From eq. (5.7), the coordinate dependence of the tidal field is given by the products $x^{\alpha_1} \cdots x^{\alpha_\ell}$, which are contracted with spacelike tensors. Therefore, all we need to know is the Fourier transform of $r^{\alpha_1} \cdots r^{\alpha_\ell}$, which reads<fnref n="17" />

$$\int d^4 x \, e^{-iq \cdot x} r^{\mu_1} \cdots r^{\mu_\ell} = (-i)^\ell \delta^{-(4)}(q) \mathbb{P}^{\mu_1 \nu_1} \cdots \mathbb{P}^{\mu_\ell \nu_\ell} \frac{\partial^\ell}{\partial q^{\nu_1} \cdots \partial q^{\nu_\ell}}. \qquad (5.14)$$

Introducing the shorthand notation

$$\partial_{\mu_1 \dots \mu_\ell}^{(q)} \equiv (-i)^\ell \frac{\partial^\ell}{\partial q^{\nu_1} \dots \partial q^{\nu_\ell}}, \qquad (5.15)$$

the tidal field in Fourier space at linear and quadratic order is given by

$$^{(1)}H_{\mu\nu}^{\pm}(q) = \delta^{-(4)}(q) \mathcal{H}_{\mu\nu|\alpha_1\alpha_2}^{\pm \alpha_1\alpha_2} \partial_{\alpha_1\alpha_2}^{(q)}, \qquad (5.16)$$

$$^{(2)}H_{\mu\nu}^{\pm}(q) = \delta^{-(4)}(q) \mathcal{H}_{\mu\nu|\alpha_1\alpha_2\alpha_3\alpha_4}^{\pm \alpha_1\alpha_2\alpha_3\alpha_4} \partial_{\alpha_1\alpha_2\alpha_3\alpha_4}^{(q)}, \qquad (5.17)$$

$$\int d^4 x \, e^{-iq \cdot x} x^{j_1} \cdots x^{j_\ell} = \delta^{-}(q^0) i^{\ell} \frac{\partial^{\ell}}{\partial q_{j_1} \cdots \partial q_{j_\ell}} \delta^{-(3)}(q). \qquad (5.13)$$

When the right-hand side acts on a function of $q^\mu$, one gets

$$\int_q \delta^{-}(q^0) i^{\ell} \frac{\partial^{\ell}}{\partial q_{i_1} \cdots \partial q_{i_{\ell}}} \delta^{-}(q) f(q) = (-i)^{\ell} \frac{\partial^{\ell} f(q)}{\partial q_{i_1} \cdots \partial q_{i_{\ell}}} \Big|_{q^{\mu}=0}.$$

Therefore, using this relation and upon covariantization of eq. (5.13) one finds eq. (5.14).

<footnote n="17">¹⁷To derive this relation, note that, in the coordinate system where $r^\alpha = (0, x^i)$, we have</footnote>
where $^{(1)}\mathcal{H}_{\mu\nu|}^{\pm \alpha_1\alpha_2} \equiv \eta^{\alpha_1\beta_1}\eta^{\alpha_1\beta_1} {}^{(1)}\mathcal{H}_{\mu\nu|\beta_1\beta_2}^{\pm}$ (see eq. (5.9)) and $^{(2)}\mathcal{H}_{\mu\nu|}^{\pm \alpha_1\dots\alpha_4} \equiv \eta^{\alpha_1\beta_1}\dots\eta^{\alpha_2\beta_2} {}^{(2)}\mathcal{H}_{\mu\nu|\beta_1\dots\beta_4}^{\pm}$ (see eqs. (5.11) and (5.12)).

## 5.2 Quadratic response

We are now ready to compute the one-point function due to the interaction between the external tidal field and the point particle, $h_{\mu\nu}(x) \equiv \langle h_{\mu\nu}(x) \rangle$. This can be expanded as

$$h_{\mu\nu} = \sum_{n=0}^{\infty} {}^{(n)}h_{\mu\nu}, \qquad (5.18)$$

where ${}^{(n)}h_{\mu\nu}$ is the response at n-th-order in the tidal field amplitude $\mathcal{E}$. Diagrammatically, ${}^{(n)}h_{\mu\nu}$ can be represented as

$$ {}^{(n)}h_{\mu\nu}(x) = M \bullet \mathcal{A} \underset{\substack{\varepsilon \xrightarrow{\hspace{1.5em}} \hspace{1.5em} \varepsilon \\ \mathbf{n}}} \overset{x}{\underset{\varepsilon \xrightarrow{\hspace{1.5em}} \hspace{1.5em} \mathbf{n}}}{\mu\nu} + \lambda \bullet \mathcal{A} \underset{\substack{\varepsilon \xrightarrow{\hspace{1.5em}} \hspace{1.5em} \varepsilon \\ \mathbf{n}}} \overset{x}{\underset{\varepsilon \xrightarrow{\hspace{1.5em}} \hspace{1.5em} \mathbf{n}}}{\mu\nu}, \qquad (5.19) $$

where the blob stands for all possible interaction vertices between $H_{\mu\nu}$ and $h_{\mu\nu}$, at a specific order in perturbation theory. Thus, this diagram distinguishes bewteen two types of nonlinearities. Those that arise from higher-order vertices in the bulk action (see eqs. (4.24) and (4.25)), corresponding to $\mathcal{O}(r_s)$ corrections to the external tidal field, which are represented by the first diagram on the right-hand side. Additionally, there are nonlinearities that arise as a genuine vertex response of the worldline, receiving $r_s$-corrections from bulk interactions. They are represented by the second diagram on the right-hand side.

In the following, we will show that the first diagram on the right-hand side alone reproduces the full GR calculation at leading order in $r_s/r$, implying that the response nonlinearities must be absent, i.e., $\lambda = 0$. The matching can be extended to higher orders in $r_s$ by considering diagrams with multiple mass insertions on the worldline and, eventually, reconstructing the full tidal field profile computed in GR in section 3. However, these corrections are not needed to derive our conclusions; see e.g. [86] for an extensive discussion in the linear case. Nonetheless, as a sanity check of the validity of our method at higher-order in $r_s$, we have verified that in the simplest case of $n = 0$, i.e., in the absence of the tidal field, multiple mass insertions reconstruct the Schwarzschild metric perturbatively up to $\mathcal{O}(r_s^2)$. The calculation is reported in details in appendix E.

We now perform the computation of $h_{\mu\nu}(x)$ up to linear order in $r_s$ and up to quadratic order in the amplitude of the tidal field $\mathcal{E}$, i.e.,

$$ h_{\mu\nu}(x) = {}^{(0)}h_{\mu\nu}(x) + {}^{(1)}h_{\mu\nu}(x) + {}^{(2)}h_{\mu\nu}(x) + \mathcal{O}(\mathcal{E}^3), \qquad (5.20) $$

where ${}^{(0)}h_{\mu\nu}(x)$ is the Schwarzschild metric expanded at $\mathcal{O}(r_s)$; see in appendix E. The diagrams to compute are represented in figure 1. The only diagram for the computation of ${}^{(1)}h_{\mu\nu}(x)$ is figure 1(a), while to compute ${}^{(2)}h_{\mu\nu}(x)$ we need 1(b,c,d).

Taking the linear tidal field to be an $\ell = 2$ mode, we have

$$ {}^{(1)}h_{\mu\nu}^{\pm}(x) = \frac{\kappa M}{2} \int_k \frac{\delta^{-}(k^0) e^{ik \cdot x}}{k^2} \int_q \delta^{-}(q) {}^{(1)}\mathcal{H}_{\rho\sigma, \alpha_1\alpha_2}^{\pm \alpha_1\alpha_2} \partial_{\alpha_1\alpha_2}^{(q)} \left( \frac{\mathcal{A}_{\mu\nu}^{(1a)}\rho\sigma(k, q)}{(k-q)^2} \right), \qquad (5.21) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Feynman diagrams needed for the computation of $h_{\mu\nu}$. Diagram (a) yields the order-$r_s$ correction to the linear tidal field solution. Diagrams (b), (c) and (d) represent instead order-$r_s$ corrections to the tidal source at second order in the external field amplitude.</figcaption></figure>

where

$$A_{\mu\nu}^{(1a)\rho\sigma}(k, q) \equiv v^{\beta_1} v^{\delta_1} P_{\beta_1\delta_1\beta_2\delta_2} V^{\rho\sigma\beta_2\delta_2\beta_3\delta_3} P_{\beta_3\delta_3\mu\nu}. \quad (5.22)$$

Performing the contraction and integrations using the explicit form of the tidal tensors in eq. (5.9), we eventually obtain the even and odd part of the linear one-point function, respectively,

$$\kappa^{(1)} h_{\mu\nu}^{+}(x) = \frac{r}{r_s} \sum_{m} \mathcal{E}_{2m}^{+} \left( 2\mathcal{Y}_{\mu\nu}^{m} + \eta_{\mu\nu} \mathcal{Y}_{\rho\sigma}^{m} \hat{r}^{\rho} \hat{r}^{\sigma} \right), \quad (5.23a)$$

$$\kappa^{(1)} h_{\mu\nu}^{-}(x) = -2 \frac{r}{r_s} \sum_{m} \mathcal{E}_{2m}^{-} v^{\alpha} v_{(\mu} \varepsilon_{\nu)} \rho_{\alpha} \sigma \mathcal{Y}_{\sigma\tau}^{m} \hat{r}^{\rho} \hat{r}^{\tau}. \quad (5.23b)$$

We can now compute $^{(2)}h_{\mu\nu}(x)$ from the last three diagrams in figure 1. They can be recasted as

$$^{(2)}h_{\mu\nu}(x)_{1b} = \frac{\kappa M}{2} \int_k \frac{\delta^{-}(k^0) e^{ik \cdot x}}{k^2} \int_q \delta^{-}(4)(q) \sum_{p=\pm} {}^{(2)}\mathcal{H}_{\rho_1 \sigma_1}^{p} \left( \frac{\alpha_1 \alpha_2 \alpha_3 \alpha_4}{\delta^{(q)}_{\alpha_1 \alpha_2 \alpha_3 \alpha_4}} \right) \left( \frac{A_{\mu\nu}^{(1a)} \rho_1 \sigma_1 (k,q)}{(k-q)^2} \right), \quad (5.24)$$

$$^{(2)}h_{\mu\nu}(x)_{1x} = \frac{\kappa M}{2} \int_k \frac{\delta^{-}(k^0) e^{ik \cdot x}}{k^2} \int_{q_1, q_2} \delta^{-}(4)(q_1) \delta^{-}(4)(q_2) \sum_{p_1=\pm} \sum_{p_2=\pm} {}^{(1)}\mathcal{H}_{\rho_1 \sigma_1}^{p_1} \left( \frac{\alpha_1 \alpha_2}{\delta^{(q_1)}_{\alpha_1 \alpha_2}} \right) {}^{(1)}\mathcal{H}_{\rho_2 \sigma_2}^{p_2} \left( \frac{\alpha_3 \alpha_4}{\delta^{(q_2)}_{\alpha_3 \alpha_4}} \right) \times \left( \frac{A_{\mu\nu}^{(2x)} \rho_1 \sigma_1 \rho_2 \sigma_2 (k, q_1, q_2)}{(k-q_{12})^2} \right), \quad x \in \{c, d\}, \quad (5.25)$$

where $q_{12} \equiv q_1 + q_2$. Clearly, figure 1(b) has the same form of figure 1(a) with the linear tidal field replaced by the quadratic solution. On the other hand, figures 1(c) and (d) are new topologies arising at quadratic order, with

$$A_{\mu\nu}^{(1c)\rho_1\sigma_1\rho_2\sigma_2}(k, q_1, q_2) \equiv -v^{\beta_1} v^{\delta_1} P_{\beta_1\delta_1\beta_2\delta_2} V_3^{\rho_1\sigma_1\beta_2\delta_2\beta_3\delta_3} \frac{P_{\beta_3\delta_3\beta_4\delta_4}}{(k-q_2)^2} V_3^{\rho_2\sigma_2\beta_4\delta_4\beta_5\delta_5} P_{\beta_5\delta_5\mu\nu}, \quad (5.26)$$

$$A_{\mu\nu}^{(1d)\rho_1\sigma_1\rho_2\sigma_2}(k, q_1, q_2) \equiv \frac{v^{\beta_1} v^{\delta_1}}{2} P_{\beta_1\delta_1\beta_2\delta_2} V_4^{\rho_1\sigma_1\rho_2\sigma_2\beta_2\delta_2\beta_3\delta_3} P_{\beta_3\delta_3\mu\nu}. \quad (5.27)$$

Notice that there is an extra factor of 1/2 in the last amplitude from the symmetry factor of the diagram.

We split once again the quadratic one-point function into even and odd contributions, $^{(2)}h_{\mu\nu}(x) = {}^{(2)}h_{\mu\nu}^{+}(x) + {}^{(2)}h_{\mu\nu}^{-}(x)$. The even contribution of eq. (5.24) comes from inserting the
quadratic even tidal tensor, eqs. (5.10) and (5.12), while the even contribution of eq. (5.25) comes from replacing the tidal tensor by either two even linear modes $^{(1)}\mathcal{H}_{\rho\sigma|}^{+\alpha\beta}$ or two linear odd modes $^{(1)}\mathcal{H}_{\rho\sigma|}^{-\alpha\beta}$, see eq. (5.9). Therefore, for each diagram, the even one-point function is the sum of a part proportional to $(\mathcal{E}^{+})^2$ (the even × even “response”), denoted by $^{(2)}h_{\mu\nu}^{(+)}(x)_{1x}$, and a part proportional to $(\mathcal{E}^{-})^2$ (the odd × odd “response”), denoted by $^{(2)}h_{\mu\nu}^{(-)}(x)_{1x}$. Here we report explicitly only the tt component, as this is the one relevant for the matching. After using the identities involving the STF tensors $\mathcal{Y}_{ij}^m$ in appendix D and after summing over the three diagrams, one finds

$$ \begin{aligned} \kappa^{(2)} h_{tt}^{+}(x) = \left(\frac{r}{r_s}\right)^3 \sum_{m_1, m_2} \left[ \begin{aligned} & \mathcal{E}_{2m_1}^{+} \mathcal{E}_{2m_2}^{+} \left( \frac{7}{4\sqrt{\pi}} I_{022}^{0m_1m_2} + \frac{7}{4} I_{222}^{(m_1+m_2)m_1m_2} \mathcal{Y}_{ij}^{m_1+m_2} \hat{r}^i \hat{r}^j \right. \\ & \qquad \left. - \frac{35}{18} I_{422}^{(m_1+m_2)m_1m_2} \mathcal{Y}_{ijkl}^{m_1+m_2} \hat{r}^i \hat{r}^j \hat{r}^k \hat{r}^l \right) \\ & + \mathcal{E}_{2m_1}^{-} \mathcal{E}_{2m_2}^{-} \left( \frac{3}{\sqrt{\pi}} I_{022}^{0m_1m_2} - \frac{3}{4} I_{222}^{(m_1+m_2)m_1m_2} \mathcal{Y}_{ij}^{m_1+m_2} \hat{r}^i \hat{r}^j \right. \\ & \qquad \left. - \frac{5}{9} I_{422}^{(m_1+m_2)m_1m_2} \mathcal{Y}_{ijkl}^{m_1+m_2} \hat{r}^i \hat{r}^j \hat{r}^k \hat{r}^l \right) \end{aligned} \right]. \end{aligned} \quad (5.28) $$

The computation of the odd modes of $h_{\mu\nu}$ is analogous. We use the odd quadratic tidal tensor for the diagram in figure 1(b) and the product of one even and one odd linear tidal tensors for the diagrams in figure 1(c) and (d). Therefore, the odd contribution to the one-point function is always proportional to $\mathcal{E}^{+}\mathcal{E}^{-}$. We report the components $ti$, that are the ones relevant for the matching. They read

$$ \kappa^{(2)} h_{ti}^{-}(x) = - \left(\frac{r}{r_s}\right)^3 \sum_{m_1, m_2} \mathcal{E}_{2m_1}^{+} \mathcal{E}_{2m_2}^{-} \epsilon_{ij}^{k} \left( 5^0 I_{222}^{(m_1+m_2)m_1m_2} \mathcal{Y}_{kl}^{m_1+m_2} - \frac{13}{5} I_{422}^{(m_1+m_2)m_1m_2} \mathcal{Y}_{klpq}^{m_1+m_2} \hat{r}^p \hat{r}^q \right) \hat{r}^j \hat{r}^l. \quad (5.29) $$

Before proceeding to do the matching, as a consistency check we have verified that the obtained $h_{\mu\nu}(x)$ satisfies the gauge condition

$$ \bar{g}^{\alpha\beta} \left( \bar{\nabla}_{\alpha} h_{\beta\mu} - \frac{1}{2} \bar{\nabla}_{\mu} h_{\alpha\beta} \right) = 0, \quad (5.30) $$

where $\bar{\nabla}_{\alpha}\bar{g}_{\mu\nu} = 0$, imposed by the gauge-fixing term (4.15) in the action. Inserting $\bar{g}_{\alpha\beta} = \eta_{\alpha\beta} + H_{\alpha\beta}$ in this equation and working at linear order in the tidal field amplitude, this equation becomes

$$ \left( \eta^{\alpha\beta} \partial_{\alpha} {}^{(1)}h_{\beta\mu} - \frac{1}{2} \partial_{\mu} {}^{(1)}h \right) - {}^{(1)}H^{\alpha\beta} \left( \partial_{\alpha} {}^{(0)}h_{\beta\mu} - \frac{1}{2} \partial_{\mu} {}^{(0)}h_{\alpha\beta} \right) + {}^{(0)}h_{\mu\alpha} \left( \partial_{\beta} {}^{(1)}H^{\beta\alpha} - \frac{1}{2} \partial^{\alpha} {}^{(1)}H \right) = \mathcal{O}(\mathcal{E}^2). \quad (5.31) $$

Since it is linear, this identity must be satisfied independently for the even and odd perturbations. Indeed, replacing in this equation the linear even solution, eq. (5.23a), the first term becomes

$$ \eta^{\alpha\beta} \partial_{\alpha} {}^{(1)}h_{\beta\mu}^{+} - \frac{1}{2} \partial_{\mu} {}^{(1)}h^{+} = \frac{\kappa}{16\pi r_s^2} \sum_m \mathcal{E}_{2m}^{+} \mathcal{Y}_{\alpha\beta}^{m} \hat{r}^{\alpha} \hat{r}^{\beta} \hat{r}_{\mu}, \quad (5.32) $$
which cancels with the second term, while the last term vanishes. For the odd solution, eq. (5.23b), each of the terms in this equation vanishes individually at $O(\mathcal{E})$, so that the above identity is trivially satisfied. Proceeding analogously at second-order, we have verified that also at that order the one-point function $^{(2)}h_{\mu\nu}(x)$ satisfies eq. (5.30), which is a non-trivial check for both even and odd solutions.

## 5.3 Matching

Here we want to compare the one-point function computed in the previous subsection with the solution derived in section (3.3). However, the former has been derived in the gauge that satisfies eq. (5.30) while the latter has been derived in RW gauge. To compare the two, we find it convenient to rewrite the EFT result in RW gauge. In the context of the EFT calculation, the field $h_{\mu\nu}$ is treated as a spin-2 field propagating on the curved background $\bar{g}_{\mu\nu} = \eta_{\mu\nu} + H_{\mu\nu}$. The linear gauge transformation $x^{\mu} \rightarrow x^{\mu} + \xi^{\mu}$ that converts $h_{\mu\nu}$ to the RW gauge is given by (see, e.g., [62])

$$ \kappa h_{\mu\nu}^{\text{RW}} = \kappa h_{\mu\nu} - \xi^{\rho} \partial_{\rho} \bar{g}_{\mu\nu} - 2 \bar{g}_{\rho(\mu} \partial_{\nu)} \xi^{\rho} + \mathcal{O}(\xi^2), \quad (5.33) $$

where we have included only terms up to linear order in $\xi$, as higher-order terms contribute only at higher order the $r_s$ expansion.<fnref n="18" />

For the even sector, we match the $tt$ component of the metric, since the other components are determined by it through Einstein's equations. Additionally, at this order in perturbation theory, the $ti$ components are parity-odd, allowing us to use them to match the odd sector. In the static limit, the transformations for the $tt$ and $ti$ components of the metric simplify as

$$ \kappa h_{tt}^{\text{RW}+} = \kappa h_{tt}^{+} - \xi^{k} \partial_{k} H_{tt}^{+}, \quad (5.34) $$

$$ \kappa h_{ti}^{\text{RW}-} = \kappa h_{ti}^{-} - \xi^{k} \partial_{k} H_{ti}^{-} - H_{kt}^{-} \partial_{i} \xi^{k}, \quad (5.35) $$

where to write the right-hand sides we have used that $H_{ti}^{+} = 0 = H_{tt}^{-}$, and that $\xi_t$ is purely even. Moreover, note that $\xi^k$ appearing in these two equations must be even because $\xi^k \partial_k H_{tt}^+$ and $\xi^k \partial_k H_{ti}^-$ must be respectively even and odd. For the gauge transformation at linear order in $\mathcal{E}$ we need $\xi^\mu = {}^{(0)}\xi^\mu$, where ${}^{(0)}\xi^\mu = \frac{r_s}{2}\hat{r}^\mu$ is the coordinate transformation of the Schwarzschild metric from de Donder to RW gauge. For the quadratic order, we insert $\xi^\mu = {}^{(0)}\xi^\mu + {}^{(1)}\xi^\mu$, where

$$ {}^{(1)}\xi^{\mu} = \sum_{m} \mathcal{E}_{2m}^{+} \frac{r^2}{2r_s} \mathcal{Y}_{\alpha\beta}^{m} (2\eta^{\beta\mu} - \hat{r}^{\mu} \hat{r}^{\beta}) \hat{r}^{\alpha} \quad (5.36) $$

is found by imposing that ${}^{(2)}h_{\mu\nu}^{\text{RW}+}$ is diagonal as in the RW gauge.

Finally, since $H_{\mu\nu}$ is already in RW gauge, we can compare the metric perturbation obtained with the EFT calculation after transforming to RW gauge, i.e., $\delta g_{\mu\nu}^{\text{EFT}} = H_{\mu\nu} +$

<footnote n="18">¹⁸We need to perform a second-order gauge transformation in $\mathcal{E}$ but at first order in $r_s$. Since the gauge parameter $\xi^\mu$ is at least of order $r_s$, any second-order term in eq. (5.33) contributes only to higher $r_s$-corrections that we do not compute here. Note that ${}^{(1)}\xi$ is one $r_s$-order higher than the tidal field, as seen from eqs. (5.9) and (5.36).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Feynman diagrams responsible for the (a) linear ad (b) nonlinear tidal deformation. In principle, one could have terms proportional to $\lambda_2^{(E,B)}$ in the topology represented in figure (b). These are excluded because the matching at leading order in the tidal field implies $\lambda_2^{(E,B)} = 0$. See also section 6.2.</figcaption></figure>

$\kappa h_{\mu\nu}^{\text{RW}}$,<fnref n="19" /> with the exact metric perturbation obtained in the full GR calculation, converted in Cartesian coordinates. The latter reads

$$ \delta g_{tt}^{+}(\boldsymbol{x}) = \sum_{m} f(r) \left[ {}^{(2)}H_{0}^{0m}(r) + {}^{(1)}H_{0}^{2m}(r) + {}^{(2)}H_{0}^{2m}(r) \right] \mathcal{Y}_{\mu\nu}^{m} \hat{r}^{\mu} \hat{r}^{\nu} + {}^{(2)}H_{0}^{4m}(r) \mathcal{Y}_{\mu\nu\alpha\beta}^{m} \hat{r}^{\mu} \hat{r}^{\nu} \hat{r}^{\alpha} \hat{r}^{\beta} \right], \quad (5.37) $$

$$ \delta g_{ti}^{-}(\boldsymbol{x}) = - \sum_{m} \varepsilon_{ij}{}^{k} \hat{r}^{j} \partial_{k} \left[ \left( {}^{(1)}h_{0}^{2m}(r) + {}^{(2)}h_{0}^{2m}(r) \right) \mathcal{Y}_{\mu\nu}^{m} \hat{r}^{\mu} \hat{r}^{\nu} + {}^{(2)}h_{0}^{4m}(r) \mathcal{Y}_{\mu\nu\alpha\beta}^{m} \hat{r}^{\mu} \hat{r}^{\nu} \hat{r}^{\alpha} \hat{r}^{\beta} \right], \quad (5.38) $$

respectively for the $tt$ and $ti$ components. We remind that the GR solutions are given in section 3.3, see eqs. (3.23) and (3.24) for the linear solutions and eqs. (3.25), (3.28), (3.29), (3.32) and (3.33) for the quadratic solutions.

As anticipated, we find that $\delta g_{\mu\nu}^{\text{EFT}}$ matches the full GR solution above at leading and next-to-leading order in $r_s/r$, without the inclusion of any of the tidal operators in the action (4.10). This implies that the contributions coming from the diagrams depicted in figure 2 must vanish. As a result,

$$ \lambda_2^{(E)} = \lambda_2^{(B)} = 0, \quad \lambda_{222}^{(E3)} = \lambda_{222}^{(EB2)} = 0. \quad (5.39) $$

We will extend this result to all orders in derivatives in the next section.

## 6 Higher multipoles

In section 3.1, we computed the static solutions for the metric perturbations induced at second order by two quadrupolar static fields. In particular, we have checked explicitly that the solutions for $H_0$ and $h_0$ (see eqs. (3.25) and (3.28) for the induced quadrupole, and eqs. (3.29)–(3.32) for the induced hexadecapole) take the form of simple finite polynomials with only positive powers of $r/r_s$. This was crucial in the matching performed in section 5 to conclude that the quadratic Love number couplings vanish at lowest order in the number of spatial derivatives in the worldline EFT (i.e., eq. (5.39)).

Our goal in this section is to show that this conclusion remains true to all orders in the derivative expansion. First of all, we will show in full generality that, for generic multipoles,

<footnote n="19">¹⁹We match only perturbations linear or quadratic in $\mathcal{E}$, since the lowest order matches trivially with the Schwarzschild metric.</footnote>
the static solutions induced at second order in full GR do not contain any decaying falloff of the form $1/r^n$, with $n > 1$, when expanded at large distances from the source. Then, in section 6.2, will use this result to perform the matching at generic order in the EFT and show that all quadratic Love number couplings vanish for Schwarzschild black holes in general relativity.

## 6.1 General polynomial form of second-order static solutions

We start by discussing the full GR solutions for generic multipoles. For the sake of the presentation, we will consider the even and odd induced solutions separately.

### 6.1.1 Even sector

Let us start by considering eq. (2.34) for the second-order even mode $^{(2)}H_0^{\ell m}(r)$. The source is given in eq. (2.35). For ease of notation, we will conveniently redefine the right-hand side of (2.34) simply as $\sum_{\ell_1,m_1,\ell_2,m_2} S_{H_0}^{m_1 m_2} S_{H_0}^{\ell_1 \ell_2}(r) \equiv S_{H_0}^{\ell m}$, in such a way that (2.34) now reads

$$ \partial_r^2 \mathcal{H}_0^{\ell m}(r) + \frac{2r-r_s}{r^2 f(r)} \partial_r \mathcal{H}_0^{\ell m}(r) - \frac{\lambda r^2 f(r) + r_s^2}{r^4 f(r)^2} \mathcal{H}_0^{\ell m}(r) = S_{H_0}^{\ell m}(r). \quad (6.1) $$

The solution is given by (3.21), with $H_{0,1}^{\ell}$ and $H_{0,2}^{\ell}$ provided in eqs. (3.2) and (3.3), i.e.,

$$ \begin{aligned} \mathcal{H}_0^{\ell m}(r) = \frac{r_s^2}{W_{H_0}^{\ell}} & \left[ Q_{\ell}^2(2r/r_s - 1) \int_1^{r/r_s} dy \, y(y-1) P_{\ell}^2(2y-1) S_{H_0}^{\ell m}(y) \right. \\ & \left. - P_{\ell}^2(2r/r_s - 1) \int_{r/r_s}^{r} dy \, y(y-1) Q_{\ell}^2(2y-1) S_{H_0}^{\ell m}(y) \right], \end{aligned} \quad (6.2) $$

where we have redefined the integration variable as $y \equiv r'/r_s$. The goal is to show that (6.2) does not contain any $1/r^n$ falloff, with $n > 1$, for generic $\ell$.

To this end, let us relate $P_{\ell}^2(2y-1)$ and $Q_{\ell}^2(2y-1)$ in the above expression to the Legendre functions of the first kind $P_{\ell}(2y-1)$. Using well known identities,<fnref n="20" /> we can write

$$ P_{\ell}^2(2y-1) = -y(y-1)\partial_y^2 P_{\ell}(2y-1), \quad (6.3) $$

$$ Q_{\ell}^{2}(2y-1) = \frac{1}{2} P_{\ell}^{2}(2y-1) \log \left( \frac{y}{y-1} \right) + \frac{T_{\ell}(y)}{y(y-1)}, \quad (6.4) $$

where $T_{\ell}(y)$ is a superposition of Legendre polynomials of the first kind and derivatives thereof,

$$ \begin{aligned} T_{\ell}(y) = & -\frac{1}{2}(2y-1)P_{\ell}(2y-1) + y(y-1)\partial_y P_{\ell}(2y-1) \\ & - (y-1)^2 y^2 \partial_y^2 \left( \sum_{k=1}^{\ell} \frac{1}{k} P_{k-1}(2y-1) P_{\ell-k}(2y-1) \right). \end{aligned} \quad (6.5) $$

$$ P_{\ell}^{2}(x) = (1-x^{2}) \frac{d^{2}P_{\ell}(x)}{dx^{2}}, \qquad Q_{\ell}^{2}(x) = (1-x^{2}) \frac{d^{2}Q_{\ell}(x)}{dx^{2}}, $$

$$ Q_{\ell}(x) = \frac{1}{2} P_{\ell}(x) \log \left( \frac{1+x}{x-1} \right) - R_{\ell}(x), \qquad R_{\ell}(x) \equiv \sum_{k=1}^{\ell} \frac{1}{k} P_{k-1}(x) P_{\ell-k}(x), $$

$$ Q_{\ell}^{2}(x) = \frac{1}{2} P_{\ell}^{2}(x) \log \left( \frac{1+x}{x-1} \right) + \frac{2x}{1-x^{2}} P_{\ell}(x) + 2P_{\ell}'(x) - (1-x^{2})R_{\ell}''(x). $$

<footnote n="20">²⁰We recall the identities</footnote>
The complete expression of $T_\ell(y)$ is not particularly important in what follows. For the moment, it will be enough to observe that $T_\ell(y)$ is a finite polynomial with only positive powers of $y$. In particular, $T_\ell(y)$ is generally nonzero at $y=0$ and $y=1$. A second important ingredient for our argument is the form of the source $S_{H_0}^{\ell m}$. Plugging the linear-order physical solutions for $H_0$ (eq. (3.6)) and $h_0$ (eq. (3.13)) in the source (2.35), it is easy to see that $S_{H_0}^{\ell m}$ is a rational function of the form

$$S_{H_0}^{\ell m}(r) = \frac{\text{polynomial}(r)}{(r-r_s)^2}. \qquad (6.6)$$

Then, by using eq. (6.4) in (6.2) we can write the second-order solution of eq. (6.1) as

$$^{(2)}H_0^{\ell m}(r) = \frac{r_s^2}{W_{H_0}^{\ell}} \left\{ \begin{aligned} & \frac{1}{2} P_{\ell}^2(2r/r_s - 1) \left[ \log\left(\frac{r}{r-r_s}\right) \int_1^{r/r_s} dy \mathcal{P}^{\ell m}(y) - \int_{r/r_s}^{r/r_s} dy \log\left(\frac{y-1}{y-1}\right) \mathcal{P}^{\ell m}(y) \right] \\ & + \frac{T_{\ell}(r)}{r(r-r_s)} \int_1^{r/r_s} dy \mathcal{P}^{\ell m}(y) - P_{\ell}^2(2r/r_s - 1) \int_{r/r_s}^{r/r_s} dy T_{\ell}(y) S_{H_0}^{\ell m}(y) \end{aligned} \right\}, \quad (6.7)$$

where eqs. (6.3) and (6.6) imply that the quantity

$$\mathcal{P}^{\ell m}(y) \equiv y(y-1)P_{\ell}^{2}(2y-1)S_{H_0}^{\ell m}(y) \qquad (6.8)$$

appearing in the integrands is a finite polynomial. Let us first focus on the square bracket on the first line of eq. (6.7). By writing $\mathcal{P}^{\ell m}(y)$ in the second integral as $\mathcal{P}^{\ell m}(y) = \partial_y \int_1^y dy' \mathcal{P}^{\ell m}(y')$ and integrating by parts, one can remove the term containing $\log(\frac{r}{r-r_s})$ and write the square bracket as $\int_{r/r_s}^{r/r_s} \frac{dy}{y(y-1)} \int_1^y dy' \mathcal{P}^{\ell m}(y') = -\log(r/r_s) \int_0^1 dy \mathcal{P}^{\ell m}(y) + \dots$, where the ellipses denote contributions that are in the form of a finite polynomial in $r/r_s$. Moving to the second line, since $\mathcal{P}^{\ell m}(y)$ is a finite polynomial in $r/r_s$, the integral $\int_1^{r/r_s} dy \mathcal{P}^{\ell m}(y)$ is a polynomial that vanishes at $r=r_s$, cancelling the prefactor $1/(r-r_s)$. As a result, the first term in the second line of (6.7) has the structure $\frac{1}{r} \times \text{polynomial}(\frac{r}{r_s})$. Note that, besides positive powers of $r/r_s$, this potentially contains a $1/r$ term, given by $\frac{T_\ell(0)}{rr_s} \int_0^1 dy \mathcal{P}^{\ell m}(y)$.

Finally, let us consider the last integral in the curly brackets of eq. (6.7). We are going to show that it is a finite polynomial in $r/r_s$. Given the $\frac{1}{(r-r_s)^2}$ factor in the source, see eq. (6.6), it is not immediately obvious that this term does not have any decaying tail at far distances. To show this, we observe that it is only the part of the source involving the odd modes that generates the $\frac{1}{(r-r_s)^2}$ factor, while the part of the source coming from the even modes is always a polynomial. In particular, inserting the odd linear solutions into (2.35), one can recast the part of the source from the odd modes in the form

$$S_{H_0}^{\ell m} \supset \frac{1}{(r-r_s)^2} \sum_{\ell_1,m_1,\ell_2,m_2} \mathcal{E}_{\ell_1 m_1}^{-} \mathcal{E}_{\ell_2 m_2}^{-} \sum_{i=0}^{\ell_1-2} \sum_{j=0}^{\ell_2-2} \sum_{n=0}^{4} C_{ij,\ell m}^{n} \left(\frac{r}{r_s}\right)^{i+j+n}, \quad (6.9)$$

where $i$ and $j$ run over the series representation of the linear odd solution (see eq. (3.12)), while the sum over $n$ stems from the additional $r/r_s$ terms in the source (2.35). The coefficients $C_{ij,\ell m}^n$ can be extracted from the general expression of the source (2.35). Moreover, the terms in $T_\ell(y)$ that might give rise to a tail in $^{(2)}H_0^{\ell m}$ after integration are those in the first line on
the right-hand side of eq. (6.5)—the second line is proportional to $(y-1)^2$, which cancels the singular prefactor in the source. Using the series representation of the Legendre polynomial,

$$P_{\ell}(2y-1) = (-1)^{\ell} \sum_{k=0}^{\ell} \binom{\ell}{k} \binom{\ell+k}{k} (-y)^k, \quad (6.10)$$

and (6.9), the last integral in eq. (6.7) can be recast as$^{21}$

$$\int^{r/r_s} dy T_\ell(y) S_{H_0}^{\ell m}(y) = \frac{1}{2} \log\left(\frac{r}{r_s} - 1\right) \sum_{\ell_1, m_1, \ell_2, m_2} \mathcal{E}_{\ell_1 m_1}^{-} \mathcal{E}_{\ell_2 m_2}^{-} \\ \times \sum_{i=0}^{\ell_1-2} \sum_{j=0}^{\ell_2-2} \sum_{n=0}^{4} C_{ij, \ell m}^{n} [\ell(\ell+1)-(i+j+n)-2]+\dots, \quad (6.11)$$

where, as usual, we have dropped polynomials in $r/r_s$. For fixed values of $\ell$ and $m$, using the explicit expressions for the coefficients $C_{ij,\ell m}^n$ and summing over $n$, we get that the second line of (6.11) is proportional to

$$\left\{ 3\ell_1(\ell_1+1)\ell_2(\ell_2+1)^0 I_{\ell\ell_1\ell_2}^{m\ell_1m\ell_2} + [2+\ell(\ell+1)-2\ell_1(\ell_1+1)-2\ell_2(\ell_2+1)]^1 I_{\ell\ell_1\ell_2}^{m\ell_1m\ell_2} + 2^2 I_{\ell\ell_1\ell_2}^{m\ell_1m\ell_2} \right\} \\ \times \sum_{i=0}^{\ell_1-2\ell_2-2} \sum_{j=0}^{\ell_2-2} (-1)^{i+j} \frac{(\ell_1+i+2)! (\ell_2+j+2)! 2(i+j+6)-\ell_1(\ell_1+1)-\ell_2(\ell_2+1)}{(\ell_1-i-2)! (\ell_2-j-2)! i!(i+3)!j!(j+3)!} \quad (6.12)$$

One can check that the sums on the second line of (6.12) admit an analytic result, which is identically zero for generic, unfixed, values of $\ell_1$ and $\ell_2$. We stress that the vanishing of this term is not a consequence of the selection rules enforced by the first line in (6.12) but is a consequence of the specific radial dependence of the source.

All in all, we conclude that the second-order solution in the even sector has the structure

$$(^{(2)}H_0^{\ell m}(r) = -\frac{1}{W_{H_0}^{\ell}} \left( r_s^2 P_{\ell}^2(2r/r_s - 1) \log(r/r_s) - \frac{r_s}{r} T_{\ell}(0) \right) \int_0^1 dy \mathcal{P}^{\ell m}(y) + \text{polynomial}(r/r_s). \quad (6.13)$$

This shows that the second-order even modes do not contain $1/r^n$ falloffs with $n > 1$ when expanded at large $r$. Although we cannot provide a general proof here, we find that the integral above vanishes exactly,

$$\int_0^1 dy \mathcal{P}^{\ell m}(y) = \int_0^1 dy y(y-1) P_{\ell}^2(2y-1) S_{H_0}^{\ell m}(y) = 0, \quad (6.14)$$

$$y^m - y_0^m = (y - y_0)(y^{m-1} + y^{m-2}y_0 + \dots + yy_0^{m-2} + y_0^{m-1}),$$

with $y_0 = 1$, to first rewrite

$$\frac{y^s}{y-1} = \frac{1}{y-1} + \dots, \quad \frac{(2y-1)y^s}{(y-1)^2} = \frac{1}{(y-1)^2} + \frac{s+2}{y-1} + \dots,$$

up to simple polynomials in $y$. Then, one can perform the integral in $y$: the integral $\int^{r/r_s} \frac{dy}{(y-1)^2}$ generates terms with an inverse power of $r-r_s$, but this gets canceled when multiplied by the associated Legendre polynomial $P_{\ell}^2(2r/r_s - 1)$ (see eq. (6.3)); the integral $\int^{r/r_s} \frac{dy}{y-1}$ produces instead the log in (6.11). Finally, resumming in $k$ yields the combination $\frac{1}{2} C_{ij}^n [\ell(\ell+1)-(i+j+n)-2]$.

<footnote>$^{21}$To obtain (6.11), one can plug (6.10) into (6.5) and (6.9), and use the identity</footnote>
for fixed values of $\ell$, $\ell_1$ and $\ell_2$, after using the explicit expression for $S_{H_0}^{\ell m}$. This implies that neither the $1/r$ nor the $\log(r/r_s)$ terms are present in (6.13), and that the second-order even solution $^{(2)}H_0^{\ell m}$ is indeed just a finite polynomial.

## 6.1.2 Odd sector

Let us now turn to eq. (2.38) for the second-order odd field $^{(2)}h_0^{\ell m}(r)$. The quadratic source on the right-hand side, that here we denote for simplicity as $\sum_{\ell_1,m_1,\ell_2,m_2} S_{h_0}^{m_1 m_2}(r) = S_{h_0}^{\ell m}(r)$, is of mixed type, given by the product of a linear even mode and a linear odd one (see eq. (2.39)). The solution to (2.38) is given by eq. (3.22):

$$ ^{(2)}h_0^{\ell m}(r) = \frac{r_s^2}{W_{h_0}^{\ell}} \left[ h_{0,4}^{\ell}(r) \int_1^{r/r_s} dy h_{0,5}^{\ell}(y) S_{h_0}^{\ell m}(y) - h_{0,5}^{\ell}(r) \int^{r/r_s} dy h_{0,4}^{\ell}(y) S_{h_0}^{\ell m}(y) \right], \quad (6.15) $$

where $h_{0,4}^{\ell}(r)$ and $h_{0,5}^{\ell}(r)$ can be read off from eqs. (3.10) and (3.11). It is worth noticing that the odd source $S_{h_0}^{\ell m}(r)$ computed on the linear field solutions $^{(1)}h_0^{\ell m}(r)$ and $^{(1)}H_0^{\ell m}(r)$ is a finite polynomial in $r$ and vanishes in $r = 0$: $S_{h_0}^{\ell m}(r) = r \cdot \text{polynomial}(r)$.

One important ingredient in the argument for the even case was the possibility of recasting the non-analytic terms in the quadratic solution in the form of the first line in eq. (6.7). This was a direct consequence of eq. (6.4), i.e. of the fact that all the $x$-dependence in the coefficient of the log in $Q_\ell^2$ is fully determined by the first Legendre solution $P_\ell^2$. Note that this is not special of the Legendre polynomials, but it holds more in general for solutions of generic Fuchsian-type second-order ordinary differential equations [91].$^{22}$ In the case of (2.38), one can indeed show that the independent solution $h_{0,4}^{\ell}(r)$ in eq. (3.10) has the structure [91]

$$ h_{0,4}^{\ell}(r) = \alpha_{\ell} \log \left( \frac{r}{r-r_s} \right) h_{0,5}^{\ell}(r) + \frac{\zeta_{\ell}(r)}{r}, \quad (6.16) $$

where $h_{0,5}^{\ell}(r)$ is the polynomial in eq. (3.12), while $\alpha_\ell$ are constant, $\ell$-dependent coefficients, and $\zeta_\ell(r)$ is an analytic function of $r$ in the neighbourhoods of $r = 0$ and $r = r_s$. In particular, it follows from the properties of the hypergeometric function that $\zeta_\ell(r)$ is a finite polynomial in $r$ [79].

Thus, we can follow the same logic as in the case above and write

$$ ^{(2)}h_0^{\ell m}(r) = \frac{\alpha_{\ell} r_s^2}{W_{h_0}^{\ell}} h_{0,5}^{\ell}(r) \left[ \log\left(\frac{r}{r-r_s}\right) \int_1^{r/r_s} dy h_{0,5}^{\ell}(y) S_{h_0}^{\ell m}(y) - \int^{r/r_s} dy \log\left(\frac{y}{y-1}\right) h_{0,5}^{\ell}(y) S_{h_0}^{\ell m}(y) \right] + \dots \quad (6.17) $$

where we kept only terms that could potentially give rise to inverse powers of $r$ at large $r$. From the same algebraic manipulations that we used for the first line in (6.7), we can then conclude that $^{(2)}h_0^{\ell m}(r)$ has no induced decaying $1/r^n$ falloff (with $n > 1$) at $r = \infty$. As for the second-order $^{(2)}H_0^{\ell m}$, $^{(2)}h_0^{\ell m}$ takes the form of a polynomial in $r/r_s$ with only positive integer powers, except for a $1/r$ and a $\log(r)$ term, which will not affect our argument in section 6.2. Nevertheless, the coefficients of these terms are proportional to the integral

$$ \int_0^1 dy h_{0,5}^{\ell}(y) S_{h_0}^{\ell m}(y) = 0, \quad (6.18) $$

<footnote>$^{22}$This happens whenever the roots of the indicial equation for a given regular singular point ($r = 0$ and $r = r_s$ in the case of (2.38)) differ by an integer number [91].</footnote>
which we find to be vanishing after fixing specific values of $\ell_1$, $\ell_2$ and $\ell$, although we cannot provide a general proof here.

## 6.2 Matching to all orders in derivatives

In this section, we will perform the matching between the worldline EFT and the full solutions for $^{(2)}$$H_0^{\ell m}(r)$ and $^{(2)}$$h_0^{\ell m}(r)$ discussed above, for all values of $\ell$. As opposed to section 5, where we studied the matching for $\ell = 2$ tidal fields and showed how to reconstruct the nonlinear tidal solution at the next-to-leading order in the large-$r$ expansion, here we will keep $\ell$ generic, while working at the leading order in $r_s$, i.e., neglecting the contribution from the point-particle action $S_{pp}$ in eq. (4.2). This will be enough to conclude that all quadratic Love number couplings vanish in the EFT in the case of Schwarzschild black holes.

After imposing parity invariance, the generalization of the cubic Love number action (4.10) to higher orders in derivatives is

$$ S_T = \int d\tau \left\{ \sum_{\ell=2}^{\infty} \left[ \lambda_{\ell}^{(E)} E_{i_1 \dots i_{\ell}} E^{i_1 \dots i_{\ell}} + \lambda_{\ell}^{(B)} B_{i_1 \dots i_{\ell}} B^{i_1 \dots i_{\ell}} \right] \right. \\ \left. \qquad + \sum_{n,m,l} \lambda_{nml}^{(E^3)} E_{i_1 \dots i_l j_1 \dots j_n} E^{i_1 \dots i_l} B_{k_1 \dots k_m} E^{j_1 \dots j_n} B_{k_1 \dots k_m} \right. \\ \left. \qquad + \sum_{n,m,l} \lambda_{nml}^{(EB^2)} E_{i_1 \dots i_l j_1 \dots j_n} B^{i_1 \dots i_l} B_{k_1 \dots k_m} B^{j_1 \dots j_n} B_{k_1 \dots k_m} + \mathcal{O}(E^4, E^2 B^2, B^4) \right\}, \quad (6.19) $$

where $E_{i_1 \dots i_n}$ and $B_{i_1 \dots i_n}$ are respectively defined in eqs. (4.4) and (4.5), while $n$, $m$ and $l$ are non-negative integers. The cubic terms in the second and third lines of (6.19) correspond to couplings of $\ell_1 = l+n$, $\ell_2 = l+m$ and $\ell_3 = n+m$ multipoles. Since $E_{i_1 \dots i_n}$ and $B_{i_1 \dots i_n}$ are STF tensors, in order to avoid overcounting, the sum in the second line of (6.19) is constrained to run over $l \ge n \ge m \ge 0$, which translates into $\ell_1 \ge \ell_2 \ge \ell_3$ (with $\ell_1, \ell_2, \ell_3 \ge 2$). The sum in the third line of (6.19) runs instead over non-negative integers $n, m$ and $l$, such that $l \ge n$, with again $\ell_1, \ell_2, \ell_3 \ge 2$.

To perform the matching, we will first solve the quadratic equations obtained from the EFT action at leading order in $r_s$ and compute the induced response. We already know that the linear Love numbers vanish [14–16, 20–22, 24], so we can set $\lambda_{\ell}^{(E)} = 0 = \lambda_{\ell}^{(B)}$ in (6.19). This will considerably simplify the analysis below. In particular, as we will see, it will be enough to work with the linearized components of the Weyl tensor. In the static case, they read

$$ E_{ij} = -\frac{1}{2} \partial_i \partial_j \delta g_{tt} + \mathcal{O}(\delta g^2), \quad B_{ij} = -\frac{1}{2} \varepsilon_{ti}^{kl} \partial_j \partial_k \delta g_{it} + \mathcal{O}(\delta g^2), \quad (6.20) $$

where here $\delta g_{\mu\nu} \equiv g_{\mu\nu} - \eta_{\mu\nu}$. Then, the variation of the action $S = S_{EH} + S_T$ with respect to $g_{\mu\nu}$, including only quadratic corrections, yields schematically

$$ \begin{aligned} M_{Pl}^2 G_{\mu\nu} & \sim \left[ \sum_{\ell_1, \ell_2, \ell_3} \lambda_{\ell_1 \ell_2 \ell_3}^{(E^3)} (-1)^{\ell_1} \partial_{i_1} \cdots \partial_{i_{\ell_1}} (\delta^{(3)}(\mathbf{x}) E_{L_2} E_{L_3}) \right. \\ & \qquad \left. + \sum_{\ell_1, \ell_2, \ell_3} \lambda_{\ell_1 \ell_2 \ell_3}^{(EB^2)} (-1)^{\ell_1} \partial_{i_1} \cdots \partial_{i_{\ell_1}} (\delta^{(3)}(\mathbf{x}) B_{L_2} B_{L_3}) + \text{perms} \right] \delta_{\mu}^{t} \delta_{\nu}^{t} \\ & \qquad + \left[ \sum_{\ell_1, \ell_2, \ell_3} \lambda_{\ell_1 \ell_2 \ell_3}^{(EB^2)} (-1)^{\ell_2} \varepsilon_{tj}^{i_1 i_2} \partial_{i_1} \cdots \partial_{i_{\ell_2}} (\delta^{(3)}(\mathbf{x}) E_{L_1} B_{L_3}) + \text{perms} \right] \delta_{(\mu}^{j} \delta_{\nu)}^{t} \\ & \equiv J^{+}(\mathbf{x}) \delta_{\mu}^{t} \delta_{\nu}^{t} + J_{j}^{-}(\mathbf{x}) \delta_{(\mu}^{j} \delta_{\nu)}^{t}, \end{aligned} \quad (6.21) $$
where the precise form of the contractions among the tensors on the right-hand side will not be important for our argument. We are interested in solving the previous perturbation equations on flat space at second order in the fields. Since the right-hand side is already quadratic in the Weyl tensor, as anticipated, the linear expansions (6.20) will suffice.

Equations (6.21) can be also written as

$$M_{Pl}^2 R_{\mu\nu} \sim J^+(\boldsymbol{x}) \left( \delta_{\mu}^t \delta_{\nu}^t + \frac{1}{2} \eta_{\mu\nu} \right) + J_j^-(\boldsymbol{x}) \delta_{(\mu}^j \delta_{\nu)}^t. \qquad (6.22)$$

Let us decompose the metric fluctuation as

$$\delta g_{\mu\nu} = H_{\mu\nu} + h_{\mu\nu}, \qquad (6.23)$$

where $H_{\mu\nu}$ represents the external tidal field, while $h_{\mu\nu}$ is the response. Each field is in turn expanded in perturbation theory as $H_{\mu\nu} = (1)H_{\mu\nu} + (2)H_{\mu\nu} + \dots$, and $h_{\mu\nu} = (1)h_{\mu\nu} + (2)h_{\mu\nu} + \dots$, to be plugged into (6.22). By definition, $H_{\mu\nu}$ solves the homogeneous vacuum Einstein equations i.e., $R_{\mu\nu}[h=0] = 0$. Since the linear response vanishes, $(1)h_{\mu\nu} = 0$, the left-hand side of (6.22) is simply $M_{Pl}^2(1)\delta R_{\mu\nu}[(2)h]$, where $(1)\delta R_{\mu\nu}$ denotes the first-order variation of the Ricci tensor, evaluated on the second-order perturbation $(2)h_{\mu\nu}$. Instead, the right-hand side of (6.22) depends quadratically on $(1)H_{\mu\nu}$, after the Weyl tensors are evaluated on the linear tidal field solution. All in all, (6.22) becomes at second order an inhomogeneous equation for $(2)h_{\mu\nu}$, with a fully known source, which can be solved through standard Green's function methods.

To this end, it is convenient to choose the de Donder gauge, defined by

$$\partial^{\mu} \delta g_{\mu\nu} - \frac{1}{2} \partial_{\nu} \delta g = 0, \qquad (6.24)$$

with $\delta g \equiv \eta^{\mu\nu}\delta g_{\mu\nu}$.<fnref n="23" /> In such a gauge, a certain number of simplifications occur on flat space, in particular $\delta^{(1)}R_{\mu\nu} = -\frac{1}{2}\Box\delta g_{\mu\nu}$. For the static linear tidal field — to plug in the expressions (6.20), which we will need to evaluate the sources $J^{\pm}$ in (6.22)—we shall thus take [22, 24]

$$(1) H_{tt}(\vec{x}) = \mathcal{E}_{i_1 \dots i_\ell}^{+} x^{i_1} \dots x^{i_\ell}, \qquad (1) H_{tj}(\vec{x}) = \mathcal{E}_{j|i_1 \dots i_\ell}^{-} x^{i_1} \dots x^{i_\ell}. \qquad (6.25)$$

The relevant equations in (6.22) boil down to

$$\vec{\nabla}^2 (2) h_{tt} = -M_{Pl}^{-2} J^+(x), \qquad \vec{\nabla}^2 (2) h_{tj} = -M_{Pl}^{-2} J_j^-(x), \qquad (6.26)$$

where the sources $J_j^{\pm}(x)$, evaluated on the tidal fields (6.25), are proportional to derivatives of the delta-function [22]. The solutions can be straightforwardly obtained by going to momentum following [22, 24], and are found to scale at large $r \equiv |x|$ as $(2)h_{tt} \sim 1/r^{\ell+1}$, and similarly for $(2)h_{tj}$, with $\ell$ corresponding to the multipole number of the induced second-order field and such that $\ell \ge 2$. (This scaling can be also understood by simple power counting, by looking at the $x$-dependence in the sources in (6.21).) These solutions should then be

<footnote n="23">²³Here the gauge fixing for $H_{\mu\nu}$ is slightly different from the one that we chose in section 4, where $H_{\mu\nu}$ was in RW gauge.</footnote>
compared with the expressions for $^{(2)}$H$_{0}^{lm}$ and $^{(2)}$h$_{0}^{lm}$ in eqs. (6.2) and (6.15), obtained in full general relativity. Note that in general such a comparison requires some care, as the two calculations are done in two different gauges, as we extensively discussed in section 5.3. However, as long as one is interested in the matching at leading order in the flat-space limit, one can easily check that, under the assumptions of vanishing linear response (i.e., $^{(1)}$h$_{µν}$ = 0) and static limit, the metric perturbations $^{(2)}$δg$_{tt}$ and $^{(2)}$δg$_{tj}$ in the RW and de Donder (dD) gauge coincide, i.e.,$^{24}$

$$ ^{(2)}\delta g_{tt}^{\text{RW}} = ^{(2)}\delta g_{tt}^{\text{dD}}, \quad ^{(2)}\delta g_{tj}^{\text{RW}} = ^{(2)}\delta g_{tj}^{\text{dD}}. \qquad (6.27) $$

This makes the matching of $^{(2)}$h$_{tt}$ and $^{(2)}$h$_{tj}$ from (6.26) with their counterpart solutions in the full theory straightforward. In particular, since we showed in section 6.1 that $^{(2)}$H$_{0}^{lm}$ and $^{(2)}$h$_{0}^{lm}$ take a simple form with no 1/r$^{n}$ (n > 1) falloff at large r, one then concludes that all the Love number coefficients $λ^{(E3)}$ and $λ^{(EB2)}$ in (6.19) vanish to all orders in derivatives. Subleading corrections in r$_{s}$ to the nonlinear tidal field solution are expected to be reconstructed by including contributions from the point-particle action S$_{pp}$, as we did in section 5 for an external quadrupolar tidal field.

# 7 Conclusions

In this work, we provided the first comprehensive analysis of static quadratic perturbations and quadratic Love numbers of Schwarzschild black holes in general relativity, extending the previous results of [52]. We studied both even and odd sectors, with no restriction on the angular dependence of the modes. We derived the second-order equations in Regge-Wheeler gauge and obtained an analytic expression for the sources (see section 2.3). We showed that, when evaluated on the linear solutions, the sources take a simple form, which is in turn inherited by the induced second-order metric perturbations. In full generality, we proved that the latter take the form of simple polynomials with only positive powers of r, except for potentially a 1/r and a log(r) term; although we could not provide a general proof of the absence thereof, explicit checks for fixed values of the angular momentum quantum numbers of the modes show that their coefficients generically vanish (see end of sections 6.1.1 and 6.1.2). We then used the worldline effective field theory [10, 38, 39, 55–59] to defined the nonlinear tidal deformability and quadratic Love numbers of compact objects. By performing an explicit matching with the full solutions derived in section 3, we concluded that all quadratic (even and odd) Love numbers vanish for Schwarzschild black holes.

Our results raise a series of interesting questions on the fundamental nature of black holes in general relativity. First of all, the vanishing of the quadratic Love numbers appears as a fine tuning in the worldline EFT, mirroring the case of the linear Love numbers [39, 40]. In the latter case, symmetries acting on the static sector of the linearized perturbations have been found [23, 24, 41] and proposed as a resolution to the puzzle. Our findings here are suggestive that the ladder symmetries of [23, 24, 41] admit an extension to nonlinear order,

<footnote>$^{24}$This can be seen also by using eqs. (5.34)–(5.35) and carefully taking the r$_{s}$ → 0 limit in eq. (5.36). Indeed, by dimensional analysis the scale of the tidal field is set by $E_{lm}/r_{s}^{l}$, which should be kept finite when taking this limit. In this case $^{(1)}$ξ$^{µ}$ = 0, barring residual large gauge transformations satisfying $^{(1)}$ξ$^{µ}$ = 0, which are allowed in de Donder gauge.</footnote>
which could be used to explain the vanishing of the nonlinear tidal response of black holes in general relativity [92, 93]. In addition, the simplicity of the obtained results, in particular the form of the sources and of the second-order solutions (see section 6.1) suggests that one might be able to reorganize perturbation theory in a way that takes full advantage of the symmetries of the system and simplifies the intermediate steps. One notable example of this is given by axisymmetric (even) perturbations: for this subset of modes, the full metric, including perturbations, can be recast in the Weyl form, which has been used to show that there is no static, axisymmetric induced response to all order in perturbations [21, 51, 92, 93]. This applies to the even sector and for modes with zero magnetic quantum number. On the other hand, we showed here that also the odd sector displays a very simple structure. It will be interesting to understand this further, and find a more systematic way to perform the matching, and generalize our results to cubic and higher orders in perturbations beyond axisymmetry. Finally, it will be interesting to study subleading tidal effects, in particular nonlinearities, for different black hole solutions, such as Kerr spacetime, black holes with charge, and solutions in higher spacetime dimensions [22, 27, 28, 94, 95]. We leave all these aspects for future investigations.

## Acknowledgments

It is a pleasure to thank Eric Poisson for interesting comments on our paper. The research of L.S. has been funded, in part, by the French National Research Agency (ANR) under project ANR-24-CE31-1097-01, and by the Programme National GRAM of CNRS/INSU with INP and IN2P3 co-funded by CNES.

M.M.R. is partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy — EXC 2121 "Quantum Universe" — 390833306, by the ERC Consolidator Grant "Precision Gravity: From the LHC to LISA" provided by the European Research Council (ERC) under the European Union's H2020 research and innovation programme (grant agreement No. 817791) and by the DOE DE-SC0011941.

## A Integration of tensor harmonics

In this appendix we provide details on how to perform the angular integrations of tensor harmonics, eqs. (2.26)–(2.30). We remind the reader that $\nabla_A$ is the covariant derivative on the unit 2-sphere $S^2$, $\epsilon_{AB}$ is the Levi-Civita tensor (with non-vanishing components $\epsilon_{\theta\phi} = -\epsilon_{\phi\theta} = \sin\theta$) and that we use the metric on $S^2$, $\gamma_{AB}$, to raise and lower indices.

We will use the so-called “pure-spin” harmonics approach, following ref. [63]. We construct symmetric trace-free (STF) tensor harmonics by acting on scalar spherical harmonics with covariant derivatives $\nabla_A$ and with the Levi-Civita symbol $\epsilon_{AB}$. In particular, we define even and odd STF s-index tensors of spin $s \ge 1$ ($s$ is integer), angular momentum $\ell$ and magnetic number $m$, respectively as

$$Z_{\ell A_1 \dots A_s}^m \equiv (\nabla_{A_1} \dots \nabla_{A_s} Y_{\ell}^m)_{\text{STF}}, \qquad (A.1)$$

$$X_{\ell A_1 \dots A_s}^m \equiv (\epsilon_{A_1}^B Z_{\ell BA_2 \dots A_s}^m)_{\text{STF}}, \qquad (A.2)$$
where by the index STF we denote the part. We discuss the case $s = 0$ separately below.
We can rewrite the integrals in compact form as

$$ {}^s I_{\ell\ell_1\ell_2}^{m_1m_2} \equiv \int d\Omega Y_{\ell}^m Z_{\ell_1 A_1 \dots A_s}^{m_1} Z_{\ell_2 A_1 \dots A_s}^{m_2}, \quad (A.3) $$

$$ {}^s J_{\ell\ell_1\ell_2}^{m_1m_2} \equiv \int d\Omega Y_{\ell}^m X_{\ell_1 A_1 \dots A_s}^{m_1} Z_{\ell_2 A_1 \dots A_s}^{m_2}, \quad (A.4) $$

Then, we can express the products of two spin-$s$ spherical tensors with all the indices contracted, into spherical harmonics, i.e. [63],

$$ Z_{\ell_1 A_1 \dots A_s}^{m_1} Z_{\ell_2}^{m_2, A_1 \dots A_s} = (-1)^s \sum_{\ell = |\ell_1 - \ell_2|}^{\ell_1 + \ell_2} [1 + (-1)^{\ell + \ell_1 + \ell_2}] E_{-s \ell_1 m_1 \ell}^{\ell_2 m_2} Y_{\ell}^{m_1+m_2}, \quad (A.5) $$

$$ X_{\ell_1 A_1 \dots A_s}^{m_1} Z_{\ell_2}^{m_2, A_1 \dots A_s} = -i(-1)^s \sum_{\ell = |\ell_1 - \ell_2|}^{\ell_1 + \ell_2} [1 - (-1)^{\ell + \ell_1 + \ell_2}] E_{-s \ell_1 m_1 \ell}^{\ell_2 m_2} Y_{\ell}^{m_1+m_2}, \quad (A.6) $$

where the factor $[1 \pm (-1)^{\ell+\ell_1+\ell_2}]$ ensures that only terms with a specific $(\ell + \ell_1 + \ell_2)$-parity (even for + and odd for -) are picked up in the sum. The coefficients of these expansions, $E_{-s \ell_1 m_1 \ell}^{\ell_2 m_2}$, can be written in terms of Clebsch-Gordan coefficients $\langle \ell_1 m_1 \ell_2 m_2 | \ell m \rangle$. Specifically, they are given by

$$ E_{s_1 \ell_1 m_1 \ell}^{s_2 \ell_2 m_2} = \frac{k(\ell_1, |s_1|)k(\ell_2, |s_2|)}{k(\ell, |s_1 + s_2|)} C_{\ell_1 \ell_2 \ell}^{m_1 m_2 (m_1+m_2)} C_{\ell_1 \ell_2 \ell}^{s_1 s_2 (s_1+s_2)}, \quad (A.7) $$

where $C_{\ell_1 \ell_2 \ell}^{m_1 m_2 m} \equiv \langle \ell_1 m_1 \ell_2 m_2 | \ell m \rangle$ and

$$ k(\ell, s) \equiv \sqrt{\frac{(1 + 2\ell)(\ell + s)!}{2^{s+2}\pi(\ell - s)!}}. \quad (A.8) $$

We can now solve the integrals (A.3) and (A.4). Using the expansions eqs. (A.5)-(A.7), the orthogonality condition, $\int d\Omega Y_{\ell}^{m} Y_{\ell'}^{m*} = \delta_{\ell\ell'} \delta^{mm'}$, and the relation

$$ \frac{k(\ell_1, |s|)k(\ell_2, |s|)}{k(\ell, 0)} = 2^{-|s|} \sqrt{\frac{(1 + 2\ell_1)(1 + 2\ell_2)(\ell_1 + |s|)!(\ell_2 + |s|)!}{4\pi(1 + 2\ell)(\ell_1 - |s|)!(\ell_2 - |s|)!}}, \quad (A.9) $$

we obtain,

$$ {}^s I_{\ell\ell_1\ell_2}^{m_1m_2} = \frac{1}{2} [1 + (-1)^{\ell+\ell_1+\ell_2}] {}^s I_{\ell\ell_1\ell_2}^{mm_1m_2}, \quad (A.10) $$

$$ {}^s J_{\ell\ell_1\ell_2}^{m_1m_2} = -i \frac{1}{2} [1 - (-1)^{\ell+\ell_1+\ell_2}] {}^s I_{\ell\ell_1\ell_2}^{mm_1m_2}, \quad (A.11) $$

where

$$ {}^s I_{\ell\ell_1\ell_2}^{m_1m_2} \equiv 2^{1-|s|} (-1)^s \sqrt{\frac{(1+2\ell_1)(1+2\ell_2)(\ell_1+|s|)!(\ell_2+|s|)!}{4\pi(1+2\ell)(\ell_1-|s|)!(\ell_2-|s|)!}} C_{\ell_1\ell_2\ell}^{m_1m_2m} C_{\ell_1\ell_2\ell}^{-ss0}. \quad (A.12) $$

The case $s = 0$ is well known, i.e.,

$$ {}^0 I_{\ell\ell_1\ell_2}^{m_1m_2} \equiv \sqrt{\frac{(1 + 2\ell_1)(1 + 2\ell_2)}{4\pi(1 + 2\ell)}} C_{\ell_1\ell_2\ell}^{m_1m_2m} C_{\ell_1\ell_2\ell}^{000}, \quad (A.13) $$
so the following definition,

$$ {}^{s}\mathcal{I}_{\ell_1\ell_1\ell_2}^{m_1m_2} \equiv (2-\delta_{0s})2^{-s}(-1)^s \sqrt{\frac{(1+2\ell_1)(1+2\ell_2)(l_1+|s|)!(l_2+|s|)!}{4\pi(1+2\ell)(l_1-|s|)!(l_2-|s|)!}} C_{\ell_1\ell_2\ell}^{m_1m_2m} C_{\ell_1\ell_2\ell}^{-ss0}, \quad (A.14) $$

holds for $s \ge 0$.

We note the following symmetry properties of the Clebsch-Gordan coefficients:

$$ C_{\ell_1 \ell_2 \ell}^{m_1 m_2 m} = (-1)^{\ell+\ell_1+\ell_2} C_{\ell_2 \ell_1 \ell}^{m_2 m_1 m} = (-1)^{\ell+\ell_1+\ell_2} C_{\ell_1 \ell_2 \ell}^{-m_1-m_2-m}. \quad (A.15) $$

Since $C_{\ell_1\ell_2\ell}^{-ss0} = C_{\ell_2\ell_1\ell}^{-ss0}$, ${}^s\mathcal{I}_{\ell_1\ell_1\ell_2}^{m_1m_2m}$ has the same symmetries as $C_{\ell_1\ell_2\ell}^{m_1m_2m}$, i.e., ${}^s\mathcal{I}_{\ell_1\ell_1\ell_2}^{m_1m_2m}$ is symmetric under the exchange $\ell_1m_1 \leftrightarrow \ell_2m_2$ if $\ell + \ell_1 + \ell_2$ is even, and antisymmetric if $\ell + \ell_1 + \ell_2$ is odd. Therefore, eqs. (A.10) and (A.11) can be rewritten as

$$ {}^s I_{\ell_1 \ell_2}^{m_1 m_2} = \frac{1}{2} ( {}^s \mathcal{I}_{\ell_1 \ell_2}^{m_1 m_2} + {}^s \mathcal{I}_{\ell_1 \ell_2}^{m_2 m_1} ), \quad (A.16) $$

$$ {}^s J_{\ell_1 \ell_2}^{m_1 m_2} = -\frac{i}{2} ( {}^s \mathcal{I}_{\ell_1 \ell_2}^{m_1 m_2} - {}^s \mathcal{I}_{\ell_1 \ell_2}^{m_2 m_1} ), \quad (A.17) $$

respectively. Note also that, if $\ell + \ell_1 + \ell_2$ is odd, eq. (A.15) implies ${}^s\mathcal{I}_{\ell_1\ell_2\ell}^{000} = 0$ and ${}^s\mathcal{I}_{\ell_1\ell_1\ell}^{m_1m_1m} = 0$. The same properties hold for ${}^s J_{\ell_1\ell_1\ell}^{m_1m_2m}$.

Finally, note that the $s=1$ integrals can be expressed in terms of the $s=0$ ones using integrations by parts. For instance,

$$ \begin{aligned} {}^1 I_{\ell_1 \ell_1 \ell_2}^{m_1 m_2 m} &= \int d\Omega Y_\ell^m \nabla_A Y_{\ell_1}^{m_1} \nabla^A Y_{\ell_2}^{m_2} = \frac{1}{2} \int_{S^2} Y_\ell^m \left[ \Delta_{S^2} (Y_{\ell_1}^{m_1} Y_{\ell_2}^{m_2}) - Y_{\ell_1}^{m_1} \Delta_{S^2} Y_{\ell_2}^{m_2} - Y_{\ell_2}^{m_2} \Delta_{S^2} Y_{\ell_1}^{m_1} \right] \\ &= \frac{1}{2} \left[ \ell_1(\ell_1+1) + \ell_2(\ell_2+1) - \ell(\ell+1) \right]^0 I_{\ell_1 \ell_1 \ell_2}^{m_1 m_2 m}. \end{aligned} \quad (A.18) $$

## B Coefficients

Here we provide the explicit expressions of coefficients appearing in the text, starting from those in the sources in section 2.3. The coefficients appearing in eq. (2.35) are

$$ a_2 = \frac{1}{4}(\lambda-2)(\lambda-\lambda_1-\lambda_2)\lambda, $$

$$ a_1 = -\frac{\lambda^3}{4} + \frac{\lambda^2(4+\lambda_1(4-4\lambda_1+\lambda_1\lambda_2)-\lambda_1\lambda_2)}{2(\lambda_1-2)(\lambda_2-2)} - \frac{\lambda_1\lambda(20-10\lambda_1-6\lambda_2+3\lambda_1\lambda_2)}{2(\lambda_1-2)(\lambda_2-2)} + \frac{1}{2}\lambda_1(\lambda_1-\lambda_2) $$

$$ a_0 = -\frac{\lambda_1(\lambda_1-\lambda_2)}{4} + \frac{\lambda^2(2\lambda_1^2-\lambda_1\lambda_2-4)}{4(\lambda_1-2)(\lambda_2-2)} + \frac{\lambda(8+\lambda_1(4-6\lambda_1+\lambda_1\lambda_2))}{4(\lambda_1-2)(\lambda_2-2)} $$

$$ b_0 = -\frac{(\lambda_1 - \lambda_2)(2 + \lambda_1 - \lambda_2)}{4} + \frac{\lambda(\lambda_1^2 + \lambda_1 + 2\lambda_2 - 6)}{2(\lambda_1 - 2)} - \frac{\lambda^2(3\lambda_1 + 2\lambda_2 - 6)}{4(\lambda_1 - 2)}, $$

$$ c_4 = \frac{1}{2}\lambda_1\lambda_2(\lambda-2)\lambda, $$

$$ c_3 = \frac{1}{4} \{ \lambda^3 - 4\lambda^2 - 12\lambda + 2\lambda [(\lambda_1 + \lambda_2)(5 - 3\lambda) + 2\lambda_1\lambda_2(6 - \lambda)] + (\lambda_1^2 + \lambda_2^2)(-6 + 11\lambda) \\ - 6(\lambda_1^3 + \lambda_2^3) + 12(\lambda_1 + \lambda_2)^2 \}, $$
$$c_2 = \frac{1}{8} \left\{ \begin{aligned} & -5\lambda^3 + 44\lambda^2 + 92\lambda + (\lambda_1 + \lambda_2)(-92 - 120\lambda + 21\lambda^2) - \lambda_1\lambda_2(132 + 56\lambda - 4\lambda^2) \\ & + (\lambda_1^2 + \lambda_2^2)(76 - 35\lambda) + 19(\lambda_1^3 + \lambda_2^3) - \lambda_1\lambda_2(\lambda_1 + \lambda_2) \end{aligned} \right\},$$

$$c_1 = \frac{1}{8} \left\{ \begin{aligned} & 3\lambda^3 - 58\lambda^2 - 112\lambda + (\lambda_1 + \lambda_2)(112 + 164\lambda - 9\lambda^2) + 13(\lambda_1^2 + \lambda_2^2)\lambda + \lambda_1\lambda_2(124 + 16\lambda) \\ & - 7(\lambda_1^3 + \lambda_2^3) - 106(\lambda_1^2 + \lambda_2^2) + \lambda_1\lambda_2(\lambda_1 + \lambda_2) \end{aligned} \right\},$$

$$c_0 = \frac{1}{8} \left\{ 19\lambda^2 + 38\lambda - (\lambda_1 + \lambda_2)(38 + 58\lambda) + 39(\lambda_1^2 + \lambda_2^2) - 40\lambda_1\lambda_2 \right\},$$

$$d_2 = \frac{1}{4}(\lambda^3 + 4\lambda^2 + 20\lambda + 2\lambda_1\lambda^2 + 3\lambda_2\lambda^2 - 5\lambda_1^2\lambda - 9\lambda_2^2\lambda - 22\lambda_1\lambda - 10\lambda_1\lambda_2\lambda - 16\lambda_2\lambda + 2\lambda_1^3 + 5\lambda_2^3 \\ + 18\lambda_1^2 - 2\lambda_1\lambda_2^2 + 12\lambda_2^2 - 20\lambda_1 + 3\lambda_1^2\lambda_2 - 26\lambda_1\lambda_2 - 20\lambda_2),$$

$$d_1 = -\frac{\lambda^3}{4} + \frac{1}{4}(-2\lambda_1 - 3\lambda_2 - 12)\lambda^2 + \frac{\lambda}{4}(5\lambda_1^2 + 10\lambda_2\lambda_1 + 48\lambda_1 + 9\lambda_2^2 + 42\lambda_2 - 36) + \\ + \frac{1}{4}(-2\lambda_1^3 - 3\lambda_2\lambda_1^2 - 36\lambda_1^2 + 2\lambda_2^2\lambda_1 + 46\lambda_2\lambda_1 + 36\lambda_1 - 5\lambda_2^3 - 30\lambda_2^2 + 36\lambda_2),$$

$$d_0 = \frac{11}{4}\lambda^2 + \frac{\lambda}{4}\{22 - 32(\lambda_1 + \lambda_2)\} + \frac{1}{4}\{21(\lambda_1^2 + \lambda_2^2) - 22(\lambda_1 + \lambda_2) - 20\lambda_1\lambda_2\},$$

$$e_2 = -\frac{\lambda^3}{4} + \frac{\lambda^2}{8}(2\lambda_1 + 2\lambda_2 + 4) + \frac{\lambda}{8}(-4\lambda_1 - 4\lambda_2),$$

$$e_1 = \frac{\lambda^3}{4} + \frac{\lambda^2}{8}(-2\lambda_1 - 2\lambda_2 - 12) + \frac{\lambda}{8}(26\lambda_1 + 26\lambda_2 - 16) + \frac{1}{8}(-14\lambda_1^2 + 12\lambda_2\lambda_1 + 16\lambda_1 - 14\lambda_2^2 + 16\lambda_2)$$

$$e_0 = \frac{1}{2}\lambda^2 + \frac{\lambda}{8}(-13\lambda_1 - 13\lambda_2 + 8) + \frac{1}{8}(9\lambda_1^2 - 10\lambda_2\lambda_1 - 8\lambda_1 + 9\lambda_2^2 - 8\lambda_2)$$

$$f_3 = 6 - \lambda + 6\lambda_1,$$

$$f_2 = \frac{1}{2}(3\lambda - 9\lambda_1 - 46),$$

$$f_1 = \frac{1}{2}(56 - \lambda + 3\lambda_1),$$

$$g_2 = -10 + \lambda - \lambda_1 - 3\lambda_2,$$

$$g_1 = 18 - \lambda + \lambda_1 + 3\lambda_2,$$

The coefficients appearing in eq. (2.39) are

$$\begin{aligned} \tilde{a}_2 &= \lambda - \lambda_1 - \lambda_2 + \lambda_1\lambda_2, & \tilde{a}_1 &= \frac{\lambda(2 - \lambda_1 + 2\lambda_2)}{\lambda_1 - 2} + \lambda_1 + \lambda_2 - \lambda_1\lambda_2, & \tilde{a}_0 &= -\frac{\lambda(2 - \lambda_1 + 2\lambda_2)}{2(\lambda_1 - 2)} - \frac{1}{2}(\lambda_1 - \lambda_2), \end{aligned} \qquad (B.1)$$

$$\tilde{b}_1 = -\lambda + \lambda_1 - \lambda_2, \qquad \tilde{b}_0 = -\frac{1}{2} \left( \frac{\lambda(6 - 3\lambda_1 - 2\lambda_2)}{\lambda_1 - 2} + 3(\lambda_1 - \lambda_2) \right),$$

The numerical coefficients appearing in eq. (5.12) are given by

$$c_{\eta}^{(0,+)} = 1, \quad c_{v}^{(0,+)} = \frac{1}{2}, \quad c_{r}^{(0,+)} = \frac{5}{2}, \quad c_{\eta}^{(0,-)} = -15, \quad c_{v}^{(0,-)} = -\frac{21}{2}, \quad c_{r}^{(0,-)} = -\frac{117}{2}, \qquad (B.2a)$$

$$c_{\eta}^{(2,+)} = -\frac{1}{8}, \quad c_{v}^{(2,+)} = -\frac{5}{8}, \quad c_{r}^{(2,+)} = -\frac{7}{8}, \quad c_{\eta}^{(2,-)} = \frac{63}{8}, \quad c_{v}^{(2,-)} = \frac{99}{8}, \quad c_{r}^{(2,-)} = \frac{105}{8}, \qquad (B.2b)$$

$$c_{\eta}^{(4,+)} = \frac{97}{36}, \quad c_{v}^{(4,+)} = \frac{85}{18}, \qquad c_{\eta}^{(4,-)} = -\frac{7}{36}, \quad c_{v}^{(4,-)} = -\frac{73}{18}, \qquad (B.2c)$$

$$d^{(2)} = -\frac{1}{2}, \quad d^{(4)} = \frac{8}{5}. \qquad (B.2d)$$
C From spherical to Cartesian coordinates

The goal of this appendix is to express the tidal field in the RW gauge in a Lorentz-covariant form, making it easy to express it in Cartesian coordinates, which is the most convenient form for Feynman diagrams. As discussed in section 5.1, the tidal field represents a perturbation around Minkowski spacetime, i.e., $H_{\mu\nu} = \bar{g}_{\mu\nu} - \eta_{\mu\nu}$. We can decompose this perturbation into even and odd modes, denoted as $H_{\mu\nu}^{+}$ and $H_{\mu\nu}^{-}$, respectively. In the RW gauge and in spherical coordinates, these take the following form:

$$H_{\mu\nu}^{+} = \begin{pmatrix} \bar{H}_0 & \bar{H}_1 & 0 & 0 \\ \bar{H}_1 & \bar{H}_2 & 0 & 0 \\ 0 & 0 & r^2 \bar{K} & 0 \\ 0 & 0 & 0 & r^2 \sin^2 \theta \bar{K} \end{pmatrix}, \qquad (C.1)$$

$$H_{\mu\nu}^{-} = \begin{pmatrix} 0 & 0 & -\frac{\partial_{\phi}\bar{h}_0}{\sin\theta} & \sin\theta\partial_{\theta}\bar{h}_0 \\ 0 & 0 & -\frac{\partial_{\phi}\bar{h}_1}{\sin\theta} & \sin\theta\partial_{\theta}\bar{h}_1 \\ -\frac{\partial_{\phi}\bar{h}_0}{\sin\theta} & \sin\theta\partial_{\theta}\bar{h}_1 & 0 & 0 \\ -\frac{\partial_{\phi}\bar{h}_0}{\sin\theta} & \sin\theta\partial_{\theta}\bar{h}_1 & 0 & 0 \end{pmatrix}. \qquad (C.2)$$

To rewrite them in covariant form, we remind the reader that $v^{\mu}$ is the four-velocity of the massive object, which defines our time-like direction. In its rest frame, where we work, $v^{\mu} = (1, 0, 0, 0)$. The radial direction is defined by the spatial vector $r^{\mu} \equiv (\delta^{\mu}_{\nu} + v^{\mu}v_{\nu})x^{\nu}$, or by the unit vector parallel to it, $\hat{r}^{\mu} = r^{\mu}/\sqrt{r^{\rho}r_{\rho}} = (0, 1, 0, 0)$, where the indices are raised and lowered with the Minkowski metric. Finally, projecting orthogonally to these two directions, we obtain the metric on the unit 2-sphere, i.e.,

$$\Omega_{\mu\nu} = \eta_{\mu\nu} + v_{\mu}v_{\nu} - \hat{r}_{\mu}\hat{r}_{\nu}. \qquad (C.3)$$

Therefore, the even part of the tidal field simply reads

$$H_{\mu\nu}^{+} = \bar{H}_{0}v_{\mu}v_{\nu} + 2\bar{H}_{1}v_{(\mu}\hat{r}_{\nu)} + \bar{K}\Omega_{\mu\nu} \\ = (\bar{H}_{0} + \bar{K})v_{\mu}v_{\nu} + 2\bar{H}_{1}v_{(\mu}\hat{r}_{\nu)} + (\bar{H}_{2} - \bar{K})\hat{r}_{\mu}\hat{r}_{\nu} + \bar{K}\eta_{\mu\nu}, \qquad (C.4)$$

where for the last equality we have used eq. (C.3).

To obtain the odd part of the tidal field, it is convenient to introduce the two basis vectors on the unit 2-sphere

$$e_A^\mu = \frac{\partial \hat{r}^\mu}{\partial x^A}, \qquad x^A \in \{\theta, \phi\}, \qquad (C.5)$$

which are obviously orthogonal among each other and to the other two vectors, and satisfy

$$\eta_{\mu\nu} e_A^\mu e_B^\nu = \Omega_{AB}. \qquad (C.6)$$

Since by eq. (C.2) only $H_{0A}^{-} = v^{\mu}e_{\nu}^{\mu}H_{\mu\nu}^{-}$ and $H_{rA}^{-} = \hat{r}^{\mu}e_{\nu}^{\mu}H_{\mu\nu}^{-}$ are nonzero, we can write

$$H_{\mu\nu}^{-} = -2v_{(\mu}e_{\nu)}^{A}H_{0A}^{-} + 2\hat{r}_{(\mu}e_{\nu)}^{A}H_{rA}^{-} = 2v_{(\mu}e_{\nu)}^{A}\epsilon_{AB}\Omega^{BC}\partial_{C}\bar{h}_{0}/-2\hat{r}_{(\mu}e_{\nu)}^{A}\epsilon_{AB}\Omega^{BC}\partial_{C}\bar{h}_{1}. \qquad (C.7)$$

Replacing $\epsilon_{AB}$ on the right-hand side using the identity $\epsilon_{AB} = v^{\alpha}\hat{r}^{\beta}e_{A}^{\rho}e_{B}^{\sigma}\epsilon_{\alpha\beta\rho\sigma}$, we finally obtain

$$H_{\mu\nu}^{-} = -2v^{\alpha}\epsilon_{\alpha\rho\sigma}(\mu v_{\nu})\hat{r}^{\rho}\partial^{\sigma}\bar{h}_{0} + 2v^{\alpha}\epsilon_{\alpha\rho\sigma}(\mu\hat{r}_{\nu})\hat{r}^{\rho}\partial^{\sigma}\bar{h}_{1}. \qquad (C.8)$$
D Symmetric trace-free tensor relations

Here we derive some useful relations involving the STF tensors $\mathcal{Y}_{i_1...i_\ell}^m$ introduced in section 5.1.
We recall their definition,

$$ Y_{\ell}^{m}(\theta, \phi) = \mathcal{Y}_{i_1 \dots i_{\ell}}^{m} \hat{r}^{i_1} \dots \hat{r}^{i_{\ell}} = \mathcal{Y}_{i_L}^{m} \hat{r}^{i_L}, \qquad (\text{D.1}) $$

where in the second equality we use the multi-index notation introduced in section 4.1.

Using this, we rewrite the normalization relation of two spherical harmonics with the same $\ell$ in terms of the STF tensors,

$$ \mathcal{Y}_{i_L}^{m*} \mathcal{Y}_{j_L}^{m'} \int d\Omega \hat{r}^{i_2 L} = \delta^{mm'}. \qquad (\text{D.2}) $$

Solving the angular integral,<fnref n="25" /> we obtain

$$ \mathcal{Y}_{i_L}^{m*} \mathcal{Y}_{j_L}^{m'} (\delta_{i_1 i_2} \cdots \delta_{i_{2\ell-1} i_{2\ell}} + \text{perms}) = \frac{(2\ell+1)!!}{4\pi} \delta^{mm'}, \qquad (\text{D.3}) $$

from which we can derive

$$ \mathcal{Y}_{ijkl}^{m*} \mathcal{Y}_{ijkl}^{m'} = \frac{315}{32\pi} \delta^{mm'}, \qquad (\text{D.4}) $$

which is needed for the computation in section 5.2.

We can follow a similar procedure for the products of three STF tensors. In particular, rewriting the definition of $^0I_{\ell\ell_1\ell_2}^{m_1m_2}$ (eq. (2.26)) in terms of the STF tensors and solving the angular integral for $\ell + \ell_1 + \ell_2$ even, we obtain

$$ \mathcal{Y}_{i_L}^{m*} \mathcal{Y}_{j_L}^{m_1} \mathcal{Y}_{k_L}^{m_2} (\delta_{i_1 i_2} \cdots \delta_{i_{2\ell-1} i_{2\ell}} + \text{perms}) = \frac{(\ell + \ell_1 + \ell_2 + 1)!!}{4\pi} I_{\ell\ell_1\ell_2}^{m_1m_2}. \qquad (\text{D.5}) $$

For the computation in section 5.2 we need the following contractions

$$ \mathcal{Y}^0 \mathcal{Y}_{ij}^{m_1} \mathcal{Y}_{ij}^{m_2} = \frac{15}{8\pi} I_{022}^{0m_1m_2}, \qquad (\text{D.6}) $$

$$ \mathcal{Y}_{ij}^{m*} \mathcal{Y}_{jk}^{m_1} \mathcal{Y}_{ki}^{m_2} = \frac{105}{32\pi} I_{222}^{mm_1m_2}, \qquad (\text{D.7}) $$

$$ \mathcal{Y}_{ijkl}^{m*} \mathcal{Y}_{ij}^{m_1} \mathcal{Y}_{kl}^{m_2} = \frac{315}{32\pi} I_{422}^{mm_1m_2}. \qquad (\text{D.8}) $$

Finally, in order to simplify the final expression for the one-point function computed in section 5.2, one needs to decompose the following tensor $\mathcal{Y}_{ij}^{m_1} \mathcal{Y}_{kl}^{m_2}$. For symmetry reasons, one necessarily has

$$ \mathcal{Y}_{ij}^{m_1} \mathcal{Y}_{kl}^{m_2} = c_1 \delta_{ij} \delta_{kl} + 2c_2 \delta_{i(k} \delta_{l)j} + c_3 \mathcal{Y}_{ij}^{m} \delta_{kl} + c_4 \mathcal{Y}_{kl}^{m} \delta_{ij} + 2c_5 (\mathcal{Y}_{i(k}^{m} \delta_{l)j} + \mathcal{Y}_{j(k}^{m} \delta_{l)i}) + c_6 \mathcal{Y}_{ijkl}^{m}. \qquad (\text{D.9}) $$

The coefficients $c_n$ can now be found by contracting the above expression with all possible combinations of $\delta_{ij}$, $\mathcal{Y}_{ij}^{m*}$ and $\mathcal{Y}_{ijkl}^{m*}$. One eventually obtains

$$ \mathcal{Y}_{ij}^{m_1} \mathcal{Y}_{kl}^{m_2} = \frac{I_{022}^{0m_1m_2}}{4\sqrt{\pi}} (3\delta_{i(k} \delta_{l)j} - \delta_{ij} \delta_{kl}) - I_{222}^{mm_1m_2} \left( \mathcal{Y}_{ij}^{m} \delta_{kl} + \mathcal{Y}_{kl}^{m} \delta_{ij} - \frac{3}{2} \mathcal{Y}_{i(k}^{m} \delta_{l)j} - \frac{3}{2} \mathcal{Y}_{j(k}^{m} \delta_{l)i} \right) + I_{422}^{mm_1m_2} \mathcal{Y}_{ijkl}^{m}. \qquad (\text{D.10}) $$

<footnote n="25">²⁵We use $\int d\Omega \hat{r}^{i_1} \hat{r}^{i_2} \dots \hat{r}^{i_{2\ell-1}} \hat{r}^{i_{2\ell}} = \frac{4\pi}{(2\ell+1)!!} (\delta_{i_1 i_2} \dots \delta_{i_{2\ell-1} i_{2\ell}} + \text{perms})$.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3. Feynman diagrams that reconstruct the Schwarzschild metric up to order $r_s^2$.</figcaption></figure>

## E Recovering the Schwarzschild metric

In this appendix we show how one can consistently reconstruct, from the EFT worldline approach, the Schwarzschild metric in de Donder coordinates. We perform the explicit calculation up to $O(r_s^2)$, in which case the metric takes the form (see e.g. [96])

$$ds^2 = - \left(1 - \frac{r_s}{r} + \frac{1}{2} \frac{r_s^2}{r^2}\right) dt^2 + \left(1 + \frac{r_s}{r} + \frac{5}{4} \frac{r_s^2}{r^2}\right) d\mathbf{x} \cdot d\mathbf{x} + \frac{7}{4} \frac{r_s^2}{r^2} \frac{(\mathbf{x} \cdot d\mathbf{x})^2}{r^2} + O(r_s^3). \quad (E.1)$$

See also [55, 96, 97] for previous results and [98, 99] for a recent reconstruction of the Schwarzschild metric to all orders. For the calculation we will assume dimensional regularisation in $d = D + 1$ dimensions.

Since for this computation $H_{\mu\nu} = 0$, we only need the two diagrams in figure 3. At leading order we only have the point-particle sourcing $h_{\mu\nu}$, i.e. the diagram 3(a). This gives

$$ \begin{aligned} \kappa^{(0)}h_{\mu\nu}(x) &= -8\pi r_s \int dt v^\rho v^\sigma \int_k \frac{e^{-ik \cdot (x(t)-x)}}{k^2} P_{\rho\sigma\mu\nu} + \mathcal{O}((r_s/r)^2) \\ &= 8\pi r_s \left(v_\mu v_\nu + \frac{1}{2}\eta_{\mu\nu}\right) \int_k \frac{e^{ik \cdot x}}{k^2} + \mathcal{O}((r_s/r)^2) \\ &= \frac{r_s}{r} (\eta_{\mu\nu} + 2v_\mu v_\nu) + \mathcal{O}((r_s/r)^2), \end{aligned} \quad (E.2) $$

where $r = \sqrt{\mathbf{x}^2} = \sqrt{\delta_{ij}x^i x^j}$ and we used $\int_k \frac{e^{ik \cdot x}}{(k^2)^\alpha} = \frac{1}{(4\pi)^{D/2}} \frac{\Gamma(\frac{D}{2}-\alpha)}{\Gamma(\alpha)} \left(\frac{\mathbf{x}^2}{4}\right)^{\alpha-\frac{D}{2}}$ for the last equality. One can check that $g_{\mu\nu} = \eta_{\mu\nu} + \kappa^{(0)}h_{\mu\nu}(x)$ reproduces the Schwarzschild metric in de Donder gauge (E.1) expanded up to $O(r_s)$.

At next-to-leading order we have to compute the diagram 3(b). Using the rules defined in the previous section, we explicitly get

$$ \begin{aligned} \kappa^{(0)} h_{\mu\nu}(x)_{1b} &= -32\pi^2 r_s^2 \int dt_1 dt_2 v^{\rho_1} v^{\sigma_1} v^{\rho_2} v^{\sigma_2} \int_{q,k} \frac{e^{-iq \cdot (x(t_1)-x(t_2)) - ik \cdot (x(t_1)-x)}}{q^2 k^2 (k-q)^2} \\ &P_{\rho_1 \sigma_1 \alpha_1 \beta_1} P_{\rho_2 \sigma_2 \alpha_2 \beta_2} V^{\alpha_1 \beta_1 \alpha_2 \beta_2 \alpha_3 \beta_3} P_{\alpha_3 \beta_3 \mu\nu} \\ &= 2\pi^2 r_s^2 \int_k \delta^{-}(k^0) \frac{e^{ik \cdot x}}{k^2} (7k_\mu k_\nu + 3k^2 \eta_{\mu\nu} - k^2 v_\mu v_\nu) \int_q \frac{\delta^{-}(q^0)}{q^2 (k-q)^2}, \end{aligned} \quad (E.3) $$

where the second equality follows from integrating in $t_1$ and $t_2$ and from tensorial reduction
of the integrands. Solving the integrals in the last line<sup>26</sup> and adding this contribution to the leading order one, we finally obtain

$$ \kappa^{(0)}h_{\mu\nu}(x) = \frac{r_s}{r} (\eta_{\mu\nu} + 2v_\mu v_\nu) + \frac{r_s^2}{4r^2} (5\eta_{\mu\nu} + 3v_\mu v_\nu - 7\hat{r}_\mu \hat{r}_\nu) + \mathcal{O}((r_s/r)^3), \quad (\text{E.4}) $$

where we recall that the normalised vector in the radial direction is $\hat{r}^\mu = (0, x^i/r)$. This result reproduces the Schwarzschild metric in de Donder gauge (E.1) expanded up to $\mathcal{O}(r_s^2)$ and agrees with refs. [96, 97].

**Data Availability Statement.** This article has no associated data or the data will not be deposited.

**Code Availability Statement.** This article has no associated code or the code will not be deposited.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

## References

[1] KAGRA et al. collaborations, GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run, Phys. Rev. X 13 (2023) 041039 [arXiv:2111.03606] [INSPIRE].

[2] L. Barack et al., Black holes, gravitational waves and fundamental physics: a roadmap, Class. Quant. Grav. 36 (2019) 143001 [arXiv:1806.05195] [INSPIRE].

[3] B.S. Sathyaprakash et al., Extreme Gravity and Fundamental Physics, arXiv:1903.09221 [INSPIRE].

[4] ET collaboration, Science Case for the Einstein Telescope, JCAP 03 (2020) 050 [arXiv:1912.02622] [INSPIRE].

[5] E. Barausse et al., Prospects for Fundamental Physics with LISA, Gen. Rel. Grav. 52 (2020) 81 [arXiv:2001.09793] [INSPIRE].

[6] V. Kalogera et al., The Next Generation Global Gravitational Wave Observatory: The Science Book, arXiv:2111.06990 [INSPIRE].

[7] E. Berti et al., Snowmass2021 Cosmic Frontier White Paper: Fundamental Physics and Beyond the Standard Model, in the proceedings of the Snowmass 2021, Seattle, U.S.A., July 17–26 (2022) [arXiv:2203.06240] [INSPIRE].

$$ \int_k \frac{1}{[\mathbf{q}^2]^a [(\mathbf{k}-\mathbf{q})^2]^b} = \frac{(\mathbf{q}^2)^{\frac{D}{2}-a-b}}{(4\pi)^{\frac{D}{2}}} \frac{\Gamma(a+b-\frac{D}{2})}{\Gamma(a)\Gamma(b)} \frac{\Gamma(\frac{D}{2}-a)\Gamma(\frac{D}{2}-b)}{\Gamma(D-a-b)}, \\ \int_k \frac{k^a k^b k^c}{(\mathbf{k}^2)^\alpha} e^{ik \cdot x} = \frac{i}{(4\pi)^{D/2}} \frac{\Gamma(\frac{D}{2}-\alpha+2)}{\Gamma(\alpha)} \left(\frac{x^2}{4}\right)^{\alpha-\frac{D}{2}-2} \left\{ \left(\alpha-\frac{D}{2}-2\right) \frac{x^a x^b x^c}{2x^2} + \frac{\delta^{ab} x^c}{4} + \frac{\delta^{c(a} x^{b)}}{2} \right\}. $$

<footnote n="26"><sup>26</sup>These two integrals are useful</footnote>
[8] L. Blanchet, *Post-Newtonian Theory for Gravitational Waves*, Living Rev. Rel. 17 (2014) 2 [arXiv:1310.1528] [INSPIRE].

[9] A. Buonanno et al., *Snowmass White Paper: Gravitational Waves and Scattering Amplitudes*, in the proceedings of the *Snowmass 2021*, Seattle, U.S.A., July 17–26 (2022) [arXiv:2204.05194] [INSPIRE].

[10] W.D. Goldberger, *Effective field theories of gravity and compact binary dynamics: A Snowmass 2021 whitepaper*, in the proceedings of the *Snowmass 2021*, Seattle, U.S.A., July 17–26 (2022) [arXiv:2206.14249] [INSPIRE].

[11] E.E. Flanagan and T. Hinderer, *Constraining neutron star tidal Love numbers with gravitational wave detectors*, Phys. Rev. D 77 (2008) 021502 [arXiv:0709.1915] [INSPIRE].

[12] G. Kälin, Z. Liu and R.A. Porto, *Conservative Tidal Effects in Compact Binary Systems to Next-to-Leading Post-Minkowskian Order*, Phys. Rev. D 102 (2020) 124025 [arXiv:2008.06047] [INSPIRE].

[13] Q. Henry, G. Faye and L. Blanchet, *Tidal effects in the gravitational-wave phase evolution of compact binary systems to next-to-next-to-leading post-Newtonian order*, Phys. Rev. D 102 (2020) 044033 [Erratum ibid. 108 (2023) 089901] [arXiv:2005.13367] [INSPIRE].

[14] H. Fang and G. Lovelace, *Tidal coupling of a Schwarzschild black hole and circularly orbiting moon*, Phys. Rev. D 72 (2005) 124016 [gr-qc/0505156] [INSPIRE].

[15] T. Damour and A. Nagar, *Relativistic tidal properties of neutron stars*, Phys. Rev. D 80 (2009) 084035 [arXiv:0906.0096] [INSPIRE].

[16] T. Binnington and E. Poisson, *Relativistic theory of tidal Love numbers*, Phys. Rev. D 80 (2009) 084018 [arXiv:0906.1366] [INSPIRE].

[17] J. Vines, E.E. Flanagan and T. Hinderer, *Post-1-Newtonian tidal effects in the gravitational waveform from binary inspirals*, Phys. Rev. D 83 (2011) 084051 [arXiv:1101.1673] [INSPIRE].

[18] D. Bini, T. Damour and G. Faye, *Effective action approach to higher-order relativistic tidal interactions in binary systems and their effective one body description*, Phys. Rev. D 85 (2012) 124034 [arXiv:1202.3565] [INSPIRE].

[19] L. Baiotti and L. Rezzolla, *Binary neutron star mergers: a review of Einstein's richest laboratory*, Rept. Prog. Phys. 80 (2017) 096901 [arXiv:1607.03540] [INSPIRE].

[20] B. Kol and M. Smolkin, *Black hole stereotyping: Induced gravito-static polarization*, JHEP 02 (2012) 010 [arXiv:1110.3764] [INSPIRE].

[21] N. Gürlebeck, *No-hair theorem for Black Holes in Astrophysical Environments*, Phys. Rev. Lett. 114 (2015) 151102 [arXiv:1503.03240] [INSPIRE].

[22] L. Hui et al., *Static response and Love numbers of Schwarzschild black holes*, JCAP 04 (2021) 052 [arXiv:2010.00593] [INSPIRE].

[23] L. Hui et al., *Ladder symmetries of black holes. Implications for love numbers and no-hair theorems*, JCAP 01 (2022) 032 [arXiv:2105.01069] [INSPIRE].

[24] M. Rai and L. Santoni, *Ladder symmetries and Love numbers of Reissner-Nordström black holes*, JHEP 07 (2024) 098 [arXiv:2404.06544] [INSPIRE].

[25] A. Le Tiec and M. Casals, *Spinning Black Holes Fall in Love*, Phys. Rev. Lett. 126 (2021) 131102 [arXiv:2007.00214] [INSPIRE].

[26] A. Le Tiec, M. Casals and E. Franzin, *Tidal Love Numbers of Kerr Black Holes*, Phys. Rev. D 103 (2021) 084021 [arXiv:2010.15795] [INSPIRE].
[27] P. Charalambous, S. Dubovsky and M.M. Ivanov, On the Vanishing of Love Numbers for Kerr Black Holes, JHEP 05 (2021) 038 [arXiv:2102.08917] [INSPIRE].

[28] M.J. Rodriguez, L. Santoni, A.R. Solomon and L.F. Temoche, Love numbers for rotating black holes in higher dimensions, Phys. Rev. D 108 (2023) 084011 [arXiv:2304.03743] [INSPIRE].

[29] E. Franzin, V. Cardoso, P. Pani and G. Raposo, Testing strong gravity with gravitational waves and Love numbers, J. Phys. Conf. Ser. 841 (2017) 012035 [INSPIRE].

[30] V. Cardoso et al., Testing strong-field gravity with tidal Love numbers, Phys. Rev. D 95 (2017) 084014 [Addendum ibid. 95 (2017) 089901] [arXiv:1701.01116] [INSPIRE].

[31] P. Pani and A. Maselli, Love in Extrema Ratio, Int. J. Mod. Phys. D 28 (2019) 1944001 [arXiv:1905.03947] [INSPIRE].

[32] V. Cardoso and P. Pani, Testing the nature of dark compact objects: a status report, Living Rev. Rel. 22 (2019) 4 [arXiv:1904.05363] [INSPIRE].

[33] H.S. Chia et al., In pursuit of Love numbers: First templated search for compact objects with large tidal deformabilities in the LIGO-Virgo data, Phys. Rev. D 110 (2024) 063007 [arXiv:2306.00050] [INSPIRE].

[34] H.S. Chia, Z. Zhou and M.M. Ivanov, Bring the Heat: Tidal Heating Constraints for Black Holes and Exotic Compact Objects from the LIGO-Virgo-KAGRA Data, arXiv:2404.14641 [INSPIRE].

[35] E. Poisson and M. Sasaki, Gravitational radiation from a particle in circular orbit around a black hole. 5: Black hole absorption and tail corrections, Phys. Rev. D 51 (1995) 5753 [gr-qc/9412027] [INSPIRE].

[36] H. Tagoshi, S. Mano and E. Takasugi, PostNewtonian expansion of gravitational waves from a particle in circular orbits around a rotating black hole: Effects of black hole absorption, Prog. Theor. Phys. 98 (1997) 829 [gr-qc/9711072] [INSPIRE].

[37] T. Damour, Gravitational radiation and the motion of compact bodies, in the proceedings of the Les Houches Summer School on Gravitational Radiation, Les Houches, France, June 02–21 (1982) [INSPIRE].

[38] R.A. Porto, The effective field theorist's approach to gravitational dynamics, Phys. Rept. 633 (2016) 1 [arXiv:1601.04914] [INSPIRE].

[39] I.Z. Rothstein, Progress in effective field theory approach to the binary inspiral problem, Gen. Rel. Grav. 46 (2014) 1726 [INSPIRE].

[40] R.A. Porto, The Tune of Love and the Nature(ness) of Spacetime, Fortsch. Phys. 64 (2016) 723 [arXiv:1606.08895] [INSPIRE].

[41] L. Hui et al., Near-zone symmetries of Kerr black holes, JHEP 09 (2022) 049 [arXiv:2203.08832] [INSPIRE].

[42] J. Ben Achour, E.R. Livine, S. Mukohyama and J.-P. Uzan, Hidden symmetry of the static response of black holes: applications to Love numbers, JHEP 07 (2022) 112 [arXiv:2202.12828] [INSPIRE].

[43] R. Berens, L. Hui and Z. Sun, Ladder symmetries of black holes and de Sitter space: love numbers and quasinormal modes, JCAP 06 (2023) 056 [arXiv:2212.09367] [INSPIRE].

[44] P. Charalambous, S. Dubovsky and M.M. Ivanov, Hidden Symmetry of Vanishing Love Numbers, Phys. Rev. Lett. 127 (2021) 101101 [arXiv:2103.01234] [INSPIRE].

[45] P. Charalambous, S. Dubovsky and M.M. Ivanov, Love symmetry, JHEP 10 (2022) 175 [arXiv:2209.02091] [INSPIRE].
[46] S. Chakrabarti, T. Delsate and J. Steinhoff, *New perspectives on neutron star and black hole spectroscopy and dynamic tides*, arXiv:1304.2228 [INSPIRE].

[47] E. Poisson, *Compact body in a tidal environment: New types of relativistic Love numbers, and a post-Newtonian operational definition for tidally induced multipole moments*, Phys. Rev. D 103 (2021) 064023 [arXiv:2012.10184] [INSPIRE].

[48] M.V.S. Saketh, Z. Zhou and M.M. Ivanov, *Dynamical tidal response of Kerr black holes from scattering amplitudes*, Phys. Rev. D 109 (2024) 064058 [arXiv:2307.10391] [INSPIRE].

[49] M. Perry and M.J. Rodriguez, *Dynamical Love Numbers for Kerr Black Holes*, arXiv:2310.03660 [INSPIRE].

[50] G.S. He, *Nonlinear Optics and Photonics*, Oxford University Press (2014) [DOI:10.1093/acprof:oso/9780198702764.001.0001].

[51] E. Poisson, *Tidally induced multipole moments of a nonrotating black hole vanish to all post-Newtonian orders*, Phys. Rev. D 104 (2021) 104062 [arXiv:2108.07328] [INSPIRE].

[52] M.M. Riva, L. Santoni, N. Savić and F. Vernizzi, *Vanishing of nonlinear tidal Love numbers of Schwarzschild black holes*, Phys. Lett. B 854 (2024) 138710 [arXiv:2312.05065] [INSPIRE].

[53] E. Poisson and I. Vlasov, *Geometry and dynamics of a tidally deformed black hole*, Phys. Rev. D 81 (2010) 024029 [arXiv:0910.4311] [INSPIRE].

[54] V. De Luca, J. Khoury and S.S.C. Wong, *Nonlinearities in the tidal Love numbers of black holes*, Phys. Rev. D 108 (2023) 024048 [arXiv:2305.14444] [INSPIRE].

[55] W.D. Goldberger and I.Z. Rothstein, *An effective field theory of gravity for extended objects*, Phys. Rev. D 73 (2006) 104029 [hep-th/0409156] [INSPIRE].

[56] W.D. Goldberger and I.Z. Rothstein, *Dissipative effects in the worldline approach to black hole dynamics*, Phys. Rev. D 73 (2006) 104030 [hep-th/0511133] [INSPIRE].

[57] W.D. Goldberger and I.Z. Rothstein, *Towers of Gravitational Theories*, Gen. Rel. Grav. 38 (2006) 1537 [hep-th/0605238] [INSPIRE].

[58] M. Levi, *Effective Field Theories of Post-Newtonian Gravity: A comprehensive review*, Rept. Prog. Phys. 83 (2020) 075901 [arXiv:1807.01699] [INSPIRE].

[59] W.D. Goldberger, *Effective Field Theory for Compact Binary Dynamics*, arXiv:2212.06677 [DOI:10.1007/978-981-19-3079-9_2-1] [INSPIRE].

[60] T. Regge and J.A. Wheeler, *Stability of a Schwarzschild singularity*, Phys. Rev. 108 (1957) 1063 [INSPIRE].

[61] C.O. Nicasio, R.J. Gleiser, R.H. Price and J. Pullin, *The collision of boosted black holes: Second order close limit calculations*, Phys. Rev. D 59 (1999) 044024 [gr-qc/9802063] [INSPIRE].

[62] R.J. Gleiser, C.O. Nicasio, R.H. Price and J. Pullin, *Gravitational radiation from Schwarzschild black holes: The second order perturbation formalism*, Phys. Rept. 325 (2000) 41 [gr-qc/9807077] [INSPIRE].

[63] D. Brizuela, J.M. Martin-Garcia and G.A. Mena Marugan, *Second and higher-order perturbations of a spherical spacetime*, Phys. Rev. D 74 (2006) 044039 [gr-qc/0607025] [INSPIRE].

[64] D. Brizuela, J.M. Martin-Garcia and G.A.M. Marugan, *High-order gauge-invariant perturbations of a spherical spacetime*, Phys. Rev. D 76 (2007) 024004 [gr-qc/0703069] [INSPIRE].

[65] H. Nakano and K. Ioka, *Second Order Quasi-Normal Mode of the Schwarzschild Black Hole*, Phys. Rev. D 76 (2007) 084007 [arXiv:0708.0450] [INSPIRE].
[66] D. Brizuela, J.M. Martin-Garcia and M. Tiglio, *A Complete gauge-invariant formalism for arbitrary second-order perturbations of a Schwarzschild black hole*, Phys. Rev. D 80 (2009) 024021 [arXiv:0903.1134] [INSPIRE].

[67] M. Lagos and L. Hui, *Generation and propagation of nonlinear quasinormal modes of a Schwarzschild black hole*, Phys. Rev. D 107 (2023) 044040 [arXiv:2208.07379] [INSPIRE].

[68] P. Bourg et al., *Quadratic quasi-normal mode dependence on linear mode parity*, arXiv:2405.10270 [INSPIRE].

[69] B. Bucciotti, L. Juliano, A. Kuntz and E. Trincherini, *Quadratic quasinormal modes of a Schwarzschild black hole*, Phys. Rev. D 110 (2024) 104048 [arXiv:2405.06012] [INSPIRE].

[70] B. Bucciotti, L. Juliano, A. Kuntz and E. Trincherini, *Amplitudes and polarizations of quadratic quasi-normal modes for a Schwarzschild black hole*, JHEP 09 (2024) 119 [arXiv:2406.14611] [INSPIRE].

[71] H. Motohashi, T. Suyama and K. Takahashi, *Fundamental theorem on gauge fixing at the action level*, Phys. Rev. D 94 (2016) 124021 [arXiv:1608.00071] [INSPIRE].

[72] G. Franciolini et al., *Effective Field Theory of Black Hole Quasinormal Modes in Scalar-Tensor Theories*, JHEP 02 (2019) 127 [arXiv:1810.07706] [INSPIRE].

[73] L. Hui, A. Podo, L. Santoni and E. Trincherini, *Effective Field Theory for the perturbations of a slowly rotating black hole*, JHEP 12 (2021) 183 [arXiv:2111.02072] [INSPIRE].

[74] T. Kobayashi, H. Motohashi and T. Suyama, *Black hole perturbation in the most general scalar-tensor theory with second-order field equations II: the even-parity sector*, Phys. Rev. D 89 (2014) 084042 [arXiv:1402.6740] [INSPIRE].

[75] G. Franciolini et al., *Stable wormholes in scalar-tensor theories*, JHEP 01 (2019) 221 [arXiv:1811.05481] [INSPIRE].

[76] F.J. Zerilli, *Effective potential for even parity Regge-Wheeler gravitational perturbation equations*, Phys. Rev. Lett. 24 (1970) 737 [INSPIRE].

[77] K.S. Thorne and A. Campolattaro, *Non-Radial Pulsation of General-Relativistic Stellar Models. I. Analytic Analysis for $L \ge 2$*, Astrophys. J. 149 (1967) 591.

[78] T. Hinderer, *Tidal Love numbers of neutron stars*, Astrophys. J. 677 (2008) 1216 [Erratum ibid. 697 (2009) 964] [arXiv:0711.2420] [INSPIRE].

[79] H. Bateman and A. Erdélyi, *Higher transcendental functions*, McGraw-Hill, New York (1955) [https://cds.cern.ch/record/100233].

[80] B.S. DeWitt, *Quantum Theory of Gravity. II. The Manifestly Covariant Theory*, Phys. Rev. 162 (1967) 1195 [INSPIRE].

[81] G. 't Hooft and M.J.G. Veltman, *One loop divergencies in the theory of gravitation*, Ann. Inst. H. Poincare Phys. Theor. 20 (1974) 69 [INSPIRE].

[82] L.F. Abbott, *The Background Field Method Beyond One Loop*, Nucl. Phys. B 185 (1981) 189 [INSPIRE].

[83] W.D. Goldberger and A. Ross, *Gravitational radiative corrections from effective field theory*, Phys. Rev. D 81 (2010) 124015 [arXiv:0912.4254] [INSPIRE].

[84] Z. Bern et al., *Leading Nonlinear Tidal Effects and Scattering Amplitudes*, JHEP 05 (2021) 188 [arXiv:2010.08559] [INSPIRE].
[85] W.D. Goldberger, J. Li and I.Z. Rothstein, Non-conservative effects on spinning black holes from world-line effective field theory, JHEP 06 (2021) 053 [arXiv:2012.14869] [INSPIRE].

[86] M.M. Ivanov and Z. Zhou, Revisiting the matching of black hole tidal responses: A systematic study of relativistic and logarithmic corrections, Phys. Rev. D 107 (2023) 084030 [arXiv:2208.08459] [INSPIRE].

[87] T. Hadad, B. Kol and M. Smolkin, Gravito-magnetic polarization of Schwarzschild black hole, JHEP 06 (2024) 169 [arXiv:2402.16172] [INSPIRE].

[88] H.S. Chia, Tidal deformation and dissipation of rotating black holes, Phys. Rev. D 104 (2021) 024013 [arXiv:2010.07300] [INSPIRE].

[89] K.S. Thorne, Multipole Expansions of Gravitational Radiation, Rev. Mod. Phys. 52 (1980) 299 [INSPIRE].

[90] M. Maggiore, Gravitational Waves. Volume 1: Theory and Experiments, Oxford University Press (2007) [DOI:10.1093/acprof:oso/9780198570745.001.0001] [INSPIRE].

[91] G. Kristensson, Second Order Differential Equations: Special Functions and Their Classification, Springer New York (2010) [DOI:10.1007/978-1-4419-7020-6].

[92] O. Combaluzier-Szteinsznaider et al., Symmetries of Vanishing Nonlinear Love Numbers of Schwarzschild Black Holes, arXiv:2410.10952 [INSPIRE].

[93] A. Kehagias and A. Riotto, Black Holes in a Gravitational Field: The Non-linear Static Love Number of Schwarzschild Black Holes Vanishes, arXiv:2410.11014 [INSPIRE].

[94] P. Charalambous and M.M. Ivanov, Scalar Love numbers and Love symmetries of 5-dimensional Myers-Perry black holes, JHEP 07 (2023) 222 [arXiv:2303.16036] [INSPIRE].

[95] R.A. Rosen and L. Santoni, Black hole perturbations of massive and partially massless spin-2 fields in (anti) de Sitter spacetime, JHEP 03 (2021) 139 [arXiv:2010.00595] [INSPIRE].

[96] S. Mougiakakos and P. Vanhove, Schwarzschild-Tangherlini metric from scattering amplitudes in various dimensions, Phys. Rev. D 103 (2021) 026001 [arXiv:2010.08882] [INSPIRE].

[97] G.U. Jakobsen, Schwarzschild-Tangherlini Metric from Scattering Amplitudes, Phys. Rev. D 102 (2020) 104065 [arXiv:2006.01734] [INSPIRE].

[98] P.H. Damgaard and K. Lee, Schwarzschild Black Hole from Perturbation Theory to All Orders, Phys. Rev. Lett. 132 (2024) 251603 [arXiv:2403.13216] [INSPIRE].

[99] S. Mougiakakos and P. Vanhove, Schwarzschild Metric from Scattering Amplitudes to All Orders in GN, Phys. Rev. Lett. 133 (2024) 111601 [arXiv:2405.14421] [INSPIRE].