Multipole expansion of gravitational waves: memory effects and Bondi aspects

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.


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][2][3][4][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 Bondivan 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][12][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][17][18] or pulsar timing arrays [19,20]. 1 .
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][63][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 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][59][60]69] and provide a corrected definition of the dressed n = 3 Bondi aspect alternative to [31], which was subsequently revised in [70].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 twodimensional sphere, while lower case Latin indices from i to z will refer to threedimensional Cartesian indices.The Minkowski metric is η µν = diag(−1, +1, +1, +1).
Furthermore, L = i 1 i 2 • • • i represents a multi-index made of spatial indices.We use short-hands for: the multi-derivative operator 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 In retarded spherical coordinates (u, r, θ a ) with u = t−r, we have k µ ∂ µ = ∂ r u .We commonly employ the natural basis on the unit 2-sphere embedded in R 3 , e a = ∂ ∂θ a , with components e i a = ∂n i /∂θ a .Given the metric on the unit sphere γ ab = diag(1, sin 2 θ), we have: n i e i a = 0, ∂ i θ a = r −1 γ ab e i b , γ ab = δ ij e i a e j b and γ ab e i a e j b =⊥ ij , where ⊥ ij = δ ij −n i n j is the projector onto the sphere tangent bundle.We also use the notation e i a e j b = e i (a e j b) − 1 2 γ ab ⊥ ij for the trace-free product of basis vectors, where round brackets denote symmetrization e i (a e j b) = 1 2 (e i a e j b + e i b e j a ).On the other hand, square brackets denote anti-symmetrization, e.g., e i [a e j b] = 1 2 (e i a e j b − e i b e j a ).The transverse-trace-free (TT) projection operator is defined as Introducing the covariant derivative D a compatible with the sphere metric, D a γ bc = 0, we can write . An arbitrary rank 2 STF tensor over the sphere may be decomposed as where ab ≡ e i a e j b n k kij is the unit sphere volume form and the modal decomposition of T ± = 2 T ± L n L , with the T ± L being STF over the multi-index L, contains only 2 harmonics without loss of generality.Note the properties D a e i b = 0, D a (e i a e j b ) = −2n (i e j) b , ∆T ± = − ( + 1)n L T ± L , and (∆ + 2)(e i a e j b T + ij ) = 0. We define the unitnormalized integral over the sphere as S 1 = 1.
Given a general manifold, harmonic/de Donder coordinates are specified by using a tilde: xµ = ( t, x) or ( t, r, θa ).The metric components are gµν (x).Asymptotically flat spacetimes admit, as a background structure, the Minkowskian outgoing vector kµ = (1, ñi ), the basis on the sphere ẽi a = ∂ ñi /∂ θa , etc.We define the retarded time ũ in harmonic coordinates as ũ = t − r, so that kµ = − ∂µ ũ.NU coordinates are denoted as x µ = (u, r, θ a ) with θ a = (θ, ϕ).The metric components in those coordinates are g µν (x), with all other notations, such as the natural basis on the sphere e i a and the metric γ ab , as previously.
2 Quadratic memory from the MPM formalism

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 ), with each PM coefficient generated in the form of a multipole expansion starting from the linearized vacuum metric in harmonic coordinates [71][72][73][74] (setting c = 1) ) 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 ũ = t − r.For the following, it is convenient to write the dominant 1/r piece of h µν 1 (when r → +∞ with fixed ũ) as 2 where M is the constant Arnowitt-Deser-Misner (ADM) energy, S i is the constant ADM angular momentum and P i ≡ M (1) i 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 2 This tensorial quantity is presented in the form Note that kν z µν 1 = 0 as the consequence of the harmonic gauge condition ∂ν h µν 1 = 0, where we have introduced the Minkowski outgoing null vector kµ = (1, ñi ).At quadratic order, h µν 2 obeys the flat wave equation ˜ h µν 2 = N µν 2 (together with ∂ν h µν 2 = 0), where the source term N µν 2 is quadratic in h 1 as well as its first and second spacetime 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 antiderivatives 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/r 2 piece in N µν 2 , and that this piece takes the form 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 kµ = (1, ñi ).The quantity Π is directly given by a quadratic product of the multipole expansion (2.3): The GW energy flux emitted in the solid angle d Ω around the direction ñ is 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/r components of h µν 2 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/r 2 are easily computed from (2.29-30) of that paper.We can write . (2.8) The tail term comes directly from the first term in eq.(2.5); at the dominant order 1/r, it involves a logarithmic kernel, (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/r: 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 r components of the radiative metric to quadratic order, they belong, like the NU coordinates, to the large class of radiative coordinate systems [47][48][49].Although these coordinates are not harmonic either (we may call them modified harmonic coordinates), we will still use tildes to denote them, xµ = (ũ, r, ñ).The leading term in 1/r is just a time anti-derivative (sometimes called "semihereditary").It is a combination of the memory and also the secular mass loss ∝ Π 0 in the 00 component of K µν .To define it, we decompose the quantity Π into STF multipolar pieces as where Then, the components of K µν are explicitly given by3 Notice the useful properties, posing K ≡ η µν K µν : kν To order 1/r, the memory contribution to the transverse-traceless (TT) asymptotic waveform reads since all terms in K kl proportional to either ñk , ñl or δ kl are killed by the TT projection.The global minus sign arises from the relation g TT ij = −h TT ij , where h ij is the perturbation of the gothic metric (following the convention of [61]).
As for the subleading 1/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 ũ −T ).In this situation, all anti-derivatives are well defined.
From eq. (2.7), the angular average Π 0 = d Ω 4π Π is proportional to the total energy flux carried by gravitational waves: Similarly, Π i = 3 d Ω 4π ñi Π is proportional to the total linear momentum flux.We also introduced the notation Σ i for the angular momentum flux and Λ i for the center-of-mass flux: (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 [73,75], except for the one associated with the position of the center-of-mass, obtained only recently from the asymptotic properties of the radiation field [76,77], using the traditional PN as well as the Bondi approaches [28,60,65,66].The complete multipole expansions of the fluxes to quadratic order (∝ G 2 ) are given by4 S iL . (2.17d) The center-of-mass flux formula Λ 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 [76,77].In [66], the following simpler alternative center-of-mass flux was derived: It is associated with a different center-of-mass vector Gi .Since Λi − Λ i = ∂ ũ∆G i is a total ũ derivative, which vanishes on average for periodic systems, the two formulations . These definitions correspond to α = 1 = β in the parametrization given in Eq. (3.18) of [60].The Kerr black hole has E = M and J z = M a with the standard orientation trθφ = 1.
of the flux-balance law are physically equivalent.One can relate the corresponding center-of-mass vectors by Gi = G i + ∆G i and the flux-formulae differ due to this shift at fixed retarded time ũ.A disadvantage of the shifted Gi , 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 Λ i henceforth.
Notice that one only needs the leading term (2.3) in the 1/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/r 2 (see [66] for details).On the other hand, the Poincaré fluxes become exact expressions when rewritten in terms of radiative multipole moments

Given a metric perturbation
) in an arbitrary (not necessarily harmonic) coordinate system {ũ, r, θa }, we constructed perturbatively in [61] the gauge transformation {ũ, r, θa } −→ {u, r, θ a } required to reach the Newman-Unti (NU) gauge.Denoting the coordinate transformation as 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 kµ is the Minkowski null outgoing vector defined in the original coordinates) and next, the quadratic order equations (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 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 involves the mass-loss and memory terms.Moreover, the differential operator kµ ∂µ = ∂r | ũ is simply solved by integration over r at constant ũ.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 We now implement our algorithm by solving the system (2.22), focusing on the memory and mass-loss terms.Reminding the property kν K µν = 0, we see from eq. (2.10) that the first equation (2.22a) reduces to kµ ∂µ U 2 = O(r −2 ) at leading order, which is directly solved to U 2 = U 0 2 (ũ, ñ) + O(r −1 ).After imposing asymptotic flatness and setting the BMS transformation at second order to zero, we get U 0 2 (ũ, ñ) = 0. Therefore, the only contributions to U 2 come from the 1/r 2 term in h µν 2 as given by eq.(2.10).This term is easily integrated and we obtain Σ q . (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, θ a } using the inverse of eqs.(2.19).In particular, this entails re-expanding the 2-sphere metric as where Γ a bc is the Christoffel symbol associated with the covariant derivative on the sphere.This brings a memory correction coming from Θ a 2 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 The above metric is entirely expressed in terms of the NU coordinates {u, r, θ 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][58][59]61] ) ) together with g rr = g ra = 0 and g ur = −1.Here, m is the Bondi mass aspect, C ab is the shear,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 There is no memory contribution in the combination Ṅ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 [78] 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 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 ×M ij in terms of the canonical mass quadrupole M ij .The corresponding quadrupole-quadrupole metric in harmonic coordinates was computed in [78].From eqs. (4.5) of [78], for quadratic products of two mass quadrupoles, the only non-zero multipolar coefficients are Π 0 = 4 5 (3) (3) M kl . (2.30) This yields M j k − 4 5r (2) M ik (v) , (2.31b) (3) M kl (v) . (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 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 i a e j b H TT ij [see e.g., [61] and eq.(2.14) above].Using eq.(2.7) and eq.(2.11), we obtain where the infinite multipole series therein can be summed up in closed form as (see below for the proof) 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, 79]: where the factor 1/(1 − n • 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 [80,81]).
Proof of the formula (2.34).We consider the following TT projection with respect to the unit vector n = (n i ): where we recall that the TT projection reads We expand the denominator in eq.(2.36) as a power series in n • n (which is convergent as soon as n • n < 1), hence (2.37) The multi-index L − 2 contains − 2 indices, namely a 1 • • • a −2 .Using eq.(A.21a) in [74], we transform the ordinary product of unit vectors n ijL−2 into a sum of STF products ) The operation over indices {} is defined as the un-normalized sum over the smallest set of permutations of i 1 • • • i which makes the object symmetrical in with, say, i = a and j = a −1 ) of which 2k are displayed onto the product of k Kronecker symbols denoted δ 2K .As an example we have The point is that, among all the terms composing δ {2K n L−2K} , we can discard all those which contain either δ ij , δ iap or δ jaq , since such terms will be cancelled by the TT projection.This is obvious for δ ij ; in the two other cases, this results from the fact that, after multiplication by n L−2 in eq.(2.38a), δ iap or δ jaq 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: (2.39) The next step is to notice from eq. (A.19) in [74] that the number of terms composing the object δ {2K n ijL−2−2K} is When contracted with n L−2 , all these terms will merge into a single one for each values of and k.Thus, we have We change into + 2k and rewrite the previous expression as introducing the coefficient S which is given as an infinite series over all integer values of k.However, we find, using the expression of the coefficients α +2k k deduced from eq. (2.38b), that this series can actually be re-summed in closed analytic form: Hence we obtain the simple result Recalling that e i a e j b ⊥ TT ijkl = e k a e l b , we have therefore proved the formula (2.34).
We now obtain the complete NU metric corresponding to the quadrupole-quadrupole interaction M ij × M ij and tails M × 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 × M ij interaction, we start from the metric in harmonic coordinates as given in [78], instead of the modified harmonic coordinates described in Sec 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. 7ur final results for the NU metric (see Sec. 3.2 for the link with the Bondi metric) read ) ) . (2.44c) Here, H TT ij is the radiation field, which is a free data relating to the Bondi shear as C ab = e i a e j b H TT ij .It can be decomposed into 2 mass/electric and current/magnetic radiative multipoles U L , V L as The nonzero radiative multipole moments are given by (displaying only the terms that are linear in M ij or that correspond to the ) where time derivatives of moments are denoted, e.g., M (q) ij ≡ (q) M ij , 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) ) jk M (3) ik M (3) jk M (3) (2.47b) The trace of the angular part of the metric (2.44c) is ik M (2) 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 (2) (2) kl ik M (1) and where the superscript T refers to transverse projection, i.e., X T i ≡ ⊥ ij X j .For the subleading terms in the angular part of the metric, we find jp + npq 6M (1) = G 2 npq 173 20 We recall that the TT projection is X TT ij ≡ ⊥ ijmn TT X mn , with ⊥ ijmn TT =⊥ m(i ⊥ j)n − 1 2 ⊥ ij ⊥ 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 × 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 (−T ), where −T is the finite instant before which the multipole moment is constant, and a dynamical part M ij (u) − M ij (−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 Furthermore, we provide here the regular expansions of the integrals appearing in the ua and ab components: 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 3. Explicitly, the M × M ij contributions to m, N i , H T T ij and (n) E ij are given by ) (2) (2.53f)

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,84].

Newman-Penrose charges
The Newman-Penrose (NP) charges [53,67,68,85] 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 = 0, 1 multipole moments: M, M i ≡ 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 ), 0, S L (−T ), 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 in terms of a null tetrad, consisting of the incoming and outgoing real null vectors and the complex null vector m and its complex conjugate m, with In the above equations, β, F, U a , ζ a and ω are functionals of the Bondi data that are fixed by the normalization property of the tetrad (see e.g., [33,34,59,86]).The pair (ζ a , ζa ) forms a null dyad on the sphere normalized as γ ab ζ a ζb = 1.The bulk extension of m is fixed by requiring that it is parallelly transported along the outgoing null vector, i.e. ν ∇ ν m µ = 0, as in [68].One may choose the dyad to be adapted to the stereographic coordinates (z, z) on the sphere, in which the metric reads ds 2 = 4 (1+z z) 2 dzdz, so that ζ a ∂ a = 1+z z 2 ∂ z .Residual transformations of the tetrad only include global Type III rotations that covariantly transform ζ a . 9Computing the asymptotic form of this Weyl scalar, we find that Ψ 0 = Ψ 0 0 r −5 + Ψ 1 0 r −6 + O(r −7 ), where the leading and subleading coefficients are The (complex-valued) Newman-Penrose charges are defined, for m = −2, −1, 0, 1, 2, as where s Y lm denote the spin weighted spherical harmonics of spin weight s (the bar indicates the complex conjugate).They are given for 0 s l in terms of the Geroch-Held-Penrose operator ð10 by and −i E ab , where, for any STF tensor X ab , we pose X ab = ac X c b .Now, we make a change of basis from spherical harmonics Y 2m to STF harmonics nij , and we adjust the normalization to define STF NP charges Q ij as . (3.8) Given that D a n i = e a i and D a D b n i = −γ ab n i , it can be checked that D a D b nij = 2e a i e b j .Therefore, we arrive at −i E ab − where we have resorted to the Cartesian presentation of the Bondi data ≡ e a i e b j E ab ≡ e a i e b j E ab = n m min E nj .
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) E ij contain monopole terms.Thus, there is no mass quadrupole-quadrupole contribution to the NP charges.The mass monopolequadrupole contribution to the NP charges is a direct consequence of eqs.(2.53) obtained from eq. (4.13b) of [61].Expanding ⊥ TT ijkl in the expression of (3) E ij in eq.(2.53e) one finds 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.

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 App.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 θ a unchanged.In Bondi gauge, the spherical metric takes the form [31] . (3.12)The relation between the two radii follows from We shall parametrize the NU metric as where W = O(G 2 ) and C NU ab = O(G).Then, eq.(3.13) yields where W and C NU ab C ab NU can be evaluated either at radius r NU or r B , neglecting O(G 4 ) corrections.At large radius, it can be shown, using g B ur (u, r B , θ a ) = −∂r NU /∂r B and g B ab (u, r B , θ a ) = g NU ab (u, r NU , θ a ), combined with the 1/r B expansion of g B ur , that and therefore which reproduces eq.(61b) of [61].Including all orders in the radius, we may pose 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 (k) 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 where f = T (θ a ) + u 2 D a Y a results from a BMS transformation.In the absence of BMS transformation, we then find The Bondi mass aspect can be read in Bondi gauge from g uu = −1 + 2m/r + O(r −2 ).The Bondi angular momentum aspect N a can be read in Bondi gauge from 12 Since the Bondi radius r B differs from the Newman-Unti radius r NU by terms of order are of order G, we simply have 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.

Flux-balance laws of the dressed n 3 Bondi aspects
In the Bondi expansion, Einstein's equations lead to time evolution equations for all n 0 Bondi aspects.The Bondi aspects can be further "dressed" by adding suitable nonlinear 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 masscurrent) duality covariant mass [40,41,45] 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 ) 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 [86,87] is equal to N a here, and N a as defined in [37] is equal to ∂ u E ab where we use the compact notation C 2 ≡ C ab C ab , (CN ) ≡ C ab N ab , N 2 ≡ N ab N ab , and where D 0 ≡ − 1 4 (∆ + 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,88]).
The n = 0 Bondi aspect m is just the Bondi mass aspect that permits defining the supermomenta i.e., the canonical charges associated with supertranslations T = T L nL .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 = 0, the energy is defined as , where E 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 as13 They obey the flux-balance laws The n = 1 dressed aspect allows defining the super-Lorentz conserved charges The dressed n = 3 Bondi aspect, which is conserved in non-radiative regions, was introduced in eq.(4.43) of [31]. 14It can be written in terms of the duality-covariant quantity m ab as it should come from gravitational electric-magnetic duality [84].The (corrected) dressed n = 3 Bondi aspect reads as It obeys the flux-balance law where the flux (3) F ab can be written, using the flux-balance laws (3.24) and the identity C c(a N c b) = 1 2 γ ab (CN ), as As expected, it vanishes when the news does.We cross-checked with a computer code that the flux-balance laws for (2) E ab and (3) E ab are indeed obeyed for the M ij × 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] ≡ 0, referred to as dimensional identities, is not explicitly "apparent".
In our problem, we meet expressions like M kl nmL e p a e q b 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 [• • • ].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 kl nmn e p] a e q b δ ij δ lm δ pq δ kn = 0 . (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 +2 (0,4) 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 (1,3) (0,4) I ab +23 (0,3) I ab +23 (0,3) I ab +3u (2,2) 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.

The n 3 celestial charges at G 2 order
In the linear theory, the dressed n 0 Bondi aspects have been shown to provide the complete set of conserved charges for asymptotically flat spacetimes admitting a Bondi expansion [84].In this section, we will explicitly write the dressed n 3 Bondi aspects in the post-Minkowskian approximation, with terms linear in M and M ij and M × M ij and M ij × M kl interactions, and compute the corresponding charges.
For convenience, we shall express the Bondi data in terms of the corresponding Cartesian transverse tensors 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 and Q ± 3,L 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 [84] by For all terms considered, J L = 0 = Q − n,L for all L and n 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 × M ij , the explicit results are

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 × M ij and quadrupole-quadrupole M ij × 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)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 App. 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 3 Bondi charges, i.e., supermomenta, super-Lorentz and n = 2, 3 celestial charges, in eqs.(3.38)

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 × M ij harmonic metric.At quadratic order, they all take the form15  Here Q m (x), m ∈ N, is the Legendre function of the second kind (with branch cut singularity from −∞ to 1), and F (u) is a quadratic product of time derivatives of quadrupole moments, which vanishes when u −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) ∼ − 1 2 ln(x − 1) when x → 1 + .The perturbation equations (2.22) we must solve in our approach are of the type k µ ∂ µ 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 ∂J m (u + r, r) ∂r 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 → +∞ on both sides of eq.(A.13a), since Q m (x) tends to zero like 1/x m+1 when x → +∞ (see e.g., eqs.(A3) in [90]).We conclude that this polynomial must be identically zero: R m−2 (x) ≡ 0, hence our formula (A.3) is proven.Moreover, we see that (A.3) is equivalent to c mj W j−1 (x) . (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 τ = r(x − 1) as Since the integration range is actually a bound interval, we can substitute to Q m (x) its asymptotic expansion near x → 1 + for large enough r.This expansion is derived from the following suitable representation of the Legendre function: (m + j)! (m − j)!j! 2 x − 1 2 Expanding ln( x+1 2 ) when x → 1 + , commuting the sums, and resorting to standard resummation techniques including the identification of hypergeometric functions, 16 we get x − 1 2 j + 1 2 x − 1 2 m+1 +∞ j=0 (−) j j!(2m + j + 1)! (m + j + 1)! 2 x − 1 2 j .
(A.17) 16 In this instance, one may resort to the identity An alternative form of this expansion is given by eq.(4.6) in [91].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 → +∞ at u = const, For the expansion of I m in the near zone, i.e., r → 0 with u or t fixed, we refer the reader to App.A of [90] .
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, • • • , m.Those logarithms, however, have to disappear from the NU metric components by construction.Moreover, for the quadrupole-quadrupole interaction M ij × M ij , it turns out that all hereditary 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 × 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 × M ij interactions
Including only mass monopole and quadrupoles we have

.6b) Upon using the decomposition ζ a ζb = 1 2 δ a b + i 2 ab
and integrating by parts, we find 28) for arbitrary Diff(S 2 ) generators Y a = − ab ∂ b nL + γ ab ∂ b nL .The n = 2 dressed aspect is associated with the n = 2 "celestial" charges S (2) E ab (D a D b S + + ac D b D c S − ), defined for arbitrary scalars S ± (θ a ) on the 2-sphere.

8 ( 1
.34b) Finally, subtracting the right-hand sides from the left-hand sides of the balance equations (3.27b) and (3.30) for the dressed aspects (2) E ab and (3) E ab , respectively, lead to 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. ackowledge support from the Partenariat Hubert Curien within the Barrande mobility programme (project number 46771VC).

17
) ij + O(G 2 ) , (B.1a) D c C bc = −4Ge i b n j M (2) ij + O(G 2 ) , (B.1b) D a D c C bc = 4G n i n j M .(B.1c) is symmetric in (ab), hence the identity D [a D c C b]c = O(G 2).We then deduce from this 17 D a C 2 = −8G 2 e i a M It is useful to observe that, from C ab = e i a e j b H T T ij [see (2.45)], one has C ab = e k a e l b 2U kl − n m n p mn(k V l)np + 1 6 n p n q U klpq , .36)Summing up the different contributions in eqs.(3.26)and(3.29)[withthe undressed quantities therein given by eqs.(2.47b), (2.51a), and (2.51b), respectively], we obtain through computer calculation the explicit expressions for N i , (2) E ij and (3) E ij .Note that (n) E ij is traceless with respect to ⊥ ij because (n) E ab is traceless with respect to γ ab .Given the length of the resulting expressions for (3.36), we only provide their expressions for quadrupole-quadrupole interactions in Appendix B.
[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).