# Multipole expansion of gravitational waves: memory effects and Bondi aspects

Luc Blanchet,ª Geoffrey Compère,ᵇ Guillaume Faye,ª Roberto Oliveriᶜ and Ali Serajᵈ

E-mail: luc.blanchet@iap.fr, gcompere@ulb.be, faye@iap.fr,
roberto.oliveri@obspm.fr, a.seraj@qmul.ac.uk

**ABSTRACT:** In our previous work, we proposed an algorithm to transform the metric of an isolated matter source in the multipolar post-Minkowskian approximation in harmonic (de Donder) gauge to the Newman-Unti gauge. We then applied this algorithm at linear order and for specific quadratic interactions known as quadratic tail terms. In the present work, we extend this analysis to quadratic interactions associated with the coupling of two mass quadrupole moments, including both instantaneous and hereditary terms. Our main result is the derivation of the metric in Newman-Unti and Bondi gauges with complete quadrupole-quadrupole interactions. We rederive the displacement memory effect and provide expressions for all Bondi aspects and dressed Bondi aspects relevant to the study of leading and subleading memory effects. Then we obtain the Newman-Penrose charges, the BMS charges as well as the second and third order celestial charges defined from the known second order and novel third order dressed Bondi aspects for mass monopole-quadrupole and quadrupole-quadrupole interactions.

**KEYWORDS:** Classical Theories of Gravity, Space-Time Symmetries

ARXIV ePRINT: 2303.07732

<footnote>ªGReCO, Institut d'Astrophysique de Paris, UMR 7095, CNRS & Sorbonne Université,
98bis boulevard Arago, 75014 Paris, France</footnote>

<footnote>ᵇUniversité Libre de Bruxelles, International Solvay Institutes,
CP 231, B-1050 Brussels, Belgium</footnote>

<footnote>ᶜLUTH, Laboratoire Univers et Théories, Observatoire de Paris,
CNRS, Université PSL, Université Paris Cité,
5 place Jules Janssen, 92190 Meudon, France</footnote>

<footnote>ᵈSchool of Mathematical Sciences, Queen Mary University of London,
Mile End Road, E1 4NS, U.K.</footnote>
# Contents

1 **Introduction** 1

2 **Quadratic memory from the MPM formalism** 4

2.1 Original formulation: modified harmonic coordinates 4

2.2 Quadratic memory in Newman-Unti coordinates 9

2.3 $M_{ij} \times M_{ij}$ and $M \times M_{ij}$ asymptotic data in NU gauge 14

3 **Bondi aspects and charges** 18

3.1 Newman-Penrose charges 18

3.2 From Newman-Unti to Bondi coordinates 20

3.3 Flux-balance laws of the dressed $n \le 3$ Bondi aspects 22

3.4 The $n \le 3$ celestial charges at $G^2$ order 25

4 **Conclusions** 27

A **Treatment of hereditary integrals** 28

A.1 Radial integration at fixed retarded time 28

A.2 Far-zone expansion 31

B **Dressed Bondi aspects for $M_{ij} \times M_{ij}$ interactions** 32

## 1 Introduction

The displacement memory effect, i.e., the permanent change in the wave amplitude from before to after the passage of a gravitational wave strain, is a definite prediction of general relativity in the non-linear regime [1–5]. The effect was initially computed using post-Newtonian (PN) theory in [1], pointing out that it is dominantly of order 2.5PN in contrast to tails arising at 1.5PN order. This was published later, together with physical interpretation, in [5]. The effect can be interpreted as due to a finite accumulation of the electric (or mass type) part of the shear during the emission of gravitational radiation, which effectively induces a transition between the initial and final asymptotic Bondi-van der Burg-Metzner-Sachs (BMS) [6, 7] frames related by a supertranslation [8, 9]. The amount of displacement memory caused by a compact binary has been computed analytically in the inspiral phase [3, 5, 10] using post-Newtonian/post-Minkowskian (PN/PM) methods [11–13] and, more recently, in the merger phase using numerical relativity with adapted asymptotic BMS frames [14, 15]. Quite remarkably, it could be observed in the coming years with either ground-based gravitational-wave detectors [16–18] or pulsar timing arrays [19, 20].<fnref n="1" />

<footnote n="1">¹Besides the non-linear memory effect of interest in this paper, the displacement memory is also caused at linear order in $G$ by a change of Bondi mass aspect (e.g., a change between initial and final compact</footnote>
Additional classes of memory effects, i.e., gravitational-wave observables that measure persistent changes between initial and final non-radiative states, have been constructed [27–35] and their relationship to extended symmetry structures of asymptotically flat spacetimes [36–43] has been partially established [44–46]. Such memory effects are subleading in magnitude but could be potentially observed from space-borne gravitational-wave detectors [28].

Memory effects, especially the subleading ones, are most easily described in radiative gauges (see [47–49] for a general definition), such as Bondi gauge [6, 7] or closely related Newman-Unti (NU) gauge [50, 51]. In this gauge-fixed formulation, the formal asymptotic Einstein solution can be obtained systematically to derive all flux-balance laws resulting from Einstein’s constraint equations in terms of an asymptotic expansion [31, 37, 52–56]. Such flux-balance laws precisely encode the memory effects and their associated conserved charges [14, 31, 57–60].

In our joint recent work [61], building upon the formulation of general radiative gauges with polynomial expansions [12], we constructed the algorithm to transform the metric of an isolated matter source in the multipolar post-Minkowskian (MPM) approximation in harmonic/de Donder gauge to NU gauge and we applied this algorithm at linear order and for specific quadratic interactions known as quadratic tail terms [13]. The main technical objective of the present work is to derive the NU metric of such isolated matter source for the quadratic interaction associated with the coupling of two mass quadrupole moments. It is an intermediate milestone for achieving the larger aim of deriving the NU metric that entirely encodes the 2.5PN radiation from binary compact mergers. This metric only requires tails, the quadrupole-quadrupole interactions which will be investigated here, and spin-quadrupole interactions [62–64].

The quadrupole-quadrupole interaction provides the leading contribution to displacement memory, so this will allow us to discuss this effect as well as subleading memories within the NU metric setup for this interaction. In addition, we shall derive (starting from the result in harmonic coordinates) the complete quadrupole-quadrupole NU metric including the memory and all instantaneous (i.e., “non-hereditary”) terms. In particular, we shall explicitly verify the cancellation of all far-zone logarithms associated with harmonic coordinates, which confirms the soundness of the general algorithm proposed in [61].

The content of the paper is as follows. In section 2, we review the presentation of hereditary terms in the MPM formalism following [5] and we discuss the Poincaré flux-balance laws using the updated results of [60, 65, 66]. We then apply the algorithm developed in [61] to memory and mass-loss hereditary terms arising at quadratic order (the remaining tail hereditary terms were treated previously in [61]), specializing to mass quadrupole-quadrupole interactions only when presenting explicit expressions. The result is valid to any order in the distance to the source. We find the corresponding Newman-Unti metric and Bondi data. We proceed with rederiving the Christodoulou formula for the displacement memory [2] using the energy flux-balance law. We finally compute the

<footnote>object velocities) [21, 22] and by matter radiation reaching null infinity [9, 23–25]. These effects are usually called ordinary and null memory, respectively [26].</footnote>
full quadratic NU metric corresponding to mass monopole-quadrupole and quadrupole-quadrupole interactions. In section 3, we discuss various aspects of gravitational charges in the multipolar expansion. We first extend to the radiative case (but for the specific multipolar interactions considered) the classic result [53, 67] establishing that the Newman-Penrose charges [67, 68] can be expressed in the center-of-mass frame in terms of the product of the ADM mass and the initial mass quadrupole moment. Second, we extend the dictionary between NU gauge and Bondi gauge initiated in [61] to quadratic order in $G$. We revise the flux-balance laws, the BMS and conserved celestial charges following recent developments [28, 31, 40, 41, 45, 58–60, 69] and provide a corrected definition of the dressed $n=3$ Bondi aspect alternative to [31], which was subsequently revised in the erratum of [31]. We proceed by explicitly computing all lower order $n=0, 1, 2, 3$ Bondi charges in the presence of mass monopole-quadrupole and quadrupole-quadrupole interactions. We conclude in section 4. Appendix A is devoted to the technical treatment of hereditary integrals. Appendix B gives more information on dressed Bondi aspects for quadratic interactions.

**Notation and conventions.** We adopt units in which the speed of light $c$ is set to 1. The Newton gravitational constant $G$ is kept explicit to bookmark post-Minkowskian (PM) orders. Lower case Latin indices from $a$ to $h$ will refer to indices on the two-dimensional sphere, while lower case Latin indices from $i$ to $z$ will refer to three-dimensional Cartesian indices. The Minkowski metric is $\eta_{\mu\nu} = \text{diag}(-1, +1, +1, +1)$.

We denote Cartesian coordinates as $x^\mu = (t, \mathbf{x})$ and spherical ones as $(t, r, \theta^a)$. More precisely, the radial coordinate $r$ is defined as $r = |\mathbf{x}|$, while $\theta^a = (\theta, \varphi)$ with $a, b, \dots = \{1, 2\}$. The unit directional vector is denoted as $n^i = n^i(\theta^a) = x^i/r$. Euclidean spatial indices $i, j, \dots = \{1, 2, 3\}$ are raised and lowered with the Kronecker metric $\delta_{ij}$.

Furthermore, $L = i_1 i_2 \cdots i_\ell$ represents a multi-index made of $\ell$ spatial indices. We use short-hands for: the multi-derivative operator $\partial_L = \partial_{i_1} \cdots \partial_{i_\ell}$ with $\partial_i = \partial/\partial x^i$, the product of vectors $n_L = n_{i_1} \cdots n_{i_\ell}$, as well as $x_L = x_{i_1} \cdots x_{i_\ell} = r^\ell n_L$. The multipole moments $M_L$ and $S_L$ are symmetric-trace-free (STF) tensors over $L$. Time derivatives are indicated by superscripts $(q)$ or by dots.

We define the Minkowskian outgoing vector $k^\mu \partial_\mu = \partial_t + n^i \partial_i$, or, in components, $k^\mu = (1, n^i)$ and $k_\mu = (-1, n^i)$. In retarded spherical coordinates $(u, r, \theta^a)$ with $u = t - r$, we have $k^\mu \partial_\mu = \partial_r|_u$. We commonly employ the natural basis on the unit 2-sphere embedded in $\mathbb{R}^3$, $e_a = \frac{\partial}{\partial \theta^a}$, with components $e_a^i = \partial n^i / \partial \theta^a$. Given the metric on the unit sphere $\gamma_{ab} = \text{diag}(1, \sin^2 \theta)$, we have: $n^i e_a^i = 0$, $\partial_i \theta^a = r^{-1} \gamma^{ab} e_b^i$, $\gamma_{ab} = \delta_{ij} e_a^i e_b^j$ and $\gamma^{ab} e_a^i e_b^j = \perp^{ij}$, where $\perp^{ij} = \delta^{ij} - n^i n^j$ is the projector onto the sphere tangent bundle. We also use the notation $e_{\langle a}^i e_{b\rangle}^j = e_{(a}^i e_{b)}^j - \frac{1}{2} \gamma_{ab} \perp^{ij}$ for the trace-free product of basis vectors, where round brackets denote symmetrization $e_{(a}^i e_{b)}^j = \frac{1}{2}(e_a^i e_b^j + e_b^i e_a^j)$. On the other hand, square brackets denote anti-symmetrization, e.g., $e_{[a}^i e_{b]}^j = \frac{1}{2}(e_a^i e_b^j - e_b^i e_a^j)$. The transverse-trace-free (TT) projection operator is defined as $\perp_{TT}^{ijkl} = \perp^{k(i} \perp^{j)l} - \frac{1}{2} \perp^{ij} \perp^{kl}$. Introducing the covariant derivative $D_a$ compatible with the sphere metric, $D_a \gamma_{bc} = 0$, we can write $D_a e_b^i = D_b e_a^i = D_a D_b n^i = -\gamma_{ab} n^i$. It is also convenient to denote $\Delta = D_c D^c$. An arbitrary
rank 2 STF tensor over the sphere may be decomposed as

$$T_{ab} = -2D_{(a}D_{b)}T^{+} + 2\epsilon_{c(a}D_{b)}D^{c}T^{-}, \qquad (1.1)$$

where $\epsilon_{ab} \equiv e_a^i e_b^j n_k \epsilon_{kij}$ is the unit sphere volume form and the modal decomposition of $T^{\pm} = \sum_{\ell \ge 2} T_L^{\pm} n_L$, with the $T_L^{\pm}$ being STF over the multi-index $L$, contains only $\ell \ge 2$ harmonics without loss of generality. Note the properties $D_{(a}e_{b)}^i = 0$, $D^a(e_{(a}^i e_{b)}^j) = -2n^{(i}e_{b)}^j$, $\Delta T^{\pm} = -\ell(\ell+1)n_L T_L^{\pm}$, and $(\Delta + 2)(e_{(a}^i e_{b)}^j T_{ij}^+) = 0$. We define the unit-normalized integral over the sphere as $\oint_S 1 = 1$.

Given a general manifold, harmonic/de Donder coordinates are specified by using a tilde: $\tilde{x}^{\mu} = (\tilde{t}, \tilde{\mathbf{x}})$ or $(\tilde{t}, \tilde{r}, \tilde{\theta}^a)$. The metric components are $\tilde{g}_{\mu\nu}(\tilde{x})$. Asymptotically flat spacetimes admit, as a background structure, the Minkowskian outgoing vector $\tilde{k}^{\mu} = (1, \tilde{n}^i)$, the basis on the sphere $\tilde{e}_a^i = \partial\tilde{n}^i/\partial\tilde{\theta}^a$, etc. We define the retarded time $\tilde{u}$ in harmonic coordinates as $\tilde{u} = \tilde{t} - \tilde{r}$, so that $\tilde{k}_{\mu} = -\tilde{\partial}_{\mu}\tilde{u}$.

NU coordinates are denoted as $x^{\mu} = (u, r, \theta^a)$ with $\theta^a = (\theta, \varphi)$. The metric components in those coordinates are $g_{\mu\nu}(x)$, with all other notations, such as the natural basis on the sphere $e_a^i$ and the metric $\gamma_{ab}$, as previously.

# 2 Quadratic memory from the MPM formalism

## 2.1 Original formulation: modified harmonic coordinates

Following the MPM formalism, the metric outside an isolated matter system is written as the formal post-Minkowskian (PM) expansion $h^{\mu\nu} = \sqrt{|\tilde{g}|} \tilde{g}^{\mu\nu} - \eta^{\mu\nu} = G h_1^{\mu\nu} + G^2 h_2^{\mu\nu} + \mathcal{O}(G^3)$, with each PM coefficient generated in the form of a multipole expansion starting from the linearized vacuum metric in harmonic coordinates [70–73] (setting $c = 1$)

$$h_1^{00} = -4 \sum_{\ell=0}^{+\infty} \frac{(-)^{\ell}}{\ell!} \tilde{\partial}_L \left( \frac{M_L(\tilde{u})}{\tilde{r}} \right), \qquad (2.1a)$$

$$h_1^{0i} = 4 \sum_{\ell=1}^{+\infty} \frac{(-)^{\ell}}{\ell!} \left[ \tilde{\partial}_{L-1} \left( \frac{M_{iL-1}(\tilde{u})}{\tilde{r}} \right) + \frac{\ell}{\ell+1} \tilde{\partial}_{pL-1} \left( \frac{\epsilon_{ipq} S_{qL-1}(\tilde{u})}{\tilde{r}} \right) \right], \qquad (2.1b)$$

$$h_1^{ij} = -4 \sum_{\ell=2}^{+\infty} \frac{(-)^{\ell}}{\ell!} \left[ \tilde{\partial}_{L-2} \left( \frac{M_{ijL-2}(\tilde{u})}{\tilde{r}} \right) + \frac{2\ell}{\ell+1} \tilde{\partial}_{pL-2} \left( \frac{\epsilon_{pq(i} S_{j)qL-2}(\tilde{u})}{\tilde{r}} \right) \right]. \qquad (2.1c)$$

Here $M_L$ and $S_L$ are the symmetric-trace-free (STF) canonical mass and current multipole moments, which depend on the harmonic coordinate retarded time $\tilde{u} = \tilde{t} - \tilde{r}$. For the following, it is convenient to write the dominant $1/\tilde{r}$ piece of $h_1^{\mu\nu}$ (when $\tilde{r} \to +\infty$ with
fixed $\tilde{u}$) as<fnref n="2" />

$$h_1^{\mu\nu} = \frac{1}{\tilde{r}} \left( \begin{array}{c} -4(M + \tilde{n}_i P^i) + z_1^{00}(\tilde{u}, \tilde{n}) \\ -4P^i + z_1^{0i}(\tilde{u}, \tilde{n}) \\ z_1^{ij}(\tilde{u}, \tilde{n}) \end{array} \right) + \mathcal{O}\left(\frac{1}{\tilde{r}^2}\right), \qquad (2.2)$$

where $M$ is the constant Arnowitt-Deser-Misner (ADM) energy, $S_i$ is the constant ADM angular momentum and $P^i \equiv M_i^{(1)}$ the constant ADM linear momentum of the system, $M_i$ being the mass dipole, identical for the conservative dynamics to the vector $G_i$ that defines the center-of-mass position (after division by $M$) and reduces to zero in the center-of-mass frame. We have introduced the quantities

$$z_1^{00} = -4 \sum_{\ell=2}^{+\infty} \frac{1}{\ell!} \tilde{n}_L M_L^{(\ell)}(\tilde{u}), \qquad (2.3a)$$

$$z_1^{0i} = -4 \sum_{\ell=2}^{+\infty} \frac{1}{\ell!} \left( \tilde{n}_{L-1} M_{iL-1}^{(\ell)}(\tilde{u}) + \frac{\ell}{\ell+1} \tilde{n}_{pL-1} \epsilon_{ipq} S_{qL-1}^{(\ell)}(\tilde{u}) \right), \qquad (2.3b)$$

$$z_1^{ij} = -4 \sum_{\ell=2}^{+\infty} \frac{1}{\ell!} \left( \tilde{n}_{L-2} M_{ijL-2}^{(\ell)}(\tilde{u}) + \frac{2\ell}{\ell+1} \tilde{n}_{pL-2} \epsilon_{pq(i} S_{j)qL-2}^{(\ell)}(\tilde{u}) \right). \qquad (2.3c)$$

Note that $\tilde{k}_\nu z_1^{\mu\nu} = 0$ as the consequence of the harmonic gauge condition $\tilde{\partial}_\nu h_1^{\mu\nu} = 0$, where we have introduced the Minkowski outgoing null vector $\tilde{k}^\mu = (1, \tilde{n}^i)$. At quadratic order, $h_2^{\mu\nu}$ obeys the flat wave equation $\tilde{\Box}h_2^{\mu\nu} = N_2^{\mu\nu}$ (together with $\tilde{\partial}_\nu h_2^{\mu\nu} = 0$), where the source term $N_2^{\mu\nu}$ is quadratic in $h_1$ as well as its first and second space-time derivatives.

Hereditary terms with respect to the canonical moments $M_L, S_L$ are defined as retarded time integrals involving those moments; for instance tail terms are hereditary terms with typically some logarithmic kernel, and memory terms are given by time anti-derivatives of product of moments. It is known [5] that, at quadratic order, hereditary terms with respect to the canonical moments $M_L, S_L$ are generated only from the $1/\tilde{r}^2$ piece in $N_2^{\mu\nu}$,

$$N_2^{\mu\nu} = \frac{1}{\tilde{r}^2} Q_2^{\mu\nu}(\tilde{u}, \tilde{n}) + \mathcal{O}\left(\frac{1}{\tilde{r}^3}\right), \qquad (2.4)$$

and that this piece takes the form

$$Q_2^{\mu\nu}(\tilde{u}, \tilde{n}) = 4(M + \tilde{n}_i P_i) \frac{d^2 z^{\mu\nu}}{d\tilde{u}^2} + \tilde{k}^\mu \tilde{k}^\nu \Pi(\tilde{u}, \tilde{n}). \qquad (2.5)$$

The first term will generate the tails while the second term is responsible for the displacement memory. The latter takes the form of the stress-energy tensor of gravitons propagating along the null direction $\tilde{k}^\mu = (1, \tilde{n}^i)$. The quantity $\Pi$ is directly given by a quadratic product of the multipole expansion (2.3):

$$\Pi = \frac{1}{2} \frac{dz^{\mu\nu}}{d\tilde{u}} \frac{dz_{\mu\nu}}{d\tilde{u}} - \frac{1}{4} \frac{dz_{\mu}^{\mu}}{d\tilde{u}} \frac{dz_{\nu}^{\nu}}{d\tilde{u}}. \qquad (2.6)$$

<footnote n="2">²This tensorial quantity is presented in the form $h^{\mu\nu} = \begin{pmatrix} h^{00} \\ h^{0i} \\ h^{ij} \end{pmatrix}$.</footnote>
The GW energy flux emitted in the solid angle $d\tilde{\Omega}$ around the direction $\tilde{n}$ is

$$ \frac{dE^{\text{GW}}}{d\tilde{u}d\tilde{\Omega}} = \frac{G}{16\pi} \Pi(\tilde{u}, \tilde{n}) + \mathcal{O}(G^2). \qquad (2.7) $$

The systematic study of hereditary terms (either tails or memory) at quadratic order was done in [5] with such formalism. In addition to the tail and memory terms, there are other hereditary terms in the harmonic gauge at this order. They are associated with GW losses of energy, linear and angular momenta, as well as center-of-mass position. The leading $1/\tilde{r}$ components of $h_2^{\mu\nu}$ for multipole-multipole interactions in a radiative metric are given in eq. (2.39) of [5]. Moreover, the hereditary terms in the subleading piece $1/\tilde{r}^2$ are easily computed from (2.29-30) of that paper. We can write

$$ h_2^{\mu\nu}|_{\text{hered}} = h_2^{\mu\nu}|_{\text{tail}} + h_2^{\mu\nu}|_{\text{mem}}. \qquad (2.8) $$

The tail term comes directly from the first term in eq. (2.5); at the dominant order $1/\tilde{r}$, it involves a logarithmic kernel,

$$ h_2^{\mu\nu}|_{\text{tail}} = \frac{2(M + \tilde{n}_i P_i)}{\tilde{r}} \int_{-\infty}^{\tilde{u}} dv \ln \left( \frac{\tilde{u} - v}{2b_0} \right) \frac{d^2 z^{\mu\nu}}{d\tilde{u}^2} (v, \tilde{n}) + \mathcal{O}\left(\frac{1}{\tilde{r}^2}\right). \qquad (2.9) $$

It will not be considered here since it has already been investigated in the previous paper [61]. The memory term actually contains both contributions mentioned above, from the non-linear memory strictly speaking and from the GW secular losses. For simplicity, we will just refer to these two types of terms as “memory”. The memory is then given at quadratic order by the “exact expression”, i.e., valid at any order in $1/\tilde{r}$:

$$ h_2^{\mu\nu}|_{\text{mem}} = \frac{1}{\tilde{r}} \int_{-\infty}^{\tilde{u}} dv K^{\mu\nu}(v, \tilde{n}) + \frac{1}{\tilde{r}^2} \begin{pmatrix} \tilde{n}_i \left[ \frac{1}{3} \frac{\Pi_i}{\tilde{r}^2} + 4 \frac{\Lambda_i}{\tilde{r}^2} \right] \\ -2\epsilon_{ipq} \tilde{n}_p \frac{\Sigma_q}{\tilde{r}^2} \\ 0 \end{pmatrix}. \qquad (2.10) $$

The radiative coordinates used in this formula are defined in [12] and in eqs. (2.17)-(2.19) of [5]; they are not Newman-Unti coordinates but, since they are designed to remove all the $\ln \tilde{r}$ components of the radiative metric to quadratic order, they belong, like the NU coordinates, to the large class of radiative coordinate systems [47–49]. Although these coordinates are not harmonic either (we may call them modified harmonic coordinates), we will still use tildes to denote them, $\tilde{x}^{\mu} = (\tilde{u}, \tilde{r}, \tilde{n})$.

The leading term in $1/\tilde{r}$ is just a time anti-derivative (sometimes called “semi-hereditary”). It is a combination of the memory and also the secular mass loss $\propto \Pi_0$ in the $00$ component of $K^{\mu\nu}$. To define it, we decompose the quantity $\Pi$ into STF multipolar pieces as

$$ \Pi(\tilde{u}, \tilde{n}) = \sum_{\ell=0}^{+\infty} \tilde{n}_\ell \Pi_L(\tilde{u}), \qquad (2.11a) $$

$$ \text{where } \Pi_L(\tilde{u}) = \frac{(2\ell + 1)!!}{\ell!} \int \frac{d\tilde{\Omega}}{4\pi} \tilde{n}_{\langle L \rangle} \Pi(\tilde{u}, \tilde{n}). \qquad (2.11b) $$
Then, the components of $K^{\mu\nu}$ are explicitly given by<fnref n="3" />

$$K^{00}(\tilde{u}, \tilde{n}) = \frac{1}{2}\Pi_0 + \frac{1}{3}\tilde{n}_i \Pi_i, \qquad (2.12a)$$

$$K^{0i}(\tilde{u}, \tilde{n}) = \frac{1}{3}\Pi_i + \frac{1}{2}\sum_{\ell \ge 0} \frac{1}{\ell+1}\tilde{n}_{iL}\Pi_L - \frac{1}{2}\sum_{\ell \ge 1} \frac{1}{\ell+1}\tilde{n}_{L-1}\Pi_{iL-1}, \qquad (2.12b)$$

$$K^{ij}(\tilde{u}, \tilde{n}) = \sum_{\ell \ge 0} \frac{1}{\ell + 2} \tilde{n}_{ijL} \Pi_L - 2 \sum_{\ell \ge 2} \frac{1}{(\ell + 1)(\ell + 2)} \tilde{n}_{L-2} \Pi_{ijL-2} \\ - \sum_{\ell \ge 1} \frac{1}{(\ell + 1)(\ell + 2)} \left[ \delta_{ij} \tilde{n}_L \Pi_L + (\ell - 2) \tilde{n}_{L-1(i} \Pi_{j)L-1} \right]. \qquad (2.12c)$$

Notice the useful properties, posing $K \equiv \eta_{\mu\nu}K^{\mu\nu}$:

$$\tilde{k}_{\nu} K^{\mu\nu} = 0, \qquad K = -\frac{1}{3} \tilde{n}_i \Pi_i. \qquad (2.13)$$

To order $1/\tilde{r}$, the memory contribution to the transverse-traceless (TT) asymptotic wave-form reads

$$H_{ij}^{\text{TT}}|_{\text{mem}} = -\perp_{ijkl}^{\text{TT}} \int_{-\infty}^{\tilde{u}} dv K^{kl}(v, \tilde{n}) = \perp_{ijkl}^{\text{TT}} \sum_{\ell \ge 2} \frac{2\tilde{n}_{L-2}}{(\ell+1)(\ell+2)} \int_{-\infty}^{\tilde{u}} dv \Pi_{klL-2}(v), \quad (2.14)$$

since all terms in $K^{kl}$ proportional to either $\tilde{n}^k$, $\tilde{n}^l$ or $\delta^{kl}$ are killed by the TT projection. The global minus sign arises from the relation $g_{ij}^{\text{TT}} = -h_{ij}^{\text{TT}}$, where $h_{ij}$ is the perturbation of the gothic metric (following the convention of [61]).

As for the subleading $1/\tilde{r}^2$ piece in eq. (2.10), it is associated with GW secular losses; it contains both simple and double anti-derivatives indicated by the superscripts $(-n)$. We assume that the metric was stationary before some remote date in the past (for any $\tilde{u} \le -T$). In this situation, all anti-derivatives are well defined.

From eq. (2.7), the angular average $\Pi_0 = \int \frac{d\tilde{\Omega}}{4\pi} \Pi$ is proportional to the total energy flux carried by gravitational waves:

$$\frac{dE^{\text{GW}}}{d\tilde{u}} = \frac{1}{4} G \Pi_0 + \mathcal{O}(G^2). \qquad (2.15)$$

Similarly, $\Pi_i = 3 \int \frac{d\tilde{\Omega}}{4\pi} \tilde{n}_i \Pi$ is proportional to the total linear momentum flux. We also introduced the notation $\Sigma_i$ for the angular momentum flux and $\Lambda_i$ for the center-of-mass flux:

$$\frac{dP_i^{\text{GW}}}{d\tilde{u}} = \frac{G \Pi_i}{12} + \mathcal{O}(G^2), \quad \frac{dS_i^{\text{GW}}}{d\tilde{u}} = G \Sigma_i + \mathcal{O}(G^2), \quad \frac{dG_i^{\text{GW}}}{d\tilde{u}} - P_i^{\text{GW}} = G \Lambda_i + \mathcal{O}(G^2). \qquad (2.16)$$

We can write flux-balance equations equating the fluxes to losses in the source; see eqs. (2.29)-(2.30) of [5]. The multipole decompositions of the fluxes are well known [72, 74], except for the one associated with the position of the center-of-mass, obtained only recently from the asymptotic properties of the radiation field [75, 76], using the traditional PN as

<footnote n="3">³We corrected a sign typo in the last term of eq. (2.40b) of [5].</footnote>
well as the Bondi approaches [28, 60, 65, 66]. The complete multipole expansions of the fluxes to quadratic order ($\propto G^2$) are given by<fnref n="4" />

$$ \Pi_0 = 4 \sum_{\ell=2}^{+\infty} \left\{ \frac{(\ell+1)(\ell+2)}{(\ell-1)\ell\ell!(2\ell+1)!!} M_L^{(\ell+1)} M_L^{(\ell+1)} + \frac{4\ell(\ell+2)}{(\ell-1)(\ell+1)!(2\ell+1)!!} S_L^{(\ell+1)} S_L^{(\ell+1)} \right\}, \quad (2.17a) $$

$$ \Pi_i = 12 \sum_{\ell=2}^{+\infty} \left\{ \frac{2(\ell+2)(\ell+3)}{\ell(\ell+1)!(2\ell+3)!!} M_{iL}^{(\ell+2)} M_L^{(\ell+1)} + \frac{8(\ell+3)}{(\ell+1)!(2\ell+3)!!} S_{iL}^{(\ell+2)} S_L^{(\ell+1)} \right. \\ \left. + \frac{8(\ell+2)}{(\ell-1)(\ell+1)!(2\ell+1)!!} \epsilon_{ijk} M_{jL-1}^{(\ell+1)} S_{kL-1}^{(\ell+1)} \right\}, \quad (2.17b) $$

$$ \Sigma_i = \epsilon_{ijk} \sum_{\ell=2}^{+\infty} \left\{ \frac{(\ell+1)(\ell+2)}{(\ell-1)\ell!(2\ell+1)!!} M_{jL-1}^{(\ell)} M_{kL-1}^{(\ell+1)} \right. \\ \left. + \frac{4\ell^2(\ell+2)}{(\ell-1)(\ell+1)!(2\ell+1)!!} S_{jL-1}^{(\ell)} S_{kL-1}^{(\ell+1)} \right\}, \quad (2.17c) $$

$$ \Lambda_i = \sum_{\ell=2}^{+\infty} \left\{ \frac{(\ell+2)(\ell+3)}{\ell\ell!(2\ell+3)!!} \left( \begin{array}{cc} (\ell+1) & (\ell+1) \\ M_{iL} & M_L \end{array} - \begin{array}{cc} (\ell) & (\ell+2) \\ M_L & M_{iL} \end{array} \right) \right. \\ \left. + \frac{4(\ell+3)}{\ell!(2\ell+3)!!} \left( \begin{array}{cc} (\ell+1) & (\ell+1) \\ S_{iL} & S_L \end{array} - \begin{array}{cc} (\ell) & (\ell+2) \\ S_L & S_{iL} \end{array} \right) \right\}. \quad (2.17d) $$

The center-of-mass flux formula $\Lambda_i$ was derived in [60] and matches the formula obtained by Nichols in terms of spherical harmonics [65] after erratum. It also reproduces the low harmonic results of [75, 76]. In [66], the following simpler alternative center-of-mass flux was derived:

$$ \tilde{\Lambda}_i = \sum_{\ell=2}^{+\infty} \left\{ \frac{2(\ell+2)(\ell+3)}{\ell\ell!(2\ell+3)!!} M_{iL}^{(\ell+1)} M_L^{(\ell+1)} + \frac{8(\ell+3)}{\ell!(2\ell+3)!!} S_{iL}^{(\ell+1)} S_L^{(\ell+1)} \right\}. \quad (2.18) $$

It is associated with a different center-of-mass vector $\tilde{G}_i$. Since $\tilde{\Lambda}_i - \Lambda_i = \partial_{\tilde{u}} \Delta G_i$ is a total $\tilde{u}$ derivative, which vanishes on average for periodic systems, the two formulations of the flux-balance law are physically equivalent. One can relate the corresponding center-of-mass vectors by $\tilde{G}_i = G_i + \Delta G_i$ and the flux-formulae differ due to this shift at fixed retarded time $\tilde{u}$. A disadvantage of the shifted $\tilde{G}_i$, as was shown in [60], is that it cannot be written by means of a covariant formula over the sphere in Bondi gauge. We will use the definition $G_i$ and the corresponding flux $\Lambda_i$ henceforth.

Notice that one only needs the leading term (2.3) in the $1/\tilde{r}$ expansion of the linearized metric in order to compute the energy and linear momentum fluxes, whereas in the case of the angular momentum and center-of-mass fluxes, one also needs to include the next-to-leading correction $1/\tilde{r}^2$ (see [66] for details). On the other hand, the Poincaré fluxes become exact expressions when rewritten in terms of radiative multipole moments $U_L = M_L^{(\ell)} + \mathcal{O}(G)$, $V_L = S_L^{(\ell)} + \mathcal{O}(G)$ [60].

<footnote n="4">⁴Fluxes are defined with the opposite sign as the corresponding ADM quantities. Accordingly, the fluxes of Poincaré charges $\mathcal{E}$, $\mathcal{P}_i$, $\mathcal{J}_i$ and $\mathcal{G}_i$ defined from canonical methods [60] are given in the notations used here as $\dot{\mathcal{E}} = -\frac{d\mathcal{E}^{GW}}{d\tilde{u}}$, $\dot{\mathcal{P}}_i = -\frac{d\mathcal{P}_i^{GW}}{d\tilde{u}}$, $\dot{\mathcal{J}}_i = -\frac{d\mathcal{S}_i^{GW}}{d\tilde{u}}$ and $\dot{\mathcal{G}}_i = \mathcal{P}_i - \frac{d\mathcal{G}_i^{GW}}{d\tilde{u}} + \mathcal{P}_i^{GW}$. These definitions correspond to $\alpha = 1 = \beta$ in the parametrization given in eq. (3.18) of [60]. The Kerr black hole has $\mathcal{E} = M$ and $\mathcal{J}_z = Ma$ with the standard orientation $\epsilon_{tr\theta\phi} = 1$.</footnote>
## 2.2 Quadratic memory in Newman-Unti coordinates

Given a metric perturbation $h^{\mu\nu} = G h_1^{\mu\nu} + G^2 h_2^{\mu\nu} + \mathcal{O}(G^3)$ in an arbitrary (not necessarily harmonic) coordinate system $\{\tilde{u}, \tilde{r}, \tilde{\theta}^a\}$, we constructed perturbatively in [61] the gauge transformation $\{\tilde{u}, \tilde{r}, \tilde{\theta}^a\} \rightarrow \{u, r, \theta^a\}$ required to reach the Newman-Unti (NU) gauge. Denoting the coordinate transformation as

$$u = \tilde{u} + G U_1(\tilde{u}, \tilde{r}, \tilde{\theta}^a) + G^2 U_2(\tilde{u}, \tilde{r}, \tilde{\theta}^a) + \mathcal{O}(G^3), \quad (2.19a)$$

$$r = \tilde{r} + G R_1(\tilde{u}, \tilde{r}, \tilde{\theta}^a) + G^2 R_2(\tilde{u}, \tilde{r}, \tilde{\theta}^a) + \mathcal{O}(G^3), \quad (2.19b)$$

$$\theta^a = \tilde{\theta}^a + G \Theta_1^a(\tilde{u}, \tilde{r}, \tilde{\theta}^a) + G^2 \Theta_2^a(\tilde{u}, \tilde{r}, \tilde{\theta}^a) + \mathcal{O}(G^3), \quad (2.19c)$$

For the transformed metric to satisfy NU gauge $g_{rr} = 0 = g_{ra}$ and $g_{ur} = -1$, we have to solve successively the linear order equations (where $\tilde{k}^\mu$ is the Minkowski null outgoing vector defined in the original coordinates)

$$\tilde{k}^{\mu} \tilde{\partial}_{\mu} U_{1} = \frac{1}{2} \tilde{k}_{\mu} \tilde{k}_{\nu} h_{1}^{\mu\nu}, \quad (2.20a)$$

$$\tilde{k}^{\mu} \tilde{\partial}_{\mu} R_{1} = \frac{1}{2} h_{1} + \tilde{n}_{i} \left( \tilde{\partial}_{i} U_{1} - \tilde{k}_{\mu} h_{1}^{\mu i} \right), \quad (2.20b)$$

$$\tilde{k}^{\mu} \tilde{\partial}_{\mu} \Theta_{1}^{a} = \frac{\tilde{\epsilon}_{i}^{a}}{\tilde{r}} \left( \tilde{\partial}_{i} U_{1} - \tilde{k}_{\mu} h_{1}^{\mu i} \right), \quad (2.20c)$$

and next, the quadratic order equations

$$\tilde{k}^{\mu} \tilde{\partial}_{\mu} U_{2} = \frac{1}{2} \tilde{k}_{\mu} \tilde{k}_{\nu} h_{2}^{\mu\nu} + \left( \frac{1}{2} \tilde{\partial}^{\mu} U_{1} - \tilde{k}_{\nu} h_{1}^{\mu\nu} \right) \tilde{\partial}_{\mu} U_{1}, \quad (2.21a)$$

$$\tilde{k}^{\mu} \tilde{\partial}_{\mu} R_{2} = \frac{1}{8} h_{1}^{2} - \frac{1}{4} h_{1}^{\mu\nu} h_{1\mu\nu} + \frac{1}{2} h_{2} + \tilde{n}_{i} \left( \tilde{\partial}_{i} U_{2} - \tilde{k}_{\mu} h_{2}^{\mu i} + (\tilde{\partial}_{\mu} U_{1}) h_{1}^{\mu i} \right) \\ + \left( \tilde{\partial}^{\mu} U_{1} - \tilde{k}_{\nu} h_{1}^{\mu\nu} \right) \tilde{\partial}_{\mu} R_{1}, \quad (2.21b)$$

$$\tilde{k}^{\mu} \tilde{\partial}_{\mu} \Theta_{2}^{a} = \frac{\tilde{\epsilon}_{i}^{a}}{\tilde{r}} \left( \tilde{\partial}_{i} U_{2} - \tilde{k}_{\mu} h_{2}^{\mu i} + (\tilde{\partial}_{\mu} U_{1}) h_{1}^{\mu i} \right) + \left( \tilde{\partial}^{\mu} U_{1} - \tilde{k}_{\nu} h_{1}^{\mu\nu} \right) \tilde{\partial}_{\mu} \Theta_{1}^{a}. \quad (2.21c)$$

In this section, we compute the memory-type and mass-loss-type hereditary terms in the NU metric, following the previous algorithm. The tail-type hereditary terms (in the quadrupole case) have been investigated in [61] and we will only report the result. Note that the first order perturbation $h_1^{\mu\nu}$ does not contain hereditary integrals (it is instantaneous in terms of the canonical moments $M_L$ and $S_L$). Only the second order metric $h_2^{\mu\nu}$ involves the mass-loss and memory terms. Moreover, the differential operator $\tilde{k}^\mu \tilde{\partial}_\mu = \tilde{\partial}_r|_{\tilde{u}}$ is simply solved by integration over $\tilde{r}$ at constant $\tilde{u}$. Therefore, no hereditary integrals can be generated from instantaneous (non-hereditary) terms on the right-hand side of the equations (2.21) and, in order to control them, it is sufficient to solve the linear equations

$$\tilde{k}^{\mu} \tilde{\partial}_{\mu} U_{2} = \frac{1}{2} \tilde{k}_{\mu} \tilde{k}_{\nu} h_{2}^{\mu\nu}, \quad (2.22a)$$

$$\tilde{k}^{\mu} \tilde{\partial}_{\mu} R_{2} = \frac{1}{2} h_{2} + \tilde{n}_{i} \left( \tilde{\partial}_{i} U_{2} - \tilde{k}_{\mu} h_{2}^{\mu i} \right), \quad (2.22b)$$

$$\tilde{k}^{\mu} \tilde{\partial}_{\mu} \Theta_{2}^{a} = \frac{\tilde{\epsilon}_{i}^{a}}{\tilde{r}} \left( \tilde{\partial}_{i} U_{2} - \tilde{k}_{\mu} h_{2}^{\mu i} \right). \quad (2.22c)$$
We now implement our algorithm by solving the system (2.22), focusing on the memory and mass-loss terms. Reminding the property $\tilde{k}_{\nu}K^{\mu\nu} = 0$, we see from eq. (2.10) that the first equation (2.22a) reduces to $\tilde{k}^{\mu}\tilde{\partial}_{\mu}U_2 = \mathcal{O}(\tilde{r}^{-2})$ at leading order, which is directly solved to $U_2 = U_2^0(\tilde{u}, \tilde{n}) + \mathcal{O}(\tilde{r}^{-1})$. After imposing asymptotic flatness and setting the BMS transformation at second order to zero, we get $U_2^0(\tilde{u}, \tilde{n}) = 0$. Therefore, the only contributions to $U_2$ come from the $1/\tilde{r}^2$ term in $h_2^{\mu\nu}$ as given by eq. (2.10). This term is easily integrated and we obtain

$$U_2|_{\text{mem}} = \frac{\tilde{n}_i}{\tilde{r}} \left( -\frac{1}{6} \Pi_i^{(-2)} + 2 \Lambda_i^{(-1)} \right). \qquad (2.23)$$

We recall that $\Pi_i$ is the flux of linear momentum while $\Lambda_i$ is the flux of center-of-mass [see eqs. (2.16)-(2.17)], all quantities being evaluated at $\tilde{u} \equiv \tilde{t} - \tilde{r}$ and $\tilde{x}^a$. Continuing the algorithm, we find that $R_2$ does not contain memory/mass-loss terms whereas $\Theta_2^a$ receives contributions from $\Pi_i$, $\Lambda_i$, as well as the angular momentum flux $\Sigma_i$:

$$R_2|_{\text{mem}} = 0, \qquad (2.24a)$$

$$\Theta_2^a|_{\text{mem}} = \frac{e_i^a}{\tilde{r}^2} \left( \frac{1}{12} \Pi_i^{(-2)} + \Lambda_i^{(-1)} + \epsilon_{ipq} \tilde{n}_p \Sigma_q^{(-1)} \right). \qquad (2.24b)$$

We refer to the previous paper [61] for the computation of the tail contributions due to mass quadrupole in the quantities (2.23)-(2.24).

We proceed with the next steps as described in [61]: we successively compute the contravariant components of the metric, $g^{rr}$, $g^{ra}$ and $g^{ab}$, and deduce its covariant components in the NU gauge as $g_{ua} = g^{rb}g_{ab}$, $g_{uu} = -g^{rr} + g^{ra}g_{ua}$ and $g_{ab} = (g^{ab})^{-1}$. In the end, we re-express the metric in terms of the NU coordinates $\{u, r, \theta^a\}$ using the inverse of eqs. (2.19). In particular, this entails re-expanding the 2-sphere metric as

$$\tilde{r}^2 \tilde{\gamma}_{ab} = r^2 \left[ \gamma_{ab} - 2G^2 \left( r^{-1} R_2 \gamma_{ab} + \Theta_2^c \Gamma_{c(a}^e \gamma_{b)e} \right) \right] + \mathcal{O}(G^3), \qquad (2.25)$$

where $\Gamma_{bc}^a$ is the Christoffel symbol associated with the covariant derivative on the sphere. This brings a memory correction coming from $\Theta_2^a$ by virtue of (2.24). We find that the Christoffel symbols cancel out so that our final metric is covariant with respect to diffeomorphisms acting on the sphere. We finally get

$$g_{uu}|_{\text{mem}} = -1 - G^2 \left[ \frac{1}{r} \left( \frac{1}{2} \Pi_0^{(-1)} + \frac{1}{2} n_i \Pi_i^{(-1)} \right) + \frac{n_i}{r^2} \left( \frac{1}{6} \Pi_i^{(-2)} + 2 \Lambda_i^{(-1)} \right) \right] + \mathcal{O}(G^3), \qquad (2.26a)$$

$$g_{ua}|_{\text{mem}} = -G^2 e_a^i \left[ \frac{1}{2} \sum_{\ell \ge 2} \frac{1}{\ell+1} n_{L-1} \Pi_{iL-1}^{(-1)} + \frac{1}{r} \left( \frac{1}{6} \Pi_i^{(-2)} + 2 \Lambda_i^{(-1)} + 2 \epsilon_{ipq} n_p \Sigma_q^{(-1)} \right) \right] + \mathcal{O}(G^3), \qquad (2.26b)$$

$$g_{ab}|_{\text{mem}} = r^2 \gamma_{ab} + 2G^2 r e_{(a}^i e_{b)}^j \sum_{\ell \ge 2} \frac{1}{(\ell+1)(\ell+2)} n_{L-2} \Pi_{ijL-2}^{(-1)} + \mathcal{O}(G^3). \qquad (2.26c)$$

The above metric is entirely expressed in terms of the NU coordinates $\{u, r, \theta^a\}$. Now, the general asymptotically flat solution of interest to the vacuum Einstein's equations in
Newman-Unti coordinates, up to $O(r^{-3})$ corrections, reads [33, 50, 51, 57–59, 61]

$$g_{uu} = -1 + \frac{2(m + \frac{1}{8}C_{ab}N^{ab})}{r} - \frac{D_a N^a}{3r^2} + \mathcal{O}(r^{-3}), \qquad (2.27a)$$

$$g_{ua} = \frac{D^b C_{ab}}{2} + \frac{2N_a}{3r} + \mathcal{O}(r^{-2}), \qquad (2.27b)$$

$$g_{ab} = r^2 \gamma_{ab} + r C_{ab} + \frac{C^{cd}C_{cd}}{8} \gamma_{ab} + \mathcal{O}(r^{-1}), \qquad (2.27c)$$

together with $g_{rr} = g_{ra} = 0$ and $g_{ur} = -1$. Here, $m$ is the Bondi mass aspect, $C_{ab}$ is the shear,<fnref n="5" /> $N_{ab}$ is the news and $N_a$ is the angular momentum aspect.

We can thus immediately read off the memory terms in the Bondi mass aspect $m$, the angular momentum aspect $N_a$ and the Bondi shear $C_{ab}$, as

$$m|_{\text{mem}} = -G^2 \left( \frac{1}{4} \Pi_0^{(-1)} + \frac{1}{4} n^i \Pi_i^{(-1)} \right) + \mathcal{O}(G^3), \qquad (2.28a)$$

$$N_a|_{\text{mem}} = -G^2 e_a^i \left( \frac{1}{4} \Pi_i^{(-2)} + 3 \Lambda_i^{(-1)} + 3 \epsilon_{ipq} n_p \Sigma_q^{(-1)} \right) + \mathcal{O}(G^3), \qquad (2.28b)$$

$$C_{ab}|_{\text{mem}} = 2G^2 e_{\langle a}^i e_{b \rangle}^j \sum_{\ell \ge 2} \frac{1}{(\ell+1)(\ell+2)} n_{L-2} \Pi_{iL-2}^{(-1)} + \mathcal{O}(G^3). \qquad (2.28c)$$

There is no memory contribution in the combination $\dot{N}_a - D_a m$, consistently with one of the Einstein field equations. The shear (2.28c) is actually traceless and agrees with eq. (5.6) of [77] after expressing the linearized gothic metric as minus the standard linearized metric. In fact, comparing the metric (2.27) in Newman-Unti coordinates up to $O(r^{-3})$ with our results (2.26), we find that they are perfectly consistent. Namely, we obtain from eq. (2.28) the divergences

$$D_a N^a|_{\text{mem}} = G^2 n^i \left( \frac{1}{2} \Pi_i^{(-2)} + 6 \Lambda_i^{(-1)} \right) + \mathcal{O}(G^3), \qquad (2.29a)$$

$$D^b C_{ab}|_{\text{mem}} = -G^2 e_a^i \sum_{\ell \ge 2} \frac{1}{\ell+1} n_{L-1} \Pi_{iL-1}^{(-1)} + \mathcal{O}(G^3), \qquad (2.29b)$$

and we see that the first expression agrees with our result for the sub-dominant term $1/r^2$ in $g_{uu}$ while the second one agrees with twice the $r^0$ term in $g_{ua}$.

We can also express the NU metric corresponding to the memory terms of the quadrupole-quadrupole interaction $M_{ij} \times M_{ij}$ in terms of the canonical mass quadrupole $M_{ij}$. The corresponding quadrupole-quadrupole metric in harmonic coordinates was computed in [77]. From eqs. (4.5) of [77], for quadratic products of two mass quadrupoles, the only non-zero multipolar coefficients are

$$\Pi_0 = \frac{4}{5} M_{ij}^{(3)} M_{ij}^{(3)}, \quad \Pi_{ij} = -\frac{24}{7} M_{k\langle i}^{(3)} M_{j\rangle k}^{(3)}, \quad \Pi_{ijkl} = M_{\langle ij}^{(3)} M_{kl\rangle}^{(3)}. \qquad (2.30)$$

<footnote n="5">⁵The shear $C_{ab}$ may be set traceless by a suitable choice of origin for the radial coordinate [51].</footnote>
This yields

$$g_{uu}|_{\text{mem}} = -1 - \frac{2G^2}{5r} \int_{-\infty}^{u} dv \, M_{ij}^{(3)}(v) M_{ij}^{(3)}(v), \qquad (2.31a)$$

$$g_{ua}|_{\text{mem}} = G^2 e_a^i n^j \int_{-\infty}^{u} dv \left[ \frac{9}{8} M_{ki}^{(3)} M_{j\bar{k}}^{(3)} - \frac{4}{5r} \left( M_{ik}^{(2)} M_{jk}^{(3)} - M_{jk}^{(2)} M_{ik}^{(3)} \right) \right] (v), \qquad (2.31b)$$

$$g_{ab}|_{\text{mem}} = r^2 \left[ \gamma_{ab} + \frac{G^2}{r} e_{\langle a}^i e_{b\rangle}^j \int_{-\infty}^{u} dv \left( -\frac{4}{7} M_{ki}^{(3)} M_{jk}^{(3)} + \frac{1}{15} n_{kl} M_{\langle ij}^{(3)} M_{kl\rangle}^{(3)} \right) (v) \right]. \qquad (2.31c)$$

In terms of physical fluxes (2.17), the memory contributions to the Bondi mass and angular momentum aspects take the simple and explicit form

$$m|_{\text{mem}} = -G(E^{\text{GW}} + 3n_i P_i^{\text{GW}}) + \mathcal{O}(G^3), \qquad (2.32a)$$

$$N_a|_{\text{mem}} = -3G e_a^i (G_i^{\text{GW}} + \epsilon_{ipq} n_p S_q^{\text{GW}}) + \mathcal{O}(G^3). \qquad (2.32b)$$

Let us also derive, to end with, an interesting alternative expression for the shear, which is related to the asymptotic waveform by $C_{ab} = e_{\langle a}^i e_{b\rangle}^j H_{ij}^{\text{TT}}$ [see e.g., [61] and eq. (2.14) above]. Using eq. (2.7) and eq. (2.11), we obtain

$$C_{ab}|_{\text{mem}} = 8e_{\langle a}^i e_{b\rangle}^j \sum_{\ell \ge 2} \frac{(2\ell + 1)!!}{(\ell + 2)!} n_{L-2} \int d\Omega' \hat{n}'_{ijL-2} \frac{dE^{\text{GW}}}{d\Omega'}(u, \mathbf{n}'), \qquad (2.33)$$

where the infinite multipole series therein can be summed up in closed form as (see below for the proof)

$$e_{\langle a}^i e_{b\rangle}^j \sum_{\ell \ge 2} \frac{(2\ell + 1)!!}{(\ell + 2)!} n_{L-2} \hat{n}'_{ijL-2} = \frac{1}{2} e_{\langle a}^i e_{b\rangle}^j \frac{n'_i n'_j}{1 - \mathbf{n} \cdot \mathbf{n}'} \qquad (2.34)$$

This leads to the following elegant result, which constitutes the best interpretation of the non-linear memory effect as due to the re-radiation of GWs by gravitons [2, 4, 78]:<fnref n="6" />

$$C_{ab}|_{\text{mem}} = 4e_{\langle a}^i e_{b\rangle}^j \int d\Omega' \frac{n'_i n'_j}{1 - \mathbf{n} \cdot \mathbf{n}'} \frac{dE^{\text{GW}}}{d\Omega'}(u, \mathbf{n}'), \qquad (2.35)$$

where the factor $1/(1-n\cdot n')$ is reminiscent of the Liénard-Wiechert potentials in the case of massless gravitons. It exactly reproduces the result of Christodoulou [2] and Thorne [4] (see also the earlier results [79, 80]).

**Proof of the formula (2.34).** We consider the following TT projection with respect to the unit vector $\mathbf{n} = (n_i)$:

$$A_{ij}^{\text{TT}} \equiv \left[ \frac{n'_i n'_j}{1 - \mathbf{n} \cdot \mathbf{n}'} \right]^{\text{TT}} \equiv \perp_{ijkl} \frac{n'_k n'_l}{1 - \mathbf{n} \cdot \mathbf{n}'}, \qquad (2.36)$$

<footnote n="6">⁶The derivation of the result (2.35) in [2, 4, 78] is very different from the one adopted here, which is rather based on [3, 5].</footnote>
where we recall that the TT projection reads $\perp_{ijkl}^{\text{TT}} = \frac{1}{2}(\perp_{ik}\perp_{jl} + \perp_{jk}\perp_{il} - \perp_{ij}\perp_{kl})$, with the usual perpendicular operator $\perp_{ij} = \delta_{ij} - n_i n_j$. We expand the denominator in eq. (2.36) as a power series in $n \cdot n'$ (which is convergent as soon as $n \cdot n' < 1$), hence

$$A_{ij}^{\text{TT}} = \left[ \sum_{\ell=2}^{+\infty} n_{L-2} n'_{ijL-2} \right]^{\text{TT}}. \qquad (2.37)$$

The multi-index $L-2$ contains $\ell-2$ indices, namely $a_1 \cdots a_{\ell-2}$. Using eq. (A.21a) in [73], we transform the ordinary product of unit vectors $n'_{ijL-2}$ into a sum of STF products

$$A_{ij}^{\text{TT}} = \left[ \sum_{\ell=2}^{+\infty} n_{L-2} \sum_{k=0}^{[\frac{\ell}{2}]} \alpha_k^{\ell} \delta_{\{2K} \hat{n}'_{L-2K}} \right]^{\text{TT}}, \qquad (2.38a)$$

$$\text{with } \alpha_k^\ell = \frac{(2\ell - 4k + 1)!!}{(2\ell - 2k + 1)!!}. \qquad (2.38b)$$

The operation over indices {} is defined as the un-normalized sum over the smallest set of permutations of $i_1 \cdots i_\ell$ which makes the object symmetrical in $L = i_1 \cdots i_\ell$. The object $\delta_{\{2K} \hat{n}'_{L-2K}}$ contains $\ell$ indices $L = ijL-2$ (with, say, $i = a_\ell$ and $j = a_{\ell-1}$) of which $2k$ are displayed onto the product of $k$ Kronecker symbols denoted $\delta_{2K}$. As an example we have $\delta_{\{ab}n_c\} \equiv \delta_{ab}n_c + \delta_{bc}n_a + \delta_{ca}n_b$.

The point is that, among all the terms composing $\delta_{\{2K} \hat{n}'_{L-2K}}$, we can discard all those which contain either $\delta_{ij}$, $\delta_{ia_p}$ or $\delta_{ja_q}$, since such terms will be cancelled by the TT projection. This is obvious for $\delta_{ij}$; in the two other cases, this results from the fact that, after multiplication by $n_{L-2}$ in eq. (2.38a), $\delta_{ia_p}$ or $\delta_{ja_q}$ will yield some $n_i$ or $n_j$. So, we can rewrite the expression (2.38a) by excluding the indices $ij$ from the operation {} , which we indicate by underlining the two indices $ij$:

$$A_{ij}^{\text{TT}} = \left[ \sum_{\ell=2}^{+\infty} n_{L-2} \sum_{k=0}^{[\frac{\ell-2}{2}]} \alpha_k^{\ell} \delta_{\{2K} \hat{n}'_{ijL-2-2K}} \right]^{\text{TT}}. \qquad (2.39)$$

The next step is to notice from eq. (A.19) in [73] that the number of terms composing the object $\delta_{\{2K} \hat{n}'_{ijL-2-2K}}$ is $\frac{(\ell-2)!}{2^k k!(\ell-2-2k)!}$. When contracted with $n_{L-2}$, all these terms will merge into a single one for each values of $\ell$ and $k$. Thus, we have

$$A_{ij}^{\text{TT}} = \left[ \sum_{\ell=2}^{+\infty} \sum_{k=0}^{[\frac{\ell-2}{2}]} \alpha_k^{\ell} \frac{(\ell-2)!}{2^k k!(\ell-2-2k)!} n_{L-2-2K} \hat{n}'_{ijL-2-2K} \right]^{\text{TT}}. \qquad (2.40)$$

We change $\ell$ into $\ell + 2k$ and rewrite the previous expression as

$$A_{ij}^{\text{TT}} = \left[ \sum_{\ell=2}^{+\infty} \sum_{k=0}^{+\infty} \alpha_k^{\ell+2k} \frac{(\ell+2k-2)!}{2^k k!(\ell-2)!} n_{L-2} \hat{n}'_{ijL-2} \right]^{\text{TT}} = \left[ \sum_{\ell=2}^{+\infty} S_\ell n_{L-2} \hat{n}'_{ijL-2} \right]^{\text{TT}}, \qquad (2.41)$$

introducing the coefficient $S_\ell$ which is given as an infinite series over all integer values of $k$. However, we find, using the expression of the coefficients $\alpha_k^{\ell+2k}$ deduced from eq. (2.38b),
that this series can actually be re-summed in closed analytic form:

$$S_{\ell} = \frac{(2\ell + 1)!!}{(\ell - 2)!} \sum_{k=0}^{+\infty} \frac{(\ell + 2k - 2)!}{2^k k! (2\ell + 2k + 1)!!} = 2^{\ell+2} \frac{\Gamma(\ell + \frac{3}{2})}{\sqrt{\pi} \Gamma(\ell + 3)} = 2 \frac{(2\ell + 1)!!}{(\ell + 2)!}. \quad (2.42)$$

Hence we obtain the simple result

$$A_{ij}^{\text{TT}} = \left[ 2 \sum_{\ell=2}^{+\infty} \frac{(2\ell + 1)!!}{(\ell + 2)!} n_{L-2} \hat{n}'_{ijL-2} \right]^{\text{TT}}. \quad (2.43)$$

Recalling that $e_{\langle a}^i e_{b\rangle}^j \perp_{ijkl}^{TT} = e_{\langle a}^k e_{b\rangle}^l$, we have therefore proved the formula (2.34).

## 2.3 $M_{ij} \times M_{ij}$ and $M \times M_{ij}$ asymptotic data in NU gauge

We now obtain the complete NU metric corresponding to the quadrupole-quadrupole interaction $M_{ij} \times M_{ij}$ and tails $M \times M_{ij}$, i.e., including, besides the non-linear memory and mass-loss terms derived in the previous section and besides the tail terms obtained previously in [61], all the instantaneous (non-hereditary) terms.

In this complete calculation of the $M_{ij} \times M_{ij}$ interaction, we start from the metric in harmonic coordinates as given in [77], instead of the modified harmonic coordinates described in section 2.1, in which the metric concerning hereditary effects was already free of far-zone logarithms [see eq. (2.10)]. Applying our algorithm, we shall check that indeed all the hereditary terms besides the memory, and all associated far-zone logarithms, disappear in the end of the calculation. A more sophisticated treatment of hereditary terms, which is explained in the appendix A below, is however required. We shall of course recover in particular the memory terms computed in the previous section.<fnref n="7" />

Our final results for the NU metric (see section 3.2 for the link with the Bondi metric) read

$$g_{uu} = -1 + \frac{2(m + \frac{1}{8}C_{ab}N^{ab})}{r} - \frac{D_a N^a}{3r^2} + \sum_{n=3}^{6} \frac{1}{r^n} g_{uu}(n), \quad (2.44a)$$

$$g_{ua} = \frac{D^b C_{ab}}{2} + \frac{2N_a}{3r} + \sum_{n=2}^{5} \frac{1}{r^n} g_{ua}(n), \quad (2.44b)$$

$$g_{ab} = r^2 \left[ \left(1 + \frac{W}{r^2}\right) \gamma_{ab} + \frac{1}{r} e_{\langle a}^i e_{b\rangle}^j \left( H_{ij}^{\text{TT}} + \sum_{n=2}^{+\infty} \frac{1}{r^n} E_{ij}(n) \right) \right]. \quad (2.44c)$$

Here, $H_{ij}^{\text{TT}}$ is the radiation field, which is a free data relating to the Bondi shear as $C_{ab} = e_{\langle a}^i e_{b\rangle}^j H_{ij}^{\text{TT}}$. It can be decomposed into $\ell \ge 2$ mass/electric and current/magnetic radiative multipoles $U_L, V_L$ as

$$H_{ij}^{\text{TT}} = 4 \perp_{ijkl}^{\text{TT}} (\mathbf{n}) \sum_{\ell=2}^{+\infty} \frac{1}{\ell!} \left\{ \hat{n}_{L-2} U_{klL-2}(u) - \frac{2\ell}{\ell+1} \hat{n}_{pL-2} \epsilon_{pq}(k V_l) q_{L-2}(u) \right\}. \quad (2.45)$$

<footnote n="7">⁷Our practical calculation is done with the software Mathematica supplemented by the xAct package [81].</footnote>
The nonzero radiative multipole moments are given by (displaying only the terms that are linear in $M_{ij}$ or that correspond to the $M_{ij} \times M_{ij}$ and $M \times M_{ij}$ interactions)$^8$

$$ U_{ij} = GM_{ij}^{(2)} + 2G^2 M \int_0^{+\infty} dz \left[ \ln\left(\frac{z}{2b_0}\right) + \frac{11}{12} \right] M_{ij}^{(4)}(u-z) \\ + G^2 \left( -\frac{2}{7} \int_{-\infty}^{u} dv M_{k\langle i}^{(3)}(v) M_{j\rangle k}^{(3)}(v) - \frac{2}{7} M_{k\langle i}^{(2)} M_{j\rangle k}^{(3)} - \frac{5}{7} M_{k\langle i}^{(1)} M_{j\rangle k}^{(4)} + \frac{1}{7} M_{k\langle i} M_{j\rangle k}^{(5)} \right), \quad (2.46a) $$

$$ U_{ijkl} = G^2 \left( \frac{2}{5} \int_{-\infty}^{u} dv M_{\langle ij}^{(3)}(v) M_{kl\rangle}^{(3)}(v) - \frac{102}{5} M_{\langle ij}^{(2)} M_{kl\rangle}^{(3)} - \frac{63}{5} M_{\langle ij}^{(1)} M_{kl\rangle}^{(4)} - \frac{21}{5} M_{\langle ij} M_{kl\rangle}^{(5)} \right), \quad (2.46b) $$

$$ V_{ijk} = G^2 \epsilon_{pq\langle i} \left( \frac{1}{2} M_{j\bar{p}}^{(1)} M_{k\rangle q}^{(4)} - \frac{1}{10} M_{j\bar{p}} M_{k\rangle q}^{(5)} \right), \quad (2.46c) $$

where time derivatives of moments are denoted, e.g., $M_{ij}^{(q)} \equiv M_{ij}^{(q)}$, and underlined index must be regarded to be outside the STF projection. Here $b_0$ is a gauge constant related to the origin of time in radiative coordinates.

The remaining data in NU gauge are the mass aspect $m$ and the angular momentum aspects $N_i = e_i a N_a$, which read (adding also the leading linear terms)

$$ m = G(M + 3n^i P_i + 3\hat{n}^{ij} M_{ij}^{(2)}) \quad (2.47a) \\ + G^2 \left[ -\frac{1}{5} \int_{-\infty}^{u} dv M_{ij}^{(3)}(v) M_{ij}^{(3)}(v) + \hat{n}^{ij} \left( -\frac{6}{7} M_{ik}^{(2)} M_{jk}^{(3)} - \frac{15}{7} M_{ik}^{(1)} M_{jk}^{(4)} + \frac{3}{7} M_{ik} M_{jk}^{(5)} \right) \right. \\ \left. + \hat{n}^{ijkl} \left( -\frac{51}{4} M_{ij}^{(2)} M_{kl}^{(3)} - \frac{63}{8} M_{ij}^{(1)} M_{kl}^{(4)} - \frac{21}{8} M_{ij} M_{kl}^{(5)} \right) \right], $$

$$ N_i = 3G (M_i + 2n^j M_{ij}^{(1)})^T \quad (2.47b) \\ + G^2 \left[ \frac{6}{5} n^j \int_{-\infty}^{u} dv [M_{jk}^{(2)} M_{ik}^{(3)} - M_{ik}^{(2)} M_{jk}^{(3)}] (v) \right. \\ \left. - \frac{6}{35} n^j (14 M_{ik}^{(2)} M_{jk}^{(2)} + 22 M_{jk}^{(1)} M_{ik}^{(3)} + 22 M_{ik}^{(1)} M_{jk}^{(3)} + M_{jk} M_{ik}^{(4)} + M_{ik} M_{jk}^{(4)}) \right. \\ \left. - \frac{3}{2} \hat{n}^{ijkl} (9 M_{ij}^{(2)} M_{kl}^{(2)} + 10 M_{jk}^{(1)} M_{il}^{(3)} + 4 M_{ij}^{(1)} M_{kl}^{(3)} + 3 M_{jk} M_{il}^{(4)} + 4 M_{ij} M_{kl}^{(4)}) \right]^T. $$

The trace of the angular part of the metric (2.44c) is

$$ W = G^2 \left[ \frac{1}{5} M_{ij}^{(2)} M_{ij}^{(2)} - \frac{6}{7} \hat{n}^{ij} M_{ik}^{(2)} M_{jk}^{(2)} + \frac{1}{4} \hat{n}^{ijkl} M_{ij}^{(2)} M_{kl}^{(2)} \right. \\ \left. + \frac{1}{r^2} \left( \frac{3}{5} M_{ij} M_{ij}^{(2)} - \frac{18}{7} \hat{n}^{ij} M_{ik} M_{jk}^{(2)} + \frac{3}{4} \hat{n}^{ijkl} M_{ij} M_{kl}^{(2)} \right) \right. \\ \left. + \frac{1}{r^4} \left( \frac{7}{25} M_{ij} M_{ij} - \frac{6}{5} \hat{n}^{ij} M_{ik} M_{jk} + \frac{7}{20} \hat{n}^{ijkl} M_{ij} M_{kl} \right) \right]. \quad (2.48) $$

Finally, the non-zero subleading $1/r^n$ components of the metric (2.44), for the linear quadrupole metric and the quadratic quadrupole-quadrupole metric are given by (with

<footnote>$^8$For more general expressions involving current moments and higher multipoles of mass moments, see section 3.3 of [82] and references therein.</footnote>
(n)g<sub>ui</sub> ≡ e<sub>i</sub><sup>a</sup> (n)g<sub>ua</sub>):

$$
g_{(3)}^{uu} = 3GM_{ij}\hat{n}^{ij} + G^2 \left( -\frac{3}{5}M_{ij}^{(1)}M_{ij}^{(2)} - \frac{39}{7}\hat{n}^{ij}M_{ik}^{(1)}M_{jk}^{(2)} - \frac{281}{8}\hat{n}^{ijkl}M_{ij}^{(1)}M_{kl}^{(2)} \right. \\
\left. \qquad + \frac{1}{5}M_{ij}M_{ij}^{(3)} - \frac{3}{7}M_{ij}\hat{n}^{ik}M_{kj}^{(3)} - \frac{313}{8}M_{ij}\hat{n}^{ijkl}M_{kl}^{(3)} \right), \quad (2.49\text{a})
$$

$$
g_{(4)}^{uu} = G^2 \left( \frac{4}{5} M_{ij}^{(1)} M_{ij}^{(1)} + \frac{12}{7} \hat{n}^{ij} M_{ik}^{(1)} M_{jk}^{(1)} + \frac{13}{2} \hat{n}^{ijkl} M_{ij}^{(1)} M_{kl}^{(1)} - \frac{3}{5} M_{ij} M_{ij}^{(2)} \right. \\
\left. \qquad - \frac{36}{7} M_{ij} \hat{n}^{ik} M_{kj}^{(2)} - \frac{387}{8} M_{ij} \hat{n}^{ijkl} M_{kl}^{(2)} \right), \quad (2.49\text{b})
$$

$$
g_{(5)}^{uu} = G^2 \left( \frac{36}{25} M_{ij} M_{ij}^{(1)} + \frac{54}{35} M_{ij} \hat{n}^{ik} M_{kj}^{(1)} - \frac{81}{5} M_{ij} \hat{n}^{ijkl} M_{kl}^{(1)} \right), \quad (2.49c)
$$

$$
g_{(6)}^{uu} = G^2 \left( \frac{3}{5} M_{ij} M_{ij} - \frac{21}{2} M_{ij} M_{kl} \hat{n}^{ijkl} \right), \quad (2.49\text{d})
$$

and

$$
g_{(2)}^{ui} = \left\{ 3GM_{ij}n^j + G^2 \left[ \left( -\frac{33}{20} M_{pq}^{(1)} M_{iq}^{(2)} - \frac{49}{20} M_{iq}^{(1)} M_{pq}^{(2)} - \frac{159}{140} M_{pq} M_{iq}^{(3)} - \frac{159}{140} M_{iq} M_{pq}^{(3)} \right) n^p \right. \right. \\
\left. \left. + \left( -\frac{163}{8} M_{pq}^{(1)} M_{ij}^{(2)} - \frac{17}{4} M_{pi}^{(1)} M_{qj}^{(2)} - \frac{87}{8} M_{pq} M_{ij}^{(3)} - \frac{51}{4} M_{pi} M_{qj}^{(3)} \right) \hat{n}^{pqj} \right] \right\}^T, \quad (2.50\text{a})
$$

$$
g_{(3)}^{ui} = G^2 \left[ \left( \frac{24}{25} M_{iq}^{(1)} M_{pq}^{(1)} - \frac{53}{25} M_{pq} M_{iq}^{(2)} - \frac{113}{25} M_{iq} M_{pq}^{(2)} \right) n^p \right. \\
\left. + \left( \frac{42}{5} M_{pi}^{(1)} M_{qj}^{(1)} - \frac{253}{10} M_{pq} M_{ij}^{(2)} - \frac{89}{5} M_{pi} M_{qj}^{(2)} \right) \hat{n}^{pqj} \right]^T, \quad (2.50\text{b})
$$

$$
g_{(4)}^{ui} = G^2 \left[ \left( -\frac{9}{5} M_{pq} M_{iq}^{(1)} - \frac{1}{5} M_{iq} M_{pq}^{(1)} \right) n^p + \left( -\frac{91}{8} M_{pq} M_{ij}^{(1)} + \frac{11}{8} M_{pi} M_{qj}^{(1)} \right) \hat{n}^{pqj} \right]^T, \quad (2.50c)
$$

$$
g_{(5)}^{ui} = G^2 \left[ -\frac{42}{25} M_{iq} M_{pq} n^p - \frac{51}{5} M_{pi} M_{qj} \hat{n}^{pqj} \right]^T, \quad (2.50\text{d})
$$

where the superscript T refers to transverse projection, i.e., $X_i^T \equiv \perp^{ij} X_j$. For the sub-leading terms in the angular part of the metric, we find

$$
E_{ij} = 2GM_{ij} + G^2 \left[ -\frac{8}{3} M_{ip}^{(1)} M_{jp}^{(2)} - \frac{52}{63} M_{ip} M_{jp}^{(3)} + \hat{n}^{pq} \left( \frac{1}{2} M_{ij}^{(1)} M_{pq}^{(2)} - 8 M_{pi}^{(1)} M_{qj}^{(2)} \right. \right. \\
\left. \left. - 3 M_{pq}^{(1)} M_{ij}^{(2)} - \frac{11}{6} M_{ij} M_{pq}^{(3)} - \frac{22}{3} M_{pi} M_{qj}^{(3)} - \frac{4}{3} M_{pq} M_{ij}^{(3)} \right) \right]^{\text{TT}}, \quad (2.51\text{a})
$$

$$
E_{ij} = G^2 \left[ -2 M_{ip} M_{jp}^{(2)} + \hat{n}^{pq} (6 M_{pi}^{(1)} M_{qj}^{(1)} + M_{pq}^{(1)} M_{ij}^{(1)} - \frac{9}{2} M_{ij} M_{pq}^{(2)} - \frac{57}{2} M_{pi} M_{qj}^{(2)} \right. \\
\left. - \frac{15}{4} M_{pq} M_{ij}^{(2)} \right]^{\text{TT}}, \quad (2.51\text{b})
$$
$$E_{ij}^{(4)} = G^2 \left[ \hat{n}^{pq} \left( \frac{173}{20} M_{ij} M_{pq}^{(1)} - \frac{67}{5} M_{pi} M_{qj}^{(1)} - \frac{97}{20} M_{pq} M_{ij}^{(1)} \right) - \frac{64}{15} M_{ip} M_{jp}^{(1)} \right]^{\text{TT}}, \quad (2.51c)$$

$$E_{ij}^{(5)} = G^2 \left[ \hat{n}^{pq} \left( -\frac{31}{2} M_{pi} M_{qj} + \frac{3}{2} M_{pq} M_{ij} \right) - \frac{11}{3} M_{ip} M_{jp} \right]^{\text{TT}}. \quad (2.51d)$$

We recall that the TT projection is $X_{ij}^{\text{TT}} \equiv \perp_{\text{TT}}^{ijmn} X_{mn}$, with $\perp_{\text{TT}}^{ijmn} = \perp^{m(i}\perp^{j)n} - \frac{1}{2}\perp^{ij}\perp^{mn}$.

The interaction between the mass monopole $M$ and quadrupole $M_{ij}$ leads to tail integrals, which appeared already in the waveform, through the radiative moment $U_{ij}$ in eq. (2.46a), but also in the Bondi aspects $m, N_i$, as well as the subdominant terms ${}_{(n)}g_{uu}$ and ${}_{(n)}E_{ij}$. The complete NU metric for the interaction $M \times M_{ij}$ has been obtained in eqs. (4.11)–(4.13) of [61] in “exact” form, valid for $r$ outside the domain of the source. The expansion of this metric at future null infinity is regular (under our assumption of past stationarity) and can be straightforwardly computed. We must split the moment into a constant piece $M_{ij}(-\mathcal{T})$, where $-\mathcal{T}$ is the finite instant before which the multipole moment is constant, and a dynamical part $M_{ij}(u) - M_{ij}(-\mathcal{T})$, which is zero in the past. The result for the integral entering the $uu$ component of the metric can be found in eq. (4.12) of [61], namely

$$\int_0^{+\infty} dz \frac{M_{ij}(u-z)}{\left(1 + \frac{z}{2r}\right)^2} = 2r M_{ij}(-\mathcal{T}) \quad (2.52a) \\ + \sum_{p=0}^{+\infty} \frac{(-)^p (p+1)}{(2r)^p} \int_0^{+\infty} dz z^p [M_{ij}(u-z) - M_{ij}(-\mathcal{T})].$$

Furthermore, we provide here the regular expansions of the integrals appearing in the $ua$ and $ab$ components:

$$\int_0^{+\infty} dz \frac{5 + \frac{3z}{2r}}{\left(1 + \frac{z}{2r}\right)^3} M_{ij}(u-z) = 8r M_{ij}(-\mathcal{T}) \quad (2.52b) \\ + \sum_{p=0}^{+\infty} \frac{(-)^p (p+1)(p+5)}{(2r)^p} \int_0^{+\infty} dz z^p [M_{ij}(u-z) - M_{ij}(-\mathcal{T})],$$

$$\int_0^{+\infty} dz \frac{18 + \frac{8z}{r} + \frac{z^2}{r^2}}{\left(1 + \frac{z}{2r}\right)^4} M_{ij}(u-z) = 20r M_{ij}(-\mathcal{T}) \quad (2.52c) \\ + \sum_{p=0}^{+\infty} \frac{(-)^p (p+1)(p+3)(p+6)}{(2r)^p} \int_0^{+\infty} dz z^p [M_{ij}(u-z) - M_{ij}(-\mathcal{T})].$$

These expansions confirm that under the assumption of stationarity in the past, the radial expansion is regular to any order. In particular, the latter integral will contribute to the subdominant ${}_{(n)}E_{ij}$ for any $n \ge 3$. Explicitly, the $M \times M_{ij}$ contributions to $m, N_i, H_{ij}^{\text{TT}}$
and $_{(n)}$E$_{ij}$ are given by

$$ m|_{M \times M_{ij}} = 6G^2 M n^{ij} \int_0^{+\infty} dz \left[ \ln\left(\frac{z}{2b_0}\right) + \frac{11}{12} \right] M_{ij}^{(4)}(u-z), \quad (2.53a) $$

$$ N_i|_{M \times M_{ij}} = 12G^2 M \perp_{ik} n^i \int_0^{+\infty} dz \left[ \ln\left(\frac{z}{2b_0}\right) + \frac{11}{12} \right] M_{kj}^{(3)}(u-z), \quad (2.53b) $$

$$ H_{ij}^{\text{TT}}|_{M \times M_{ij}} = 4G^2 M \perp_{ijkl}^{\text{TT}} \int_0^{+\infty} dz \left[ \ln\left(\frac{z}{2b_0}\right) + \frac{11}{12} \right] M_{kl}^{(4)}(u-z), \quad (2.53c) $$

$$ E_{ij} \Big|_{M \times M_{ij}}^{(2)} = G^2 M \perp_{ijkl}^{\text{TT}} \left( M_{kl}^{(1)} + 4 \int_0^{+\infty} \left[ \ln\left(\frac{z}{2b_0}\right) + \frac{11}{12} \right] M_{kl}^{(2)}(u-z) \right), \quad (2.53d) $$

$$ E_{ij} \Big|_{M \times M_{ij}}^{(3)} = 5G^2 M \perp_{ijkl}^{\text{TT}} M_{kl}(-T), \quad (2.53e) $$

$$ E_{ij} \Big|_{M \times M_{ij}}^{(n \ge 4)} = G^2 M \perp_{ijkl}^{\text{TT}} \frac{(-)^n (n-3)(n-1)(n+2)}{2^{n-2}} \\ \times \int_0^{+\infty} dz \, z^{n-4} [M_{kl}(u-z) - M_{kl}(-T)]. \quad (2.53f) $$

## 3 Bondi aspects and charges

In this section, we will discuss gravitational charges in the multipolar expansion of the metric. We will discuss two classes of gravitational charges: 1) those strictly conserved in the generic dynamics of GR in asymptotically flat spacetimes, 2) charges that are time-independent only in non-radiative spacetimes, while they obey a flux-balance equation in the presence of radiation. For the first class, we will discuss the Newman-Penrose charges [67, 68], while for the latter, we will study the celestial charges as defined in [31, 46, 83].

### 3.1 Newman-Penrose charges

The Newman-Penrose (NP) charges [53, 67, 68, 84] are exactly conserved quantities at any u defined at null infinity. In the MPM formalism, assuming that there is no incoming radiation from past null infinity, the metric is entirely determined as a functional of the canonical multipole moments $M_L(u)$, $S_L(u)$, which can have arbitrary time dependence, except for the lowest $\ell = 0, 1$ multipole moments: $M, M_i \equiv G_i = P_i u + K_i$ and $S_i$. Under the stronger assumption of past stationarity, all multipole moments become constants before a given retarded time $u = -T$. Thus, $P_i = 0$ and, furthermore, in the center-of-mass frame that we are using, $M_i = 0$. Therefore, a fully conserved quantity is expected to be a polynomial of the set of Geroch-Hansen multipole moments at spatial infinity given in terms of $M_L(-T)$, $\ell \ge 0$, $S_L(-T)$, $\ell \ge 1$. In what follows, we evaluate the NP charges and check that they do not display quadrupole-quadrupole interactions, while they reproduce known expressions including monopole-quadrupole tail interactions.

The NP charges are defined using the Weyl scalar

$$ \Psi_0 = -C_{\mu\nu\alpha\beta} \ell^\mu m^\nu \ell^\alpha m^\beta, \quad (3.1) $$
in terms of a null tetrad, consisting of the incoming and outgoing real null vectors

$$n = e^{-2\beta}(\partial_u - F\partial_r + U^a\partial_a), \quad \ell = \partial_r, \qquad (3.2)$$

and the complex null vector $m$ and its complex conjugate $\bar{m}$, with

$$m = \frac{1}{r} \left( \zeta^a - \frac{1}{2r} C^a{}_b \zeta^b + \mathcal{O}(r^{-2}) \right) \partial_a + \omega \partial_r. \qquad (3.3)$$

In the above equations, $\beta$, $F$, $U^a$, $\zeta^a$ and $\omega$ are functionals of the Bondi data that are fixed by the normalization property of the tetrad (see e.g., [33, 34, 59, 85]). The pair $(\zeta^a, \bar{\zeta}^a)$ forms a null dyad on the sphere normalized as $\gamma_{ab}\zeta^a\bar{\zeta}^b = 1$. The bulk extension of $m$ is fixed by requiring that it is parallelly transported along the outgoing null vector, i.e. $\ell^\nu \nabla_\nu m^\mu = 0$, as in [68]. One may choose the dyad to be adapted to the stereographic coordinates $(z, \bar{z})$ on the sphere, in which the metric reads $ds^2 = \frac{4}{(1+z\bar{z})^2}dzd\bar{z}$, so that $\zeta^a\partial_a = \frac{1+z\bar{z}}{2}\partial_z$. Residual transformations of the tetrad only include global Type III rotations that covariantly transform $\zeta^a$.<fnref n="9" /> Computing the asymptotic form of this Weyl scalar, we find that $\Psi_0 = \Psi_0^0 r^{-5} + \Psi_0^1 r^{-6} + \mathcal{O}(r^{-7})$, where the leading and subleading coefficients are

$$\Psi_0^0 = 3 \left( E_{ab} - \frac{1}{16} C^2 C_{ab} \right) \zeta^a \zeta^b, \qquad \Psi_0^1 = 6 E_{ab} \zeta^a \zeta^b. \qquad (3.4)$$

The (complex-valued) Newman-Penrose charges are defined, for $m = -2, -1, 0, 1, 2$, as

$$Q_m \equiv \oint_S 2\bar{Y}_{2m} \Psi_0^1, \qquad (3.5)$$

where $sY_{lm}$ denote the spin weighted spherical harmonics of spin weight $s$ (the bar indicates the complex conjugate). They are given for $0 \le s \le l$ in terms of the Geroch-Held-Penrose operator $\partial^{10}$ by

$$s Y_{lm} = \sqrt{\frac{(l-s)!}{(l+s)!}} \partial^s Y_{lm}, \qquad (3.6a)$$

$$\text{and thus } \ _2\bar{Y}_{2m} = \frac{1}{2\sqrt{6}} \bar{\partial}^2 \bar{Y}_{2m} = \frac{1}{2\sqrt{6}} \bar{\zeta}^a \bar{\zeta}^b D_a D_b \bar{Y}_{2m}. \qquad (3.6b)$$

Upon using the decomposition $\zeta^a\bar{\zeta}_b = \frac{1}{2}\delta^a{}_b + \frac{i}{2}\epsilon^a{}_b$ and integrating by parts, we find

$$Q_m = \frac{1}{2} \sqrt{\frac{3}{2}} \oint_S \bar{Y}_{2m} D_a D_b \left( E_{(3)}^{ab} - i \tilde{E}_{(3)}^{ab} \right), \qquad (3.7)$$

$$\bar{\partial}W = \zeta^c \zeta^{a_1} \dots \zeta^{a_n} \bar{\zeta}^{b_1} \dots \bar{\zeta}^{b_m} D_c W_{a_1 \dots a_n b_1 \dots b_m},$$

$$\bar{\partial}W = \bar{\zeta}^c \zeta^{a_1} \dots \zeta^{a_n} \bar{\zeta}^{b_1} \dots \bar{\zeta}^{b_m} D_c W_{a_1 \dots a_n b_1 \dots b_m}.$$

<footnote n="9">⁹This residual symmetry was used in [35] to interpret the gyroscopic memory as a vacuum transition.</footnote>

<footnote n="10">¹⁰Recall that the $\partial$ and $\bar{\partial}$ operators act on a weighted scalar $W = W_{a_1 \dots a_n b_1 \dots b_m} \zeta^{a_1} \dots \zeta^{a_n} \bar{\zeta}^{b_1} \dots \bar{\zeta}^{b_m} D_c W_{a_1 \dots a_n b_1 \dots b_m}$ as</footnote>
where, for any STF tensor $X_{ab}$, we pose $\tilde{X}_{ab} = \epsilon_{ac}X^c_b$. Now, we make a change of basis from spherical harmonics $Y_{2m}$ to STF harmonics $\hat{n}_{ij}$, and we adjust the normalization to define STF NP charges $Q_{ij}$ as

$$Q_{ij} = \oint_S D^a D^b \hat{n}_{ij} \left( E_{ab} - i \tilde{E}_{ab} \right). \qquad (3.8)$$

Given that $D_a n^i = e_a^i$ and $D_a D_b n^i = -\gamma_{ab} n^i$, it can be checked that $D^a D^b \hat{n}_{ij} = 2e^a_{(i}e^b_{j)}$. Therefore, we arrive at

$$Q_{ij} = 2 \oint_S e^a_{(i} e^b_{j)} \left( E_{ab} - i \tilde{E}_{ab} \right) = 2 \oint_S \left( E_{ij} - i \tilde{E}_{ij} \right), \qquad (3.9)$$

where we have resorted to the Cartesian presentation of the Bondi data

$$E_{ij} \equiv e^a_i e^b_j E_{ab}, \qquad \tilde{E}_{ij} \equiv e^a_i e^b_j \tilde{E}_{ab} = n_m \epsilon_{min} E_{nj}. \qquad (3.10)$$

Evaluating the integral in eq. (3.9) amounts to computing the monopole moment of its integrand. However, by explicit computation of the TT projection in eq. (2.51b), we observe that neither $_{(3)}E_{ij}$ nor $_{(3)}\tilde{E}_{ij}$ contain monopole terms. Thus, there is no mass quadrupole-quadrupole contribution to the NP charges. The mass monopole-quadrupole contribution to the NP charges is a direct consequence of eqs. (2.53) obtained from eq. (4.13b) of [61]. Expanding $\perp_{ijkl}^{TT}$ in the expression of $_{(3)}E_{ij}$ in eq. (2.53e) one finds

$$Q_{ij} = 5G^2 M \oint_S \left[ 2M_{ij}(-T) - 2M_{il}(-T)n_j n_l - 2M_{jl}(-T)n_i n_l + n_{ikjl}M_{kl}(-T) \right] = 4G^2 M M_{ij}(-T). \qquad (3.11)$$

The NP charge is indeed conserved for any $u$. This result is consistent with the expressions found by Newman-Penrose [67, 68] and van den Burg [53]. Note that our result is written in the center-of-mass frame in which the mass dipole moment is zero.

## 3.2 From Newman-Unti to Bondi coordinates

In section 2, we obtained the radiative metric corresponding to selected multipole interactions in Newman-Unti (NU) gauge. However, many results on the asymptotic structure of asymptotically flat spacetimes are known in Bondi gauge instead. In the following, we complete the map from NU to Bondi gauge outlined in appendix A of [61].

The NU and Bondi coordinates are related by a change of radius, namely $r_{NU}$ in NU coordinates and $r_B$ in Bondi coordinates, with $u$ and angle coordinates $\theta^a$ unchanged. In Bondi gauge, the spherical metric takes the form [31]

$$g_{ab} = r_{\text{B}}^2 \sqrt{1 + \frac{C_{cd}^{\text{B}} C_{\text{B}}^{cd}}{2r_{\text{B}}^2}} \gamma_{ab} + r_{\text{B}} C_{ab}^{\text{B}}, \qquad C_{ab}^{\text{B}} \equiv e^i_{(a} e^j_{b)} \left( H_{ij}^{\text{TT}} + \sum_{n=2}^{+\infty} \frac{1}{r_{\text{B}}^n} E_{ij}^{(n)} \right). \qquad (3.12)$$

The relation between the two radii follows from

$$r_{\text{B}}^4 = \left. \frac{\det g_{ab}}{\det \gamma_{ab}} \right|_{\text{NU}}. \qquad (3.13)$$
We shall parametrize the NU metric as

$$g_{ab} = (r_{\text{NU}}^2 + W(r_{\text{NU}}, \theta^a))\gamma_{ab} + r_{\text{NU}} C_{\langle ab \rangle}^{\text{NU}}(r_{\text{NU}}, \theta^a), \quad (3.14)$$

where $W = \mathcal{O}(G^2)$ and $C_{\langle ab \rangle}^{\text{NU}} = \mathcal{O}(G)$. Then, eq. (3.13) yields

$$r_B = \left[ (r_{\text{NU}}^2 + W)^2 - \frac{r_{\text{NU}}^2}{2} C_{\langle ab \rangle}^{\text{NU}} C_{\text{NU}}^{\langle ab \rangle} \right]^{1/4} = r_{\text{NU}} + \frac{1}{2r_{\text{NU}}} \left( W - \frac{1}{4} C_{\langle ab \rangle}^{\text{NU}} C_{\text{NU}}^{\langle ab \rangle} \right) + \mathcal{O}(G^4), \quad (3.15)$$

where $W$ and $C_{\langle ab \rangle}^{\text{NU}} C_{\text{NU}}^{\langle ab \rangle}$ can be evaluated either at radius $r_{\text{NU}}$ or $r_B$, neglecting $\mathcal{O}(G^4)$ corrections. At large radius, it can be shown, using $g_{\text{ur}}^{\text{B}}(u, r_B, \theta^a) = -\partial r_{\text{NU}} / \partial r_B$ and $g_{ab}^{\text{B}}(u, r_B, \theta^a) = g_{ab}^{\text{NU}}(u, r_{\text{NU}}, \theta^a)$, combined with the $1/r_B$ expansion of $g_{\text{ur}}^{\text{B}}$, that

$$W = \frac{1}{8} C_{ab} C^{ab} + \mathcal{O}(r_{\text{NU}}^{-1}), \quad C_{\langle ab \rangle}^{\text{NU}} = C_{ab} + \mathcal{O}(r_{\text{NU}}^{-1}), \quad (3.16)$$

and therefore

$$r_B = r_{\text{NU}} - \frac{1}{16r_{\text{NU}}} C_{ab} C^{ab} + \mathcal{O}(r_{\text{NU}}^{-2}), \quad (3.17)$$

which reproduces eq. (61b) of [61]. Including all orders in the radius, we may pose

$$C_{\langle ab \rangle}^{\text{NU}} = C_{ab} + e_{\langle a}^{i} e_{b \rangle}^{j} \sum_{n=2}^{+\infty} r_{\text{NU}}^{-n} E_{ij} \quad (3.18)$$

We see that the relation (3.15) can be computed exactly at quadratic order, knowing $W$ from the quadratic NU metric, as well as the linear part of $C_{ab}$ and $\langle k \rangle E_{ij}^{11}$ given, respectively, by eqs. (3.18) and (3.15c) in [61]. Keeping only linear terms in $M$ and $M_{ij}$ together with BMS diffeomorphisms, we first get

$$C_{\langle ab \rangle}^{\text{NU}} = -2GD_{\langle a}D_{b \rangle}f + 2Ge_{\langle a}^{i}e_{b \rangle}^{j} \left( M_{ij}^{(2)} + \frac{1}{r_{\text{NU}}^{2}}M_{ij} \right), \quad (3.19)$$

where $f = T(\theta^a) + \frac{u}{2}D_a Y^a$ results from a BMS transformation. In the absence of BMS transformation, we then find

$$r_{\text{NU}} = r_B - \frac{1}{2r_B} \left[ W - G^2 \perp_{TT}^{ijkl} \left( M_{ij}^{(2)} + \frac{1}{r_B^2} M_{ij} \right) \left( M_{kl}^{(2)} + \frac{1}{r_B^2} M_{kl} \right) \right] + \mathcal{O}(G^4). \quad (3.20)$$

The Bondi mass aspect can be read in Bondi gauge from $g_{uu} = -1 + 2m/r + \mathcal{O}(r^{-2})$. The Bondi angular momentum aspect $N_a$ can be read in Bondi gauge from$^{12}$

$$g_{ua} = \frac{D^b C_{ab}}{2} + \frac{2}{3r_B} \left[ N_a - \frac{3}{32} D_a (C_{bc} C^{bc}) \right] + \mathcal{O}(r_B^{-2}). \quad (3.21)$$

<footnote>$^{11}$Note that in [61], the linear part of $\langle k \rangle E_{ij}$ is denoted by $Q_k^{ij}$ for $k \ge 2$.</footnote>

<footnote>$^{12}$Given the non-universal conventions in the literature, we provide here a dictionary to some other references: $N_a$ as defined in [31, 58] is equal to $N_a$ here; the covariant momentum $P_a$ as defined in [85, 86] is equal to $N_a$ here, and $N_a$ as defined in [37] is equal to $N_a - \frac{1}{4}C_{ab}D_c C^{cb} - \frac{3}{32}\partial_a(C_{cd}C^{cd})$ here.</footnote>
Since the Bondi radius $r_B$ differs from the Newman-Unti radius $r_{NU}$ by terms of order $G^2$ while $H_{ij}^{TT}$ and $E_{ij}$ are of order $G$, we simply have
(n)

$$E_{ij}^{B} = E_{ij} + \mathcal{O}(G^3). \qquad (3.22)$$

This completes the map of all initial data in Bondi gauge in terms of Newman-Unti initial data at order $G^2$ for the interactions considered in this work.

### 3.3 Flux-balance laws of the dressed $n \le 3$ Bondi aspects

In the Bondi expansion, Einstein's equations lead to time evolution equations for all $n \ge 0$ Bondi aspects. The Bondi aspects can be further "dressed" by adding suitable non-linear combinations of Bondi data such that the dressed Bondi aspects are identically conserved when the news vanishes. In this section, we recall the definitions of the dressed $n = 0, 1, 2$ Bondi aspects and propose a new definition of the dressed $n = 3$ Bondi aspect following the construction of [31]. We do not consider the case where $n > 3$.

It will be useful to first introduce the gravitational electric-magnetic (or mass-current) duality covariant mass [40, 41, 45]

$$m_{ab} = m \gamma_{ab} + \frac{1}{2} D_{[a} D^{c} C_{b]c}. \qquad (3.23)$$

Mass multipole moments are annihilated by the differential operator entering the second term of eq. (3.23). As a result, only current multipole moments appear in eq. (3.23). The flux-balance laws of the first $n = 0, 1, 2, 3$ Bondi aspects are given by

$$\partial_u m = -\frac{1}{8}N^2 + \frac{1}{4}D_a D_b N^{ab}, \qquad (3.24a)$$

$$\partial_u N_a = D^b m_{ab} + \frac{1}{4} (N^{bc} D_c C_{ab} + 3 C_{ab} D_c N^{bc}), \qquad (3.24b)$$

$$\partial_u E_{(2)}_{ab} = \frac{1}{4}(CN)C_{ab} + \frac{1}{2}m_{ac}C_b^c + \frac{1}{3}D_{\langle a}N_{b\rangle}, \qquad (3.24c)$$

$$\partial_u E_{(3)}_{ab} = D_0 E_{(2)}_{ab} + D^c \left[ \left( \frac{1}{4} D_e C^{de} C_{d\langle a} - \frac{3}{32} D_{\langle a} C^2 + \frac{5}{32} C^2 D_{\langle a} - \frac{1}{3} N_{\langle a} \right) C_{b\rangle c} \right], \qquad (3.24d)$$

where we use the compact notation $C^2 \equiv C_{ab}C^{ab}$, $(CN) \equiv C_{ab}N^{ab}$, $N^2 \equiv N_{ab}N^{ab}$, and where $D_0 \equiv -\frac{1}{4}(\Delta + 2)$. The $n = 1$ flux-balance law matches [58]. The evolution of the $n = 2$ Bondi aspect in vacuum was first derived in [53] and agrees with that obtained in [28]. The $n = 2, 3$ expressions are taken from [31] (see alternative derivations in [51, 53, 87]).

The $n = 0$ Bondi aspect $m$ is just the Bondi mass aspect that permits defining the supermomenta

$$\mathcal{P}_L = \oint_S m \hat{n}_L, \qquad (3.25)$$

i.e., the canonical charges associated with supertranslations $T = T_L \hat{n}_L$. Its explicit expression, in the post-Minkowskian approximation, restricted to linear terms, tails and quadrupole-quadrupole interactions can be obtained from eq. (2.47a). For the monopole
$\ell = 0$, the energy is defined as $G\mathcal{P}_0(u) \equiv \mathcal{E}(u) = M - E^{\text{GW}}(u)$, where $E^{\text{GW}}$ was defined in eq. (2.15).

The $n=1$ Bondi aspect $N_a$, usually referred to as the angular momentum aspect, can be supplemented or “dressed” with suitable terms such that its retarded time derivative vanishes when the news vanishes. The dressed $n=1$ [60, 69] and $n=2$ Bondi aspects [31] read as<fnref n="13" />

$$ \mathcal{N}_a = N_a - \frac{1}{16} D_a C^2 - \frac{1}{4} C_a^b D^c C_{bc} - u D^b m_{ab}, \qquad (3.26a) $$

$$ \mathcal{E}_{ab}^{(2)} = E_{ab}^{(2)} - \frac{u}{2} C_{(a}^{c} m_{b)c} - \frac{u}{3} D_{\langle a} N_{b\rangle} + \frac{u^2}{6} D_{\langle a} D^{c} m_{b\rangle c}. \qquad (3.26b) $$

They obey the flux-balance laws

$$ \partial_u \mathcal{N}_a = \frac{1}{4} N^{bc} D_b C_{ca} + \frac{1}{2} C_{ab} D_c N^{bc} - \frac{1}{4} N_{ab} D_c C^{bc} - \frac{1}{8} D_a (C N) + \frac{u}{8} D_a N^2 - \frac{u}{2} D_c D_{\langle a} D_{b\rangle} N^{bc}, \qquad (3.27a) $$

$$ \partial_u \mathcal{E}_{ab}^{(2)} = \frac{1}{4}(CN)C_{ab} - \frac{u}{2}\partial_u(C_a^c m_{bc}) - uD_{\langle a} \left( \frac{1}{12}N^{cd}D_{\underline{d}}C_{\underline{cb}\rangle} + \frac{1}{4}C_{b\rangle c}D_d N^{cd} \right) - \frac{u^2}{48}D_{\langle a}D_{b\rangle}N^2 + \frac{u^2}{12}\text{STF}_{ab}[D_a D_c D_{\langle b}D_{d\rangle}N^{cd}]. \qquad (3.27b) $$

The $n=1$ dressed aspect allows defining the super-Lorentz conserved charges

$$ \mathcal{J}_L = -\frac{1}{2} \oint_S \epsilon^{ab} \partial_b \hat{n}_L \mathcal{N}_a, \qquad \mathcal{K}_L = \frac{1}{2} \oint_S \gamma^{ab} \partial_b \hat{n}_L \mathcal{N}_a \qquad (3.28) $$

for arbitrary $\text{Diff}(S^2)$ generators $Y^a = -\epsilon^{ab} \partial_b \hat{n}_L + \gamma^{ab} \partial_b \hat{n}_L$. The $n=2$ dressed aspect is associated with the $n=2$ “celestial” charges $\oint_S \mathcal{E}_{ab}(D^a D^b S^+ + \epsilon_{ac} D^b D^c S^-)$, defined for arbitrary scalars $S^{\pm}(\theta^a)$ on the 2-sphere.

The dressed $n=3$ Bondi aspect, which is conserved in non-radiative regions, was introduced in eq. (4.43) of [31].<fnref n="14" /> It can be written in terms of the duality-covariant quantity $m_{ab}$ as it should come from gravitational electric-magnetic duality [83]. The (corrected) dressed $n=3$ Bondi aspect reads as

$$ \mathcal{E}_{ab}^{(3)} = E_{ab}^{(3)} - u \left\{ \mathcal{D}_0 E_{ab}^{(2)} + D^c \left[ \left( \frac{1}{4} D_e C^{de} C_{d\langle a} - \frac{3}{32} D_{\langle a} C^2 + \frac{5}{32} C^2 D_{\langle a} - \frac{1}{3} N_{\langle a} \right) C_{b\rangle c} \right] \right\} + \frac{u^2}{2} \left[ -\frac{1}{3} D^c (D^d m_{d\langle a} C_{b\rangle c}) + \frac{1}{2} \mathcal{D}_0 (m_{ac} C_b^c) + \frac{1}{3} \mathcal{D}_0 D_{\langle a} N_{b\rangle} \right] - \frac{u^3}{18} \mathcal{D}_0 D_{\langle a} D^c m_{b\rangle c}. \qquad (3.29) $$

It obeys the flux-balance law

$$ \partial_u \mathcal{E}_{ab}^{(3)} = \mathcal{F}_{ab}^{(3)} - \frac{u^3}{36} \mathcal{D}_0 \text{STF}_{ab} (D_a D_c D_{\langle b} D_{d\rangle} N^{cd}), \qquad (3.30) $$

<footnote n="13">¹³Equation (2.10) of [83] should read as eq. (3.26b).</footnote>

<footnote n="14">¹⁴However, three typos arose in writing this expression as can be deduced from logically following the algorithm described in [31]; see the erratum of [31]. In eq. (4.40c), the factor 1/12 should read 1/6; in eq. (4.39a), the factor -3/4 should be -3/8 and the factor 5/4 should be 5/8.</footnote>
where the flux $\mathcal{F}_{ab}$ can be written, using the flux-balance laws (3.24) and the identity $C_{c(a}N_{b)}^c = \frac{1}{2}\gamma_{ab}(CN)$, as

$$ \begin{aligned} \mathcal{F}_{ab} = & -\frac{u}{4}\mathcal{D}_0((CN)C_{ab}) + \frac{u}{3}D^c \left[ N_{\langle a}N_{b\rangle c} + \left(\frac{1}{4}N^{de}D_eC_{d\langle a} + \frac{3}{4}D_eN^{de}C_{d\langle a}\right)C_{b\rangle c} \right] \\ & + \frac{u}{4}\partial_u D^c \left[ \frac{3}{8}D_{\langle a}(C^2C_{b\rangle c}) - (D_eC^{de})(C_{d\langle a}C_{b\rangle c}) - C^2D_{\langle a}C_{b\rangle c} \right] \\ & + \frac{u^2}{24}\mathcal{D}_0D_{\langle a}(N^{de}D_{\underline{e}}C_{b\rangle d} + 3C_{b\rangle d}D_eN^{de}) + \frac{u^2}{2}\partial_u \left[ -\frac{1}{3}D^c(D^d m_{d\langle a}C_{b\rangle c}) + \frac{1}{2}\mathcal{D}_0(m_{ac}C_b^c) \right] \\ & + \frac{u^3}{144}\mathcal{D}_0D_{\langle a}D_{b\rangle}N^2. \end{aligned} \quad (3.31) $$

As expected, it vanishes when the news does. We cross-checked with a computer code that the flux-balance laws for (2) $\mathcal{E}_{ab}$ and (3) $\mathcal{E}_{ab}$ are indeed obeyed for the $M_{ij} \times M_{ij}$ interactions.

**Use of dimensional identities.** In order to verify the flux-balance laws in practice (and, more generally, to verify the relations between different pieces of the radial expansion of the metric), we must make use of so-called dimensional identities. Indeed, in some cases, the difference between the left and right-hand sides of the balance equations, although zero, is not manifestly zero. Indeed, given a tensor of rank greater than four in dimension three, say $T_{ijkl}$, an identity such as $T_{ijkl} \equiv 0$, referred to as dimensional identities, is not explicitly “apparent”.

In our problem, we meet expressions like $M_{ij}^{(n_1)} M_{kl}^{(n_2)} \hat{n}_{mL} e_a^p e_b^q$ which do contain at least four free indices whose antisymmetrization will not trivially yield zero, e.g., the indices $i, k, m$ and $p$ in this example. Moreover, the number of free indices can be reduced by contracting pairs lying outside the antisymmetrization operator $[\cdots]$. In particular, it is possible to construct a rank 2 dimensional identity from the latter monomial, e.g., by contracting the pairs $\{i, j\}, \{l, m\}, \{p, q\}, \{k, L = n\}$ (assuming that the length of $L$ is 1). The dimensional identities produced in that way are thus

$$ I_{ab}^{(n_1, n_2)} = -2G^2 M_{j[i}^{(n_1)} M_{k\underline{l}]}^{(n_2)} \hat{n}_{m\underline{n}} e_a^{p]} e_b^q \delta^{ij} \delta^{lm} \delta_{pq} \delta^{kn} = 0. \quad (3.32) $$

They are used in several instances of our calculations. For example, we find that the TT projection of the piece of $g_{ab}$ that is proportional to $r^0$ is

$$ 8 I_{ab}^{(1,3)} + 2 I_{ab}^{(0,4)} = 0. \quad (3.33) $$

The differences of both sides of the balance equations (3.24c) for (2) $E_{ab}$ and (3.24d) for (3) $E_{ab}$ reduce, respectively, to the identities

$$ 8 I_{ab}^{(1,3)} + 2 I_{ab}^{(0,4)} - 6 I_{ab}^{(2,2)} = 0, \quad (3.34a) $$

$$ 21 I_{ab}^{(1,2)} + 23 I_{ab}^{(0,3)} = 0. \quad (3.34b) $$
Finally, subtracting the right-hand sides from the left-hand sides of the balance equa-
tions (3.27b) and (3.30) for the dressed aspects $^{(2)}$E$_{ab}$ and $^{(3)}$E$_{ab}$ , respectively, lead to

$$
8 I_{ab}^{(1,3)} + 2 I_{ab}^{(0,4)} - 6 I_{ab}^{(2,2)} = 0, \qquad (3.35a)
$$

$$
21 I_{ab}^{(1,2)} + 23 I_{ab}^{(0,3)} - 4u I_{ab}^{(1,3)} - u I_{ab}^{(0,4)} + 3u I_{ab}^{(2,2)} = 0. \qquad (3.35b)
$$

Hence thanks to the identities (3.32), we have been able to verify the above claims. A way
to by-pass the use of dimensional identities is to specify and expand all the components of
vectors and tensors in a given chosen basis.

3.4 The $n \le 3$ celestial charges at $G^2$ order

In the linear theory, the dressed $n \ge 0$ Bondi aspects have been shown to provide the
complete set of conserved charges for asymptotically flat spacetimes admitting a Bondi
expansion [83]. In this section, we will explicitly write the dressed $n \le 3$ Bondi aspects in
the post-Minkowskian approximation, with terms linear in $M$ and $M_{ij}$ and $M \times M_{ij}$ and
$M_{ij} \times M_{kl}$ interactions, and compute the corresponding charges.

For convenience, we shall express the Bondi data in terms of the corresponding Carte-
sian transverse tensors

$$
\mathcal{N}_i \equiv e^a{}_i \mathcal{N}_a, \qquad \mathcal{E}_{ij} = e^a{}_i e^b{}_j \mathcal{E}_{ab}. \tag{3.36}
$$

Summing up the different contributions in eqs. (3.26) and (3.29) [with the undressed quantities therein given by eqs. (2.47b), (2.51a), and (2.51b), respectively], we obtain through computer calculation the explicit expressions for $\mathcal{N}_i$, $^{(2)}\mathcal{E}_{ij}$ and $^{(3)}\mathcal{E}_{ij}$. Note that $^{(n)}\mathcal{E}_{ij}$ is traceless with respect to $\perp_{ij}$ because $^{(n)}\mathcal{E}_{ab}$ is traceless with respect to $\gamma_{ab}$. Given the length of the resulting expressions for (3.36), we only provide their expressions for quadrupole-quadrupole interactions in appendix B.

We are now in the position to compute the contributions up to order $G^2$ to the Bondi supermomenta, super-angular momenta and super-center-of-mass defined in eqs. (3.25) and (3.28) as well as the two celestial charges $Q_{2,L}^{\pm}$ and $Q_{3,L}^{\pm}$ defined, respectively, from the $n=2$ and $n=3$ Bondi aspects. The latter charges, which are STF with respect to their indices $L$, are defined in [83] by

$$
Q_{n,L}^{+} \equiv \oint_S \mathcal{E}_{(n)}^{ab} D_a D_b \hat{n}_L, \qquad Q_{n,L}^{-} \equiv \oint_S \mathcal{E}_{(n)}^{ab} \epsilon_{ac} D_b D^c \hat{n}_L. \qquad (3.37)
$$

For all terms considered, $J_L = 0 = Q_{n,L}^-$ for all $L$ and $n \ge 2$ because the integrand is
parity odd. For terms involving the sole quadrupole moment $M_{ij}$ at linear and quadratic
orders and tails of the form $M \times M_{ij}$, the explicit results are

$$ \mathcal{E} = M - \frac{1}{5} G \int_{-\infty}^{u} M_{ij}^{(3)} M_{ij}^{(3)}, \qquad (3.38a) $$

$$ \mathcal{P}_{ij} = \frac{2}{5} G M_{ij}^{\text{rad}(2)} + \frac{2}{35} G^2 [-2M_{ik}^{(2)} M_{jk}^{(3)} - 5M_{ik}^{(1)} M_{jk}^{(4)} + M_{ik} M_{kj}^{(5)}]^{\text{STF}}, \qquad (3.38b) $$

$$ \mathcal{P}_{ijkl} = \frac{1}{105} G^2 [-34M_{ij}^{(2)} M_{kl}^{(3)} - 21M_{ij}^{(1)} M_{kl}^{(4)} - 7M_{ij} M_{kl}^{(5)}]^{\text{STF}}, \qquad (3.38c) $$

$$ \mathcal{K}_{ij} = \frac{6}{5} G M_{ij}^{\text{rad}(1)} + \frac{6}{35} G^2 [-M_{ik}^{(2)} M_{jk}^{(2)} - 6M_{ik}^{(1)} M_{jk}^{(3)} + M_{ik} M_{jk}^{(4)}]^{\text{STF}} \qquad (3.38d) $$

$$ \mathcal{K}_{ijkl} = -\frac{2}{21} G^2 [9M_{ij}^{(2)} M_{kl}^{(2)} + 14M_{ij}^{(1)} M_{kl}^{(3)} + 7M_{ij} M_{kl}^{(4)}]^{\text{STF}}, \qquad (3.38e) $$

$$ Q_{2,ij}^{+} = \frac{8}{5} G \left( M_{ij}^{\text{rad}} - u M_{ij}^{\text{rad}(1)} + \frac{1}{2} u^2 M_{ij}^{\text{rad}(2)} \right) \\ + \frac{8}{35} G^2 \left[ -7 M_{ik}^{(1)} M_{jk}^{(2)} + M_{ik} M_{jk}^{(3)} + u \left( 7 M_{ik}^{(2)} M_{jk}^{(2)} + 6 M_{ik}^{(1)} M_{jk}^{(3)} - M_{ik} M_{jk}^{(4)} \right) \right. \\ \left. - \frac{1}{2} u^2 \left( 2 M_{ik}^{(2)} M_{jk}^{(3)} + 5 M_{ik}^{(1)} M_{jk}^{(4)} - M_{ik} M_{jk}^{(5)} \right) \right]^{\text{STF}}, \qquad (3.38f) $$

$$ Q_{3,ij}^{+} = 4G^2 MM_{ij}(-T), \qquad (3.38g) $$

where $M_{ij}^{\text{rad}}(u)$ contains the tail interactions of the form $M \times M_{ij}$:

$$ M_{ij}^{\text{rad}}(u) = M_{ij}(u) + 2GM \int_{0}^{+\infty} dz \left[ \ln\left(\frac{z}{2b_0}\right) + \frac{11}{12} \right] M_{ij}^{(2)}(u-z) + \mathcal{O}(G^2). \qquad (3.39) $$

These charges were already computed at $\mathcal{O}(G)$ in [83]. In particular, the set of charges for $2 \le l \le n-1$, obeying memory-less flux-balance laws, were shown to be vanishing at the linear order. Including the quadratic contributions, the first non-trivial such charge $Q_{3,ij}^{+}$ gets a constant value, given by the product between the ADM mass and the quadrupole moment at the early time $-T$ when the system is assumed to be stationary. In fact, $Q_{3,ij}^{+}$ exactly matches with the Newman-Penrose charges as demonstrated in eqs. (3.8) and (3.11).

For the complementary set of charges with $l \ge n$, obeying the memory-full flux-balance laws, we have computed ${}_{(2)}\mathcal{E}^{ab}$ and ${}_{(3)}\mathcal{E}^{ab}$ in eqs. (3.26b) and (3.29) explicitly. For quadrupole-quadrupole interactions, the non-vanishing charges can only have 0, 2 or 4 free indices, since they are traceless quantities only built from $M_{ij}M_{kl}$, with $u$ derivatives, $u$-dependent factors and contractions of indices. The remaining non-vanishing charges for
quadrupole-quadrupole interactions are

$$Q_{2,ijkl}^{+} = G^2 \left[ -4(M_{ij}^{(1)}M_{kl}^{(2)} + M_{ij}M_{kl}^{(3)}) + 4u(M_{ij}^{(2)}M_{kl}^{(2)} + 2M_{ij}^{(1)}M_{kl}^{(3)} + M_{ij}M_{kl}^{(4)}) \right. \\ \left. + u^2 \left( -\frac{68}{7}M_{ij}^{(2)}M_{kl}^{(3)} - 6M_{ij}^{(1)}M_{kl}^{(4)} - 2M_{ij}M_{kl}^{(5)} \right) \right]^{\text{STF}}, \quad (3.40a)$$

$$Q_{3,ijkl}^{+} = G^2 \left[ \frac{8}{3} M_{ij}^{(1)} M_{kl}^{(1)} - 14 M_{ij} M_{kl}^{(2)} + u \left( \frac{26}{3} M_{ij}^{(1)} M_{kl}^{(2)} + 14 M_{ij} M_{kl}^{(3)} \right) \right. \\ \left. + u^2 \left( -\frac{13}{3} M_{ij}^{(2)} M_{kl}^{(2)} - 7(2M_{ij}^{(1)} M_{kl}^{(3)} + M_{ij} M_{kl}^{(4)}) \right) \right. \\ \left. + u^3 \left( \frac{34}{3} M_{ij}^{(2)} M_{kl}^{(3)} + 7M_{ij}^{(1)} M_{kl}^{(4)} + \frac{7}{3} M_{ij} M_{kl}^{(5)} \right) \right]^{\text{STF}}. \quad (3.40b)$$

The stationary limit is obtained straightforwardly.

## 4 Conclusions

The asymptotic Bondi-Sachs formalism is a convenient setup to study various gravitational-wave observables, such as fluxes of conserved quantities or the canonical structure of radiative phase space in terms of the Bondi shear. In this setup, the Bondi shear, which characterizes the gravitational-wave strain, is not constrained by Einstein's equations and is thus free data. However, to obtain the Bondi shear for a given source, one needs to combine it with other wave generation methods. This paper, following [61], continues the programme to incorporate results from the PN/MPM formalism into the Bondi-Sachs setup in a systematic way. This analysis allows us to infer several interesting properties of the spacetime geometry; in particular, 1) we highlight its global features, such as the peeling property or the late-time behaviour of the waveform, including memory and tail effects, and 2) we check the conservation of asymptotic BMS charges and their generalizations.

In section 2, we retrace the computation of hereditary terms in the MPM formalism [5]. We also apply the algorithm introduced in [61] to transform into NU gauge the sector of the 2PM metric containing non-linear memory and secular losses. This first result is presented in eqs. (2.26) whereas the asymptotic data, consisting of the Bondi mass and angular momentum aspects and the Bondi shear, are written in eqs. (2.28). The second and main result of this work is the NU metric at 2PM order, limited to multipole interactions of the type monopole-quadrupole $M \times M_{ij}$ and quadrupole-quadrupole $M_{ij} \times M_{ij}$, including tail terms (already obtained in [61]) as well as non-linear memory and mass-loss terms. The NU metric for such interactions is given in eqs. (2.44), where the Bondi mass and angular momentum aspects are reported in eqs. (2.47) [supplemented by the tail contributions, resp., of the first two eqs. (2.53)], while the Bondi shear is written in eqs. (2.45)–(2.46) [supplemented by the tail contribution of the third equation in eqs. (2.53)]. Sub-leading contributions to the $uu$, $ui$ and $ij$ components of the NU metric are displayed, resp., in eqs. (2.49)–(2.50)–(2.51), being the latter supplemented by tail contributions to the $ij$ components in eqs. (2.53). We also remark the explicit proof of the re-summation of
the infinite multipole series in eq. (2.34), which has seldomly appeared in the previous literature, and some special treatment of hereditary terms explained in appendix A.

In section 3, based on the previous section, we discuss gravitational charges and their time evolution in the presence of radiation. One non-trivial test of our results is the explicit check (limited to the specific multipole interactions considered in this work) that Newman-Penrose charges are conserved; they turn out to be proportional to the ADM mass times the initial mass quadrupole moment. A third outcome is the (corrected) derivation of the dressed $n = 3$ Bondi aspect in eq. (3.29) and its flux in eq. (3.31). The flux-balance laws for $n = 2, 3$ were checked with the help of a computer code up to quadrupole-quadrupole interactions. Finally, we explicitly write the first $n \le 3$ Bondi charges, i.e., supermomenta, super-Lorentz and $n = 2, 3$ celestial charges, in eqs. (3.38) and eqs. (3.40) in terms of the mass quadrupole interactions.

## Acknowledgments

G.C. is Senior Research Associate of the F.R.S.-FNRS and acknowledges support from the FNRS research credit J.0036.20F and the IISN convention 4.4503.15. The work of R.O. is supported by the Région Île-de-France within the DIM ACAV⁺ project SYMONGRAV (Symétries asymptotiques et ondes gravitationnelles). A.S. is supported by a Royal Society University Research Fellowship. L.B. and R.O. acknowledge support from the Partenariat Hubert Curien within the Barrande mobility programme (project number 46771VC).

# A Treatment of hereditary integrals

## A.1 Radial integration at fixed retarded time

In order to systematize the construction of the Newman-Unti metric, we develop general formulae for handling the hereditary tail integrals entering the $M_{ij} \times M_{ij}$ harmonic metric. At quadratic order, they all take the form<fnref n="15" />

$$I_m(t, r) = \int_{1}^{+\infty} dx \, Q_m(x) F(t - rx). \qquad (A.1)$$

Here $Q_m(x)$, $m \in \mathbb{N}$, is the Legendre function of the second kind (with branch cut singularity from $-\infty$ to $1$), and $F(u)$ is a quadratic product of time derivatives of quadrupole moments, which vanishes when $u \le -T$, so that the integration range is actually finite. Explicit expressions of the Legendre function are given below in (A.12) and (A.16). They show that the only possible convergence issue concerns the bound $x = 1$, but the convergence is actually guaranteed by the local behaviour $Q_m(x) \sim -\frac{1}{2} \ln(x-1)$ when $x \to 1^+$.

The perturbation equations (2.22) we must solve in our approach are of the type $k^\mu \partial_\mu J_m = r^{-k} I_m$. Let us first treat the case $k = 0$, and see later how we deal with other cases. We are thus looking for a function $J_m(t, r)$ such that

$$\left. \frac{\partial J_m(u+r, r)}{\partial r} \right|_{u=\text{const}} = I_m(t, r). \qquad (A.2)$$

<footnote n="15">¹⁵For simplicity in this appendix we denote the coordinates by $(t, r)$, and $u \equiv t - r$, although they are really meant to be the harmonic coordinates we start with in our calculation.</footnote>
The problem we meet with the hereditary term (A.1) is that the integral $J_m$ formally satisfying eq. (A.2), holding $u = \text{const}$, reads $-\int_1^{+\infty} dx F^{(-1)}(t-rx)Q_m(x)/(x-1)$, plus a possible integration constant, irrelevant for the present discussion. Due to the factor $(x-1)^{-1}$, the integral diverges since $Q_m(x)/(x-1)$ is no longer locally integrable at the bound $x = 1^+$. A way around is to isolate the non-integrable part of $Q_m(x)/(x-1)$ modulo a well-behaved function. In the previous paper [61], we tackled the problem only on a "case-by-case" basis. Here, we present a more general method, based on the following formula for the Legendre function, which we shall prove at the end of this subsection:

$$Q_m(x) - Q_0(x) + H_m = (x-1) \sum_{j=0}^{m-1} c_{mj} Q_j(x), \qquad (A.3a)$$

$$\text{with } c_{mj} \equiv (2j+1)(H_m - H_j), \qquad (A.3b)$$

where $H_m = \sum_{k=1}^m k^{-1}$ denotes the usual harmonic number and we recall that $Q_0(x) = 1/2 \ln[(x+1)/(x-1)]$. Plugging $Q_m(x)$ as deduced from (A.3a) into $I_m$, we readily obtain

$$I_m = \sum_{j=0}^{m-1} c_{mj} \int_1^{+\infty} dx (x-1) Q_j(x) F(t-rx) + \frac{1}{r} \int_0^{+\infty} d\tau \left[ \frac{1}{2} \ln \left( 1 + \frac{2r}{\tau} \right) - H_m \right] F(u-\tau). \qquad (A.4)$$

There is now an explicit factor $x-1$ in the integrand of the first term, so it can be directly integrated and the result will be a sum of hereditary integrals of the same structure as (A.1). As for the second term, where we have posed $x = 1 + \frac{\tau}{r}$, it yields upon integration a dilogarithm function $\text{Li}_2(z) \equiv -\int_0^z ds \ln(1-s)/s$. Hence we obtain

$$J_m = - \sum_{j=0}^{m-1} c_{mj} \int_1^{+\infty} dx Q_j(x) F^{(-1)}(t-rx) \\ - \int_0^{+\infty} d\tau \left[ \frac{1}{2} \text{Li}_2 \left( -\frac{2r}{\tau} \right) + H_m \ln \left( \frac{r}{r_0} \right) \right] F(u-\tau), \qquad (A.5a)$$

plus a possible constant with respect to $r$. Here, $F^{(-1)}$ is the time-antiderivative of $F$, and $r_0$ is an arbitrary integration constant which can be chosen to be the Hadamard regularization scale of the MPM formalism. Alternatively, by decomposing the logarithm arising in eq. (A.4) as $\ln[1+\tau/(2r)] - \ln[\tau/(2r)]$, coming back to the original variable $x$ in the integral of the first term and, finally, integrating with respect to $r$ for constant $u$ with the same procedure as for $\int_1^{+\infty} dx (x-1)Q_j(x)F(t-rx)$, we get

$$J_m = - \sum_{j=0}^{m-1} c_{mj} \int_1^{+\infty} dx Q_j(x) F^{(-1)}(t-rx) - \frac{1}{2} \int_0^{+\infty} \frac{d\tau}{\tau} \ln \left( 1 + \frac{\tau}{2r} \right) F^{(-1)}(u-\tau) \\ + \frac{1}{4} \int_0^{+\infty} d\tau \ln^2 \left( \frac{\tau}{2r} \right) F(u-\tau) - H_m \ln \left( \frac{r}{r_0} \right) F^{(-1)}(u). \qquad (A.5b)$$

This expression makes explicit the appearance of $\ln r$ in the far-zone expansion, assuming that $F(u)$ and $F^{(-1)}(u)$ identically vanish in the remote past. To recover those logarithms from the form (A.5a), we refer to the behaviour of $\text{Li}_2(z)$ as $z \to -\infty$, which immediately
follows from the relation $\mathrm{Li}_2(z^{-1}) = -\pi^2/6 - \mathrm{Li}_2(z) - 1/2 \ln^2(-z)$, valid for $z \notin [0, +\infty[$. An integration by parts then shows that the results (A.5a) and (A.5b) differ by a mere constant in $r$ as expected. In our calculation of the NU metric, following the algorithm (2.22), all the logarithms of $r$ and associated constants, as well as all dilogarithm functions disappear in the end.

We treat now the case of hereditary terms of the type $I_m/r^k$, more general than the one we just studied. Those terms could be handled by means of successive integrations by parts on $r$, but it is more convenient to provide a formula that allows decreasing the power $k$ in their pre-factors, until we get back to the case $k=0$.

The idea consists in integrating by parts the hereditary term (A.1) over $x$. We introduce the anti-derivative of the Legendre function of the second kind that vanishes when $x \to 1^+$. Using the known identity $\frac{d}{dx}[Q_{m+1}(x) - Q_{m-1}(x)] = (2m+1)Q_m(x)$ for the Legendre function [88], we get (for $m \ge 1$)

$$ Q_m(x) \equiv \int_1^x dz Q_m(z) = \frac{1}{2m+1} (Q_{m+1}(x) - Q_{m-1}(x) + H_{m+1} - H_{m-1}). \quad (A.6) $$

Therefore, we find that an equivalent expression for $I_m$, which is obtained by integration by parts over $x$, reads

$$ I_m = \frac{r}{2m+1} \int_1^{+\infty} dx [Q_{m+1}(x) - Q_{m-1}(x)] F^{(1)}(t-rx) + \frac{H_{m+1} - H_{m-1}}{2m+1} F(u). \quad (A.7) $$

The remaining hereditary integral is now endowed with an explicit extra factor $r$, so that, by iterating the process, we end up considering only hereditary integrals without prefactor $r^{-k}$, i.e., the simpler case studied before.

In our calculation of the NU metric, we also need to differentiate $I_m$ with respect to $r$ at $u = \text{const}$. For completeness, we give the formula for the radial derivative of $I_m$, which relies on the fact that the Legendre function satisfies the recurrence formula $(2m+1)xQ_m(x) = (m+1)Q_{m+1}(x) + mQ_{m-1}(x)$ (for $m \ge 1$) [88]:

$$ \left. \frac{\partial I_m}{\partial r} \right|_{u=\text{const}} = - \frac{1}{2m+1} \int_1^{+\infty} dx \left[ (m+1)Q_{m+1}(x) - (2m+1)Q_m(x) + mQ_{m-1}(x) \right] F^{(1)}(t-rx). \quad (A.8) $$

Since the latter integral is a radial derivative, it can be directly integrated and, indeed, we see that the combination of Legendre functions in the integrand of (A.8) is locally integrable when $x \to 1^+$.

**Proof of the formula (A.3).** The coefficients $c_{mj}$ in this formula are those entering the decomposition of the polynomial of degree $m-1$

$$ \frac{P_m(x) - 1}{x - 1} = \sum_{j=0}^{m-1} c_{mj} P_j(x), \quad (A.9) $$
into a sum of ordinary Legendre polynomials, where we recall that $P_m(1) = 1$. The coefficients are readily computed as (for $j \le m$)

$$ 
\begin{aligned}
c_{mj} &= \frac{2j+1}{2} \int_{-1}^{1} dx \frac{P_m(x)-1}{x-1} P_j(x) = \frac{2j+1}{2} \lim_{z \to 1} \int_{-1}^{1} dx \frac{P_m(x)-1}{x-z} P_j(x) \\
&= (2j+1) \lim_{z \to 1} [-P_j(z)Q_m(z) + Q_j(z)] = (2j+1)(H_m - H_j).
\end{aligned} 
\quad (A.10) $$

Notice the side results

$$ \sum_{j=0}^{m-1} c_{mj} = \frac{m(m+1)}{2}, \qquad \sum_{j=0}^{m-1} c_{mj} H_j = \frac{m+1}{2} (mH_m - m + 1), \quad (A.11) $$

which may be derived by commuting the sum on $j$ with the one originating from the definition of harmonic numbers. Next, we know that the Legendre function reads

$$ Q_m(x) = \frac{1}{2} P_m(x) \ln \left( \frac{x+1}{x-1} \right) - W_{m-1}(x), \quad (A.12) $$

where $W_{m-1}(x)$ is a polynomial of degree $m-1$ which consists of the positive and zero powers of $x$ in the expansion of $\frac{1}{2}P_m(x)\ln(\frac{x+1}{x-1})$ in descending powers of $x$, i.e., when $x \to +\infty$ [88]. Notably, we have $W_{-1}(x) = 0$ and $W_{m-1}(1) = H_m$. Using (A.12) together with the decomposition (A.9) yields

$$ \frac{Q_m(x) - Q_0(x) + H_m}{x-1} = \sum_{j=0}^{m-1} c_{mj} Q_j(x) + R_{m-2}(x), \quad (A.13a) $$

$$ \text{where } R_{m-2}(x) = \sum_{j=0}^{m-1} c_{mj} W_{j-1}(x) - \frac{W_{m-1}(x) - W_{m-1}(1)}{x-1}. \quad (A.13b) $$

The point is that $R_{m-2}(x)$ is a polynomial, with degree $m-2$. However, the existence of non-vanishing coefficients in $R_{m-2}(x)$ is incompatible with taking the limit when $x \to +\infty$ on both sides of eq. (A.13a), since $Q_m(x)$ tends to zero like $1/x^{m+1}$ when $x \to +\infty$ (see e.g., eqs. (A3) in [89]). We conclude that this polynomial must be identically zero: $R_{m-2}(x) \equiv 0$, hence our formula (A.3) is proven. Moreover, we see that (A.3) is equivalent to

$$ \frac{W_{m-1}(x) - W_{m-1}(1)}{x-1} = \sum_{j=0}^{m-1} c_{mj} W_{j-1}(x). \quad (A.14) $$

## A.2 Far-zone expansion

In order to determine the far-zone behaviour of the hereditary integral $I_m$, we write it in terms of the variable $\tau = r(x-1)$ as

$$ I_m = \frac{1}{r} \int_0^{+\infty} d\tau Q_m \left(1 + \frac{\tau}{r}\right) F(u-\tau). \quad (A.15) $$

Since the integration range is actually a bound interval, we can substitute to $Q_m(x)$ its asymptotic expansion near $x \to 1^+$ for large enough $r$. This expansion is derived from the following suitable representation of the Legendre function:

$$ Q_m(x) = \sum_{j=0}^{m} \frac{(m+j)!}{(m-j)!j!^2} \left(\frac{x-1}{2}\right)^j \left[ \frac{1}{2} \ln\left(\frac{x+1}{x-1}\right) + H_j - H_m \right]. \quad (A.16) $$
Expanding $\ln(\frac{x+1}{2})$ when $x \to 1^+$, commuting the sums, and resorting to standard resum-
mation techniques including the identification of hypergeometric functions,<fnref n="16" /> we get

$$
\begin{equation}
\begin{split}
Q_m(x) \stackrel{x \to 1^+}{\sim} \frac{1}{2} \sum_{j=0}^{m} \frac{(m+j)!}{(m-j)! j!^2} \left[ 2H_j - H_{m+j} - H_{m-j} - \ln\left(\frac{x-1}{2}\right) \right] \left(\frac{x-1}{2}\right)^j \\
\qquad + \frac{1}{2} \left(\frac{x-1}{2}\right)^{m+1} \sum_{j=0}^{+\infty} (-)^j \frac{j!(2m+j+1)!}{(m+j+1)!^2} \left(\frac{x-1}{2}\right)^j.
\end{split}
\tag{A.17}
\end{equation}
$$

An alternative form of this expansion is given by eq. (4.6) in [90]. Having inserted the expression (A.17) in eq. (A.15), we permute the sums and the integrals, perform series of integration by parts to separate “instantaneous” terms (but which involve anti-derivatives of the function F) from logarithmic kernel hereditary integrals. We finally obtain the expansion, for $r \to +\infty$ at $u = \text{const}$,

$$
\begin{equation}
\begin{split}
I_m^{r \to +\infty} \approx {} & \sum_{j=0}^{m} \frac{(m+j)!}{(2r)^{j+1}(m-j)!j!} \left[ (H_j - H_{m+j} - H_{m-j}) \stackrel{(-j-1)}{F}(u) - \int_0^{+\infty} d\tau \ln\left(\frac{\tau}{2r}\right) \stackrel{(-j)}{F}(u-\tau) \right] \\
& + \sum_{j=m+1}^{+\infty} (-)^{j+m+1} \frac{(j-m-1)!(j+m)!}{(2r)^{j+1}j!} \stackrel{(-j-1)}{F}(u).
\end{split}
\tag{A.18}
\end{equation}
$$

For the expansion of $I_m$ in the near zone, i.e., $r \to 0$ with $u$ or $t$ fixed, we refer the reader to appendix A of [89].

The appearance of arbitrarily high order anti-derivatives reflects the non-locality of
$I_m$. Note also the presence of terms proportional to $\ln r/r^{j+1}$ with $j = 0, \dots, m$. Those
logarithms, however, have to disappear from the NU metric components by construction.
Moreover, for the quadrupole-quadrupole interaction $M_{ij} \times M_{ij}$, it turns out that all hered-
itary integrals $I_m$ cancel each other once they have been transformed with the help of the
power reduction formula (A.7) so as to bear the same common pre-factor $r^{-k_0}$ for some
chosen $k_0$. The $M_{ij} \times M_{ij}$ non-local terms of the NU metric are thus of pure memory type
and their far-zone expansion involves only a finite number of terms, which contrasts with
the situation in harmonic coordinates.

B Dressed Bondi aspects for $M_{ij} \times M_{ij}$ interactions

Including only mass monopole and quadrupoles we have

$$
C_{bc} = 2G e_{\langle b}^{i} e_{c\rangle}^{j} M_{ij}^{(2)} + \mathcal{O}(G^2), \qquad (B.1a)
$$

$$
D^c C_{bc} = -4G e_b^i n^j M_{ij}^{(2)} + \mathcal{O}(G^2), \qquad (B.1b)
$$

$$
D_a D^c C_{bc} = 4G \left( n_i n_j M_{ij}^{(2)} \gamma_{ab} - e_a^i e_b^j M_{ij}^{(2)} \right) + \mathcal{O}(G^2). \quad (B.1c)
$$

$$
_4F_3(-n, 1, 1, x; 2, y, x-y-n+1; 1) = \frac{(y-1)(x-y-n)}{(n+1)(x-1)} [-\psi(y-1) + \psi(x-y+1) + \psi(y+n) - \psi(x-y-n)],
$$

valid for $n \in \mathbb{N}$, where $_4F_3$ is a generalized hypergeometric function and $\psi$ is the digamma function.

<footnote n="16"><sup>16</sup>In this instance, one may resort to the identity</footnote>
Notice that eq. (B.1c) is symmetric in (ab), hence the identity $D_{[a}D^c C_{b]c} = \mathcal{O}(G^2)$. We then deduce from this<fnref n="17" />

$$D_a (C^2) = -8G^2 e_a^i M_{ij}^{(2)} (2M_{jk}^{(2)} - n_j n_l M_{lk}^{(2)}) n_k + \mathcal{O}(G^3), \quad (B.2a)$$

$$C_a^b D^c C_{bc} = -4G^2 e_a^i M_{ij}^{(2)} (2M_{jk}^{(2)} - n_j n_l M_{lk}^{(2)}) n_k + \mathcal{O}(G^3), \quad (B.2b)$$

$$D^b D_{[a} D^c C_{b]c} = 15G^2 e_a^l n^m \epsilon_{lmj} n^i \epsilon_{pq\langle i} \left( \frac{1}{2} M_{jp\underline{l}}^{(1)} M_{k\rangle q}^{(4)} - \frac{1}{10} M_{jp\underline{l}} M_{k\rangle q}^{(5)} \right) n_k + \mathcal{O}(G^3), \quad (B.2c)$$

while $D_a m$ is immediately computed from eq. (2.47a).

Resorting to a computer algebra software, we now include all quadratic-quadratic interactions to obtain

$$\begin{aligned} \mathcal{N}_i = & G \left[ 3M_i - 3uP_i + \hat{n}^m (6M_{im}^{(1)} - 6uM_{im}^{(2)}) \right]^T + G^2 \left\{ \hat{n}^m \left[ \frac{12}{5} \int_{-\infty}^u dv M_{n|m}^{(2)}(v) M_{i|n}^{(3)}(v) \right. \right. \\ & \left. \left. - \frac{264}{35} M_{n(m}^{(1)} M_{i)n}^{(3)} - \frac{12}{35} M_{n(m} M_{i)n}^{(4)} + u \left( \frac{264}{35} M_{n(m}^{(2)} M_{i)n}^{(3)} + \frac{276}{35} M_{n(m}^{(1)} M_{i)n}^{(4)} + \frac{12}{35} M_{n(m} M_{i)n}^{(5)} \right) \right] \right. \\ & \left. + \hat{n}^{mnj} \left[ -15M_{im}^{(2)} M_{nj}^{(2)} - 15M_{mn}^{(1)} M_{ij}^{(3)} - 6M_{im}^{(1)} M_{nj}^{(3)} - \frac{9}{2} M_{mn} M_{ij}^{(4)} - 6M_{im} M_{nj}^{(4)} \right. \right. \\ & \left. \left. + u \left( 51M_{m(n}^{(2)} M_{i)j}^{(3)} + \frac{39}{2} M_{mn}^{(1)} M_{ij}^{(4)} + 12M_{im}^{(1)} M_{nj}^{(4)} + \frac{9}{2} M_{mn} M_{ij}^{(5)} + 6M_{im} M_{nj}^{(5)} \right) \right] \right\}^T, \quad (B.3) \end{aligned}$$

$$\begin{aligned} \mathcal{E}_{ij} = & G(2M_{ij} - 2uM_{ij}^{(1)} + u^2 M_{ij}^{(2)})^{\text{TT}} + G^2 \left\{ -\frac{8}{3} M_{im}^{(1)} M_{jm}^{(2)} - \frac{52}{63} M_{im} M_{jm}^{(3)} + u \left( 2M_{im}^{(2)} M_{jm}^{(2)} \right. \right. \\ & \left. \left. + \frac{92}{21} M_{im}^{(1)} M_{jm}^{(3)} + \frac{22}{21} M_{im} M_{jm}^{(4)} \right) + u^2 \left( -\frac{74}{21} M_{im}^{(2)} M_{jm}^{(3)} - \frac{19}{7} M_{im}^{(1)} M_{jm}^{(4)} - \frac{11}{21} M_{im} M_{jm}^{(5)} \right) \right. \\ & \left. + \hat{n}^{mn} \left[ -3M_{mn}^{(1)} M_{ij}^{(2)} - 8M_{im}^{(1)} M_{jn}^{(2)} + \frac{1}{2} M_{ij}^{(1)} M_{mn}^{(2)} - \frac{4}{3} M_{mn} M_{ij}^{(3)} - \frac{22}{3} M_{im} M_{jn}^{(3)} \right. \right. \\ & \left. \left. - \frac{11}{6} M_{ij} M_{mn}^{(3)} + u(9M_{im}^{(2)} M_{jn}^{(2)} + \frac{3}{2} M_{ij}^{(2)} M_{mn}^{(2)} + 5M_{mn}^{(1)} M_{ij}^{(3)} + 14M_{im}^{(1)} M_{jn}^{(3)} + 2M_{ij}^{(1)} M_{mn}^{(3)} \right. \right. \\ & \left. \left. + \frac{3}{2} M_{mn} M_{ij}^{(4)} + 7M_{im} M_{jn}^{(4)} + 2M_{ij} M_{mn}^{(4)} \right) + u^2 \left( -\frac{17}{4} M_{mn}^{(2)} M_{ij}^{(3)} - 17M_{im}^{(2)} M_{jn}^{(3)} \right. \right. \\ & \left. \left. - \frac{17}{4} M_{ij}^{(2)} M_{mn}^{(3)} - \frac{13}{4} M_{mn}^{(1)} M_{ij}^{(4)} - \frac{21}{2} M_{im}^{(1)} M_{jn}^{(4)} - 2M_{ij}^{(1)} M_{mn}^{(4)} - \frac{3}{4} M_{mn} M_{ij}^{(5)} - \frac{7}{2} M_{im} M_{jn}^{(5)} \right. \right. \\ & \left. \left. - M_{ij} M_{mn}^{(5)} \right) \right] \right\}^{\text{TT}}, \quad (B.4) \end{aligned}$$

$$C_{ab} = e_{\langle a}^k e_{b\rangle}^l \left( 2U_{kl} - n_m n_p \epsilon_{mn(k} V_{l)p} + \frac{1}{6} n_p n_q U_{klpq} \right),$$

$$D^c C_{bc} = -2n^k e_b^l \left( 2U_{kl} + \frac{1}{4} n_p n_q U_{klpq} - \frac{5}{4} n_m \epsilon_{mnl} V_{knp} n_p \right),$$

$$D_{[a} D^c C_{b]c} = 5 \left( e_{[a}^k e_{b]}^l n^m + \frac{1}{2} e_a^m e_b^l n^k \right) \epsilon_{mnl} V_{knp} n_p,$$

$$D^b D_{[a} D^c C_{b]c} = 15 e_a^l n^m \epsilon_{mnl} V_{knp} n_k n_p.$$

<footnote n="17">¹⁷It is useful to observe that, from $C_{ab} = e_{\langle a}^i e_{b\rangle}^j H_{ij}^{TT}$ [see (2.45)], one has</footnote>
$$
\begin{equation}
\begin{aligned}
\mathcal{E}_{ij} = {}& G^2 \left\{ -2M_{im}M_{jm}^{(2)} + u\hat{n}^m(-2M_mM_{ij}^{(2)} - 4M_iM_{jm}^{(2)}) + u\left(\frac{13}{3}M_{im}^{(1)}M_{jm}^{(2)} + \frac{41}{9}M_{im}M_{jm}^{(3)}\right) \right. \\
& \left. \qquad + u^2\left(-2M_{im}^{(2)}M_{jm}^{(2)} - \frac{14}{3}M_{im}^{(1)}M_{jm}^{(3)} - \frac{7}{3}M_{im}M_{jm}^{(4)}\right) + u^3\left(\frac{34}{9}M_{im}^{(2)}M_{jm}^{(3)} + \frac{7}{3}M_{im}^{(1)}M_{jm}^{(4)}\right) \right. \\
& \left. \qquad + \frac{7}{9}M_{im}M_{jm}^{(5)} \right) + \hat{n}^{mn} \left[ 6M_{im}^{(1)}M_{jn}^{(1)} + M_{ij}^{(1)}M_{mn}^{(1)} - \frac{15}{4}M_{mn}M_{ij}^{(2)} - \frac{57}{2}M_{im}M_{jn}^{(2)} - \frac{9}{2}M_{ij}M_{mn}^{(2)} \right. \\
& \left. \qquad + u\left(\frac{9}{2}M_{mn}^{(1)}M_{ij}^{(2)} + 13M_{im}^{(1)}M_{jn}^{(2)} + \frac{21}{4}M_{ij}^{(1)}M_{mn}^{(2)} + \frac{17}{3}M_{mn}M_{ij}^{(3)} + \frac{74}{3}M_{im}M_{jn}^{(3)}\right) \right. \\
& \left. \qquad + \frac{77}{12}M_{ij}M_{mn}^{(3)} \right) + u^2\left(-\frac{27}{4}M_{im}^{(2)}M_{jn}^{(2)} - \frac{37}{8}M_{ij}^{(2)}M_{mn}^{(2)} - \frac{29}{4}M_{mn}^{(1)}M_{ij}^{(3)} - \frac{49}{2}M_{im}^{(1)}M_{jn}^{(3)}\right. \\
& \left. \qquad - 5M_{ij}^{(1)}M_{mn}^{(3)} - \frac{23}{8}M_{mn}M_{ij}^{(4)} - \frac{49}{4}M_{im}M_{jn}^{(4)} - \frac{13}{4}M_{ij}M_{mn}^{(4)}\right) + u^3\left(\frac{119}{24}M_{mn}^{(2)}M_{ij}^{(3)}\right. \\
& \left. \qquad + \frac{119}{6}M_{im}^{(2)}M_{jn}^{(3)} + \frac{119}{24}M_{ij}^{(2)}M_{mn}^{(3)} + \frac{27}{8}M_{mn}^{(1)}M_{ij}^{(4)} + \frac{49}{4}M_{im}^{(1)}M_{jn}^{(4)} + \frac{11}{4}M_{ij}^{(1)}M_{mn}^{(4)}\right. \\
& \left. \qquad + \frac{23}{24}M_{mn}M_{ij}^{(5)} + \frac{49}{12}M_{im}M_{jn}^{(5)} + \frac{13}{12}M_{ij}M_{mn}^{(5)}\right) \right]^{\text{TT}}.
\end{aligned}
\tag{B.5}
\end{equation}
$$

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

**References**

[1] L. Blanchet, Contribution à l'étude du rayonnement gravitationnel émis par un système isolé (in French), habilitation thesis, Université Pierre et Marie Curie, Paris VI, Paris, France (1990).

[2] D. Christodoulou, *Nonlinear nature of gravitation and gravitational wave experiments*, Phys. Rev. Lett. **67** (1991) 1486 [INSPIRE].

[3] A.G. Wiseman and C.M. Will, *Christodoulou's nonlinear gravitational wave memory: evaluation in the quadrupole approximation*, Phys. Rev. D **44** (1991) R2945 [INSPIRE].

[4] K.S. Thorne, *Gravitational-wave bursts with memory: the Christodoulou effect*, Phys. Rev. D **45** (1992) 520 [INSPIRE].

[5] L. Blanchet and T. Damour, *Hereditary effects in gravitational radiation*, Phys. Rev. D **46** (1992) 4304 [INSPIRE].

[6] H. Bondi, M.G.J. van der Burg and A.W.K. Metzner, *Gravitational waves in general relativity*. 7. Waves from axisymmetric isolated systems, Proc. Roy. Soc. Lond. A **269** (1962) 21 [INSPIRE].

[7] R.K. Sachs, *Gravitational waves in general relativity*. 8. Waves in asymptotically flat space-times, Proc. Roy. Soc. Lond. A **270** (1962) 103 [INSPIRE].

[8] E.T. Newman and R. Penrose, *Note on the Bondi-Metzner-Sachs group*, J. Math. Phys. **7** (1966) 863 [INSPIRE].

[9] A. Strominger and A. Zhiboedov, *Gravitational memory, BMS supertranslations and soft theorems*, JHEP **01** (2016) 086 [arXiv:1411.5745] [INSPIRE].
[10] L. Blanchet, G. Faye, B.R. Iyer and S. Sinha, *The third post-Newtonian gravitational wave polarisations and associated spherical harmonic modes for inspiralling compact binaries in quasi-circular orbits*, Class. Quant. Grav. **25** (2008) 165003 [Erratum ibid. **29** (2012) 239501] [arXiv:0802.1249 [INSPIRE]].

[11] L. Blanchet and T. Damour, *Radiative gravitational fields in general relativity 1*. General structure of the field outside the source, Phil. Trans. Roy. Soc. Lond. A **320** (1986) 379 [INSPIRE].

[12] L. Blanchet, *Radiative gravitational fields in general relativity 2*. Asymptotic behaviour at future null infinity, Proc. Roy. Soc. Lond. A **409** (1987) 383 [INSPIRE].

[13] L. Blanchet and T. Damour, *Tail transported temporal correlations in the dynamics of a gravitating system*, Phys. Rev. D **37** (1988) 1410 [INSPIRE].

[14] K. Mitman et al., *Computation of displacement and spin gravitational memory in numerical relativity*, Phys. Rev. D **102** (2020) 104007 [arXiv:2007.11562] [INSPIRE].

[15] K. Mitman et al., *Fixing the BMS frame of numerical relativity waveforms*, Phys. Rev. D **104** (2021) 024051 [arXiv:2105.02300] [INSPIRE].

[16] M. Favata, *Nonlinear gravitational-wave memory from binary black hole mergers*, Astrophys. J. Lett. **696** (2009) L159 [arXiv:0902.3660] [INSPIRE].

[17] P.D. Lasky et al., *Detecting gravitational-wave memory with LIGO: implications of GW150914*, Phys. Rev. Lett. **117** (2016) 061102 [arXiv:1605.01415] [INSPIRE].

[18] L.O. McNeill, E. Thrane and P.D. Lasky, *Detecting gravitational wave memory without parent signals*, Phys. Rev. Lett. **118** (2017) 181103 [arXiv:1702.01759] [INSPIRE].

[19] J.B. Wang et al., *Searching for gravitational wave memory bursts with the Parkes Pulsar Timing Array*, Mon. Not. Roy. Astron. Soc. **446** (2015) 1657 [arXiv:1410.3323] [INSPIRE].

[20] NANOGRAV collaboration, *NANOGlav constraints on gravitational wave bursts with memory*, Astrophys. J. **810** (2015) 150 [arXiv:1501.05343] [INSPIRE].

[21] Y.B. Zel'dovich and A.G. Polnarev, *Radiation of gravitational waves by a cluster of superdense stars*, Sov. Astron. **18** (1974) 17 [INSPIRE].

[22] V.B. Braginsky and L.P. Grishchuk, *Kinematic resonance and memory effect in free mass gravitational antennas*, Sov. Phys. JETP **62** (1985) 427 [INSPIRE].

[23] M. Turner, *Gravitational radiation from point-masses in unbound orbits — Newtonian results*, Astrophys. J. **216** (1977) 610.

[24] R. Epstein, *The generation of gravitational radiation by escaping supernova neutrinos*, Astrophys. J. **223** (1978) 1037 [INSPIRE].

[25] L. Bieri, P.N. Chen and S.-T. Yau, *Null asymptotics of solutions of the Einstein-Maxwell equations in general relativity and gravitational radiation*, Adv. Theor. Math. Phys. **15** (2011) 1085 [arXiv:1011.2267] [INSPIRE].

[26] L. Bieri and D. Garfinkle, *Perturbative and gauge invariant treatment of gravitational wave memory*, Phys. Rev. D **89** (2014) 084039 [arXiv:1312.6871] [INSPIRE].

[27] S. Pasterski, A. Strominger and A. Zhiboedov, *New gravitational memories*, JHEP **12** (2016) 053 [arXiv:1502.06120] [INSPIRE].

[28] D.A. Nichols, *Center-of-mass angular momentum and memory effect in asymptotically flat spacetimes*, Phys. Rev. D **98** (2018) 064032 [arXiv:1807.08767] [INSPIRE].
[29] É.É. Flanagan, A.M. Grant, A.I. Harte and D.A. Nichols, *Persistent gravitational wave observables: general framework*, Phys. Rev. D 99 (2019) 084044 [arXiv:1901.00021] [INSPIRE].

[30] E. Himwich, Z. Mirzaiyan and S. Pasterski, *A note on the subleading soft graviton*, JHEP 04 (2021) 172 [arXiv:1902.01840] [INSPIRE].

[31] A.M. Grant and D.A. Nichols, *Persistent gravitational wave observables: curve deviation in asymptotically flat spacetimes*, Phys. Rev. D 105 (2022) 024056 [Erratum ibid. 107 (2023) 109902] [arXiv:2109.03832] [INSPIRE].

[32] A. Seraj, *Gravitational breathing memory and dual symmetries*, JHEP 05 (2021) 283 [arXiv:2103.12185] [INSPIRE].

[33] A. Seraj and B. Oblak, *Gyroscopic gravitational memory*, arXiv:2112.04535 [INSPIRE].

[34] A. Seraj and B. Oblak, *Precession caused by gravitational waves*, Phys. Rev. Lett. 129 (2022) 061101 [arXiv:2203.16216] [INSPIRE].

[35] M. Godazgar, G. Macaulay, G. Long and A. Seraj, *Gravitational memory effects and higher derivative actions*, JHEP 09 (2022) 150 [arXiv:2206.12339] [INSPIRE].

[36] G. Barnich and C. Troessaert, *Symmetries of asymptotically flat 4 dimensional spacetimes at null infinity revisited*, Phys. Rev. Lett. 105 (2010) 111103 [arXiv:0909.2617] [INSPIRE].

[37] G. Barnich and C. Troessaert, *Aspects of the BMS/CFT correspondence*, JHEP 05 (2010) 062 [arXiv:1001.1541] [INSPIRE].

[38] M. Campiglia and A. Laddha, *New symmetries for the gravitational S-matrix*, JHEP 04 (2015) 076 [arXiv:1502.02318] [INSPIRE].

[39] G. Compère, R. Oliveri and A. Seraj, *Gravitational multipole moments from Noether charges*, JHEP 05 (2018) 054 [arXiv:1711.08806] [INSPIRE].

[40] H. Godazgar, M. Godazgar and C.N. Pope, *New dual gravitational charges*, Phys. Rev. D 99 (2019) 024013 [arXiv:1812.01641] [INSPIRE].

[41] H. Godazgar, M. Godazgar and C.N. Pope, *Tower of subleading dual BMS charges*, JHEP 03 (2019) 057 [arXiv:1812.06935] [INSPIRE].

[42] G. Compère, *Infinite towers of supertranslation and superrotation memories*, Phys. Rev. Lett. 123 (2019) 021101 [arXiv:1904.00280] [INSPIRE].

[43] A. Strominger, $w_{1+\infty}$ and the celestial sphere, arXiv:2105.14346 [INSPIRE].

[44] D. Kapec, P. Mitra, A.-M. Raclariu and A. Strominger, *2D stress tensor for 4D gravity*, Phys. Rev. Lett. 119 (2017) 121601 [arXiv:1609.00282] [INSPIRE].

[45] H. Godazgar, M. Godazgar and C.N. Pope, *Dual gravitational charges and soft theorems*, JHEP 10 (2019) 123 [arXiv:1908.01164] [INSPIRE].

[46] L. Freidel, D. Pranzetti and A.-M. Raclariu, *Higher spin dynamics in gravity and $w_{1+\infty}$ celestial symmetries*, Phys. Rev. D 106 (2022) 086013 [arXiv:2112.15573] [INSPIRE].

[47] A. Papapetrou, *Coordonnées radiatives cartésiennes* (in French), Ann. Inst. Henri Poincaré A XI (1969) 251.

[48] J. Madore, *Gravitational radiation from a bounded source*. I, Ann. Inst. Henri Poincaré 12 (1970) 285.
[49] J. Madore, *Gravitational radiation from a bounded source*. II, *Ann. Inst. Henri Poincaré* **12** (1970) 365.

[50] E.T. Newman and T.W.J. Unti, *A class of null flat-space coordinate systems*, *J. Math. Phys.* **4** (1963) 1467 [INSPIRE].

[51] G. Barnich and P.-H. Lambert, *A note on the Newman-Unti group and the BMS charge algebra in terms of Newman-Penrose coefficients*, *Adv. Math. Phys.* **2012** (2012) 197385 [arXiv:1102.0589] [INSPIRE].

[52] R. Sachs, *Asymptotic symmetries in gravitational theory*, *Phys. Rev.* **128** (1962) 2851 [INSPIRE].

[53] M.G.J. van der Burg, *Gravitational waves in general relativity*. IX. *Conserved quantities*, *Proc. Roy. Soc. Lond. A* **294** (1966) 112.

[54] L.A. Tamburino and J.H. Winicour, *Gravitational fields in finite and conformal Bondi frames*, *Phys. Rev.* **150** (1966) 1039 [INSPIRE].

[55] J. Winicour, *Logarithmic asymptotic flatness*, *Found. Phys.* **15** (1985) 605.

[56] L. Freidel, D. Pranzetti and A.-M. Raclariu, *Sub-subleading soft graviton theorem from asymptotic Einstein's equations*, *JHEP* **05** (2022) 186 [arXiv:2111.15607] [INSPIRE].

[57] G. Barnich and C. Troessaert, *BMS charge algebra*, *JHEP* **12** (2011) 105 [arXiv:1106.0213] [INSPIRE].

[58] É.É. Flanagan and D.A. Nichols, *Conserved charges of the extended Bondi-Metzner-Sachs algebra*, *Phys. Rev. D* **95** (2017) 044002 [arXiv:1510.03386] [INSPIRE].

[59] G. Barnich, P. Mao and R. Ruzziconi, *BMS current algebra in the context of the Newman-Penrose formalism*, *Class. Quant. Grav.* **37** (2020) 095010 [arXiv:1910.14588] [INSPIRE].

[60] G. Compère, R. Oliveri and A. Seraj, *The Poincaré and BMS flux-balance laws with application to binary systems*, *JHEP* **10** (2020) 116 [arXiv:1912.03164] [INSPIRE].

[61] L. Blanchet et al., *Multipole expansion of gravitational waves: from harmonic to Bondi coordinates*, *JHEP* **02** (2021) 029 [arXiv:2011.10000] [INSPIRE].

[62] L. Blanchet, B.R. Iyer, C.M. Will and A.G. Wiseman, *Gravitational wave forms from inspiralling compact binaries to second post-Newtonian order*, *Class. Quant. Grav.* **13** (1996) 575 [gr-qc/9602024] [INSPIRE].

[63] K.G. Arun, L. Blanchet, B.R. Iyer and M.S.S. Qusailah, *The 2.5PN gravitational wave polarisations from inspiralling compact binaries in circular orbits*, *Class. Quant. Grav.* **21** (2004) 3771 [Erratum ibid. **22** (2005) 3115] [gr-qc/0404085] [INSPIRE].

[64] L.E. Kidder, L. Blanchet and B.R. Iyer, *Radiation reaction in the 2.5PN waveform from inspiralling binaries in circular orbits*, *Class. Quant. Grav.* **24** (2007) 5307 [arXiv:0706.0726] [INSPIRE].

[65] D.A. Nichols, *Spin memory effect for compact binaries in the post-Newtonian approximation*, *Phys. Rev. D* **95** (2017) 084048 [arXiv:1702.03300] [INSPIRE].

[66] L. Blanchet and G. Faye, *Flux-balance equations for linear momentum and center-of-mass position of self-gravitating post-Newtonian systems*, *Class. Quant. Grav.* **36** (2019) 085003 [arXiv:1811.08966] [INSPIRE].
[67] E.T. Newman and R. Penrose, 10 exact gravitationally-conserved quantities, Phys. Rev. Lett. 15 (1965) 231 [INSPIRE].

[68] E.T. Newman and R. Penrose, New conservation laws for zero rest-mass fields in asymptotically flat space-time, Proc. Roy. Soc. Lond. A 305 (1968) 175 [INSPIRE].

[69] G. Compère, A. Fiorucci and R. Ruzziconi, Superboost transitions, refraction memory and super-Lorentz charge algebra, JHEP 11 (2018) 200 [Erratum ibid. 04 (2020) 172] [arXiv:1810.00377] [INSPIRE].

[70] R. Sachs and P.G. Bergmann, Structure of particles in linearized gravitational theory, Phys. Rev. 112 (1958) 674 [INSPIRE].

[71] F.A.E. Pirani, Introduction to gravitational radiation theory, volume 1 of Brandeis summer institute in theoretical physics, Prentice-Hall, Englewood Cliffs, NJ, U.S.A. (1964), p. 249.

[72] K.S. Thorne, Multipole expansions of gravitational radiation, Rev. Mod. Phys. 52 (1980) 299 [INSPIRE].

[73] L. Blanchet and T. Damour, Radiative gravitational fields in general relativity I. General structure of the field outside the source, Phil. Trans. Roy. Soc. Lond. A 320 (1986) 379 [INSPIRE].

[74] R. Epstein and R.V. Wagoner, Post-Newtonian generation of gravitational waves, Astrophys. J. 197 (1975) 717.

[75] C.N. Kozameh and G.D. Quiroga, Center of mass and spin for isolated sources of gravitational radiation, Phys. Rev. D 93 (2016) 064050 [arXiv:1311.5854] [INSPIRE].

[76] C.N. Kozameh, J.I. Nieva and G.D. Quiroga, Spin and center of mass comparison between the post-Newtonian approach and the asymptotic formulation, Phys. Rev. D 98 (2018) 064005 [arXiv:1711.11375] [INSPIRE].

[77] L. Blanchet, Quadrupole-quadrupole gravitational waves, Class. Quant. Grav. 15 (1998) 89 [gr-qc/9710037] [INSPIRE].

[78] V.B. Braginsky and K.S. Thorne, Gravitational-wave bursts with memory and experimental prospects, Nature 327 (1987) 123.

[79] P.N. Payne, Smarr's zero frequency limit calculation, Phys. Rev. D 28 (1983) 1894 [INSPIRE].

[80] M. Ludvigsen, Geodesic deviation at null infinity and the physical effects of very long wave gravitational radiation, Gen. Rel. Grav. 21 (1989) 1205 [INSPIRE].

[81] J.M. Martín-García et al., xAct: efficient tensor computer algebra for Mathematica, http://www.xact.es/, GPL (2002)–(2012).

[82] L. Blanchet, Gravitational radiation from post-Newtonian sources and inspiralling compact binaries, Living Rev. Rel. 17 (2014) 2 [arXiv:1310.1528] [INSPIRE].

[83] G. Compère, R. Oliveri and A. Seraj, Metric reconstruction from celestial multipoles, JHEP 11 (2022) 001 [arXiv:2206.12597] [INSPIRE].

[84] A.R. Exton, E.T. Newman and R. Penrose, Conserved quantities in the Einstein-Maxwell theory, J. Math. Phys. 10 (1969) 1566 [INSPIRE].

[85] L. Freidel, R. Oliveri, D. Pranzetti and S. Speziale, The Weyl BMS group and Einstein's equations, JHEP 07 (2021) 170 [arXiv:2104.05793] [INSPIRE].
[86] L. Freidel and D. Pranzetti, Gravity from symmetry: duality and impulsive waves, JHEP 04 (2022) 125 [arXiv:2109.06342] [INSPIRE].

[87] H. Godazgar, M. Godazgar and C.N. Pope, Subleading BMS charges and fake news near null infinity, JHEP 01 (2019) 143 [arXiv:1809.09076] [INSPIRE].

[88] E. Whittaker and G. Watson, A course of modern analysis, Cambridge University Press, Cambridge, U.K. (1990).

[89] L. Blanchet, G. Faye and F. Larrouturou, The quadrupole moment of compact binaries to the fourth post-Newtonian order: from source to canonical moment, Class. Quant. Grav. 39 (2022) 195003 [arXiv:2204.11293] [INSPIRE].

[90] D. Trestini, F. Larrouturou and L. Blanchet, The quadrupole moment of compact binaries to the fourth post-Newtonian order: relating the harmonic and radiative metrics, Class. Quant. Grav. 40 (2023) 055006 [arXiv:2209.02719] [INSPIRE].