Multipole Expansion of Gravitational Waves: from Harmonic to Bondi coordinates

We implement the diffeomorphism transforming the metric of an isolated matter source in the multipolar post-Minkowskian approximation from harmonic (de Donder) coordinates to radiative Newman-Unti (NU) coordinates. At linearized order we obtain the NU metric as a functional of the mass and current multipole moments of the source, valid all-over the exterior region of the source. Imposing appropriate boundary conditions we recover the generalized Bondi-van der Burg-Metzner-Sachs residual symmetry group. To quadratic order, in the case of the mass-quadrupole interaction, we determine the contributions of gravitational-wave tails in the NU metric, and prove that the expansion of the metric in terms of the radius is regular to all orders. The mass and angular momentum aspects, as well as the Bondi shear, are read off from the metric. They are given by the radiative quadrupole moment including the tail terms.


Motivations
Gravitational waves (GWs), whose physical existence was controversial for years, were established rigorously in the seminal works of Bondi, van der Burg, Metzner and Sachs [1,2]. The Bondi-Sachs formalism describes the asymptotic structure near future null infinity of the field generated by isolated self-gravitating sources. This asymptotic structure was further elucidated thanks to the tools of the Newman-Penrose formalism [3] and conformal compactifications [4] leading to the concept of asymptotically simple spacetimes in the sense of Penrose [5]. Asymptotically simple spacetimes are now proven to follow from large sets of initial data which are stationary at spatial infinity, see e.g. the review [6].
Bondi coordinates or Bondi tetrad frames are defined from an outgoing light cone congruence with radial sections parametrized by the luminosity (areal) distance. A variant of these coordinates are the Newman-Unti (NU) coordinates whose radial coordinate is instead an affine parameter [7]. Bondi gauge and NU gauge share all essential features and can easily be mapped to each other [8,9]. Under the assumption of asymptotic simplicity, Einstein's equations admit a consistent asymptotic solution [10,11]. Such an asymptotic series is however limited to the vicinity of null infinity and it does not, in particular, resolve the source that generates the radiation.
For GW generation and applications to the data analysis of the GW events one needs the connection between the asymptotic structure of the field and explicit matter sources. This is achieved by the multipolar post-Minkowskian (MPM) expansion [35][36][37][38] which combines the multipole expansion for the field in the exterior region of the source with a nonlinearity expansion in powers of the gravitational constant G. The MPM formalism is defined in harmonic coordinates, also known as de Donder coordinates. At linear order the MPM expansion reduces to the linear metric as written by Thorne [39] and is characterized in terms of two infinite sets of canonical multipole moments, namely the mass and current multipole moments. A class of radiative coordinate systems exists such that the MPM expansion leads to asymptotically simple spacetimes for sources that are stationary before some given time in the remote past [36]. In such radiative coordinates, two infinite sets of radiative multipoles can be defined in terms of the canonical multipoles. They parametrize the asymptotic transverse-traceless waveform or, equivalently, the two polarizations of the Bondi shear.
In addition, the MPM formalism has to be matched to the post-Newtonian (PN) field in the near-zone and the interior of the source, which allows us to express the canonical multipoles in terms of the actual source multipoles and, furthermore, yields the radiationreaction forces caused by the radiation onto their sources [40][41][42]. The MPM-PN formalism was applied to compact binary systems and permitted to compute the GW phase evolution of inspiralling compact binaries to high PN order, see notably [43,44].
The main objective of this paper is to make explicit the relationship between Bondi expansions and the MPM formalism. The Bondi and NU gauges belong to the general class of radiative gauges in the sense of [36,45]. Here we will describe the construction of the explicit diffeomorphism transforming the metric in the MPM expansion from harmonic coordinates to NU coordinates. The diffeomorphism is perturbative in powers of G and, for each PM order, it is valid everywhere outside the source. After imposing standard boundary conditions, we find this diffeomorphism to be unique up to generalized BMS transformations [11-16, 24, 46-48], as we will cross-check in details. This allows us to transpose known results on the exterior MPM metric in harmonic gauge for a particular multipolar mode coupling and a given post-Minkowskian order to a metric in NU gauge, written as an exact expression to all orders in the radial expansion. As an illustration, we will explicitly derive the Bondi metric of the second-order post-Minkowskian (2PM or G 2 ) perturbation corresponding to mass-quadrupole interactions [38,49]. In particular this entails the description of GW tails within the Bondi asymptotic framework.
The rest of the paper is organized as follows. Section 1.2 is devoted to our notation and conventions. Section 2.1 recalls the harmonic-coordinates description of the metric in terms of canonical moments at linearized order. In Sec. 2.2 we present an algorithm implementing the transformation from harmonic coordinates to NU coordinates. In Sec. 3.1 we derive the NU metric as a function of the mass and current multipoles to linearized order. In Sec. 3.2 we impose standard boundary conditions for the asymptotic metric and naturally recover from our algorithm the gauge freedom associated with the BMS group. Notably, in Sec. 3.3, we obtain the Bondi mass aspect, the angular momentum aspect and the Bondi shear as multipole expansions parametrized by the canonical moments. In Sec. 4 we apply the algorithm to the quadratic metric (i.e. to 2PM order in the MPM formalism), focusing on the quadratic interaction between the mass monopole and the mass quadrupole. Explicit results on GW tails obtained in harmonic coordinates, are then conveyed into the NU metric in Sec. 4.1, to any order in the radial expansion. Finally we discuss in Sec. 4.2 the mass and angular momentum GW losses in the Bondi-NU framework at the level of the quadrupolequadrupole interaction. The paper ends with a short conclusion and perspectives in Sec. 5. Two Appendices gather technical details on the map between Bondi and NU gauges (A), and the all-order PM formulae for the coordinate change equations (B).

Notation and conventions
We adopt units with the speed of light c = 1. The Newton gravitational constant G is kept explicit to bookmark post-Minkowskian (PM) orders. We will refer to lower case Latin indices from a to h as 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 η µν = diag(−1, +1, +1, +1).
We denote Cartesian coordinates as x µ = (t, x) and spherical ones as (t, r, θ a ). Here, the radial coordinate r is defined as r = |x| and θ a = (θ, ϕ) with a, b, · · · = {1, 2}. The unit directional vector is denoted as n i = n i (θ a ) = x i /r. Euclidean spatial indices i, j, · · · = {1, 2, 3} are raised and lowered with the Kronecker metric δ ij . Furthermore, we define the Minkowskian outgoing vector k µ ∂ µ = ∂ t + n i ∂ i with ∂ i = ∂/∂x i , or, in components, k µ = (1, n i ) and k µ = (−1, n i ). In retarded spherical coordinates (u, r, θ a ) with u = t − r, we have k µ ∂ µ = ∂ r u . We employ the natural basis on the unit 2-sphere e a = ∂ ∂θ a embedded in R 3 with components e i a = ∂n i /∂θ a . Given the unit metric on the 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. 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. Introducing the covariant derivative D a compatible with the sphere metric, D a γ bc = 0, we have D a e i b = D b e i a = D a D b n i = −γ ab n i . Given a general manifold, harmonic/de Donder coordinates are specified by using a tilde: x µ = (t,x) or (t,r,θ a ). The metric tensor isg µν (x). Asymptotically flat spacetimes admit as a background structure the Minkowskian outgoing vectork µ = (1,ñ i ), the basis on the sphereẽ i a = ∂ñ i /∂θ a , etc. We define the retarded timeũ in harmonic coordinates asũ =t−r, such thatk µ = −∂ µũ .
Newman-Unti (NU) coordinates are denoted x µ = (u, r, θ a ) with θ a = (θ, ϕ). The metric tensor in NU coordinates is denoted as g µν (x), with all other notation, such as the natural basis on the sphere e i a and the metric γ ab , as previously. We denote by L = i 1 i 2 . . . i a multi-index made of spatial indices. We use short-hands for: the multi-derivative operator ∂ L = ∂ i 1 . . . ∂ i , the product of vectors n L = n i 1 . . . n i and x L = x i 1 . . . x i = r n L . The multipole moments M L and S L are symmetric and tracefree (STF). The transverse-trace-free (TT) projection operator is denoted ⊥ ijkl TT =⊥ k(i ⊥ j)l − 1 2 ⊥ ij ⊥ kl . Time derivatives are indicated by superscripts (q) or by dots.
2 From harmonic gauge to Newman-Unti gauge

Linear metric in harmonic coordinates
We work with the gothic metric deviation defined as h µν = |g|g µν − η µν and satisfying the de Donder (or harmonic) gauge condition∂ µ h µν = 0. The Einstein field equations in harmonic coordinates read as where˜ =˜ η is the flat d'Alembertian operator, and the right-hand side contains the matter stress-energy tensor T µν as well as the back-reaction from the metric itself, in the form of an infinite sum Λ µν of quadratic or higher powers of h and its space-time derivatives. We shall consider the metric generated by an isolated matter system, in the form of a non-linearity or post-Minkowskian (PM) expansion, labeled by G, Furthermore we consider the metric in the vacuum region outside the isolated matter system, and assume that each PM coefficient h µν n in Eq. (2) is in the form of a multipole expansion, parametrized by so-called canonical multipole moments. We call this the multipolar-post-Minkowskian approximation [35]. In the linearized approximation the vacuum Einstein field equations in harmonic coordinates read˜ h µν 1 =∂ ν h µν 1 = 0, whose most general retarded solution, modulo an infinitesimal harmonic gauge transformation, is [39] h 00 given in terms of symmetric-trace-free (STF) canonical mass and current multipole moments M L and S L depending on the harmonic coordinate retarded timeũ =t −r. Among these moments, the mass monopole M is the constant (ADM) total mass of the system, P i = M (1) i is the constant linear momentum and S i is the constant angular momentum. We can expand the linear metric in powers of 1/r using the formula (valid for arbitrary STF tensors M L ) A method has been proposed in [35] to compute each of the PM coefficients up to any order n, starting from the linear metric (3). Each of the PM approximation is then obtained as a functional of the canonical multipole moments M L and S L . The construction represents the most general solution of the Einstein field equations outside a matter source without any incoming flux from past null infinity. This is the so-called MPM formalism. The relation between the canonical moments and the source moments depending on actual source parameters is known [40][41][42].
In this paper we assume that the metric is stationary in the past in the sense that all the multipole moments are constant before some finite instant in the past, say M L (ũ) = const and S L (ũ) = const whenũ −T . Under this assumption all non-local integrals we shall meet (tails, . . . ) will be convergent at their bound in the infinite past. 1

Algorithm to transform harmonic to NU metrics
Consistently with the PM expansion (2), we assume that the NU coordinates are related to the harmonic coordinates by the following class of transformations where the PM coefficients U n , R n and Θ a n are functions of the harmonic coordinates (ũ,r,θ a ) to be determined, withũ =t −r.
The NU gauge 2 is defined by the following conditions: For computational reasons, it is more convenient to work with the inverse metric components, for which the NU gauge reads as The gauge is constructed such that 1) the outgoing vector k µ = −∂ µ u is null, the angular coordinates are constant along null rays k µ ∂ µ θ a = 0, and the radial coordinate is an affine 1 This assumption may be weakened to the situation where the source is initially made of free particles moving on unbound hyperbolic like orbits (initial scattering). In this case we would have M L (ũ) ∼ (−ũ) and S L (ũ) ∼ (−ũ) whenũ → −∞, and the tail integrals in the radiative moment, Eq. (46) below, would still be convergent for such initial state [50]. 2 The NU and Bondi gauges differ by a choice of the radial coordinate. See more details in [8] and in the Appendix A below. parameter on outgoing null curves, i.e. k µ ∂ µ r = 1. The strategy to construct the perturbative diffeomorphism is the following. From the NU gauge conditions (7), one finds the following constraints on the transformation laws (5), namelỹ Inserting the linear metric (3) this permits to solve for the linear corrections U 1 , R 1 and Θ a 1 , modulo an arbitrariness related in fine to BMS transformations. Then one uses to deduce g uu = −g rr + g ra g ua , g ua = g rb g ab , g ab = (g ab ) −1 to linear order. We can then read off, respectively, the Bondi mass aspect, the Bondi angular momentum aspect and the Bondi shear. To quadratic order one inserts the metric h µν 2 (x) in harmonic coordinates solving Eq. (1) to order G 2 , and obtain U 2 , R 2 , Θ a 2 and the NU metric to order G 2 . In the end we have to re-express the metric in terms of NU coordinates using the inverse of Eq. (5). This algorithm can be iterated in principle at any arbitrary order in powers of G.

Solving the NU gauge conditions
At linear order in G, the constraints (8) are equivalent to the following equations for the linear coefficients U 1 , R 1 and Θ a 1 , involving the directional derivative along the directioñ k µ = (1,ñ i ) of the Minkowski null cone: Notice that h ii 1 = 0 for the metric (3). Using the explicit form of the linearized metric (3)-(4) one readily obtains the most general solutions of those equations as where P is an irrelevant constant. We recognize the standard logarithmic deviation u = u − 2GM ln(r/P) + O(r −1 ) from harmonic to radiative coordinates; see e.g. [36]. Furthermore we have added the most general homogeneous solution of the differential equations (10) denoted by ξ µ . These are indeed the residual gauge transformations preserving the NU gauge (8), at linearized order, i.e. x µ → x µ + ξ µ with ξ µ = O(G 1 ). The linear gauge transformation, ξ = Gξ 1 takes the form where f , Q and Y a are arbitrary functions ofũ =t −r and the anglesθ a . Note that for later convenience, we made explicit into the expression of R 1 given by Eq. (11b) some constant monopolar and dipolar ( = 0, 1) contributions corresponding to a redefinition of the radial coordinate asr →r + G(M + 3ñ i P i ), thanks to the arbitrary function Q in Eq. (12). The metric in NU gauge is immediately obtained at linear order in G from the linear metric h µν 1 (and its trace h 1 = η µν h µν 1 = −h 00 1 ) as given by Eq. (3) together with the linear coefficients U 1 , R 1 and Θ a 1 as Note that the final result for the metric has been written in terms of the NU coordinates x µ . As a result the spatial metric g ab is given by a covariant tensorial expression on the sphere, involving the Lie derivative L Θ 1 γ ab = 2D (a Θ 1b) . To this end, we have written the leading contribution in the spatial metric g ab in terms of NU coordinates to linear order in G as where Γ a bc denotes the Christoffel symbol on the sphere. At linear order in G, we can equivalently replace the harmonic coordinates by the NU ones, as the correction will be at O(G 2 ). Plugging the results (11) into the metric (13), we find where we have posed e i a e j b = e i (a e j b) − 1 2 γ ab ⊥ ij . The last terms correspond to the freedom left in the metric, which is associated with the gauge vector (12), and are given by

Boundary conditions and the BMS group
A Bondi frame is defined from boundary Dirichlet gauge fixing conditions, which pick a specific foliation by constant u surfaces and a specific measure on the codimension 2 boundary. The boundary gauge fixing conditions when r → ∞ are where the first term in Eq. (17c) is the determinant of the metric on the unit sphere metric. Notice that Eq. (17c) not only fixes the measure on the sphere, but also requires that the shear which appears at order O(r 3 ) is trace-free, see the discussion around Eq. (2.5) of [8].
The metric (15)- (16) does not yet respect the boundary conditions (17). Thus one has to implement an infinitesimal transformation in order to achieve the gauge with appropriate asymptotic behavior. The first condition (17a) implies thatf = 0, henceḟ must only be a function of the angles θ a . The second condition (17b) implies thatẎ a = 0, i.e. that Y a also is only a function of the angles. To impose the last condition (17c), we note that the leading metric on the sphere γ ab already satisfies the leading behavior of (17c), i.e., its measure is that of a unit metric on the sphere. Therefore the leading term in Eq. (16) must be trace-free, thuṡ f = 1 2 D a Y a , which is consistent withf = 0. Similarly the next-to-leading term O(r 3 ) in g ab must also be trace-free, hence Q = 1 2 ∆f where ∆ = D a D a is the Laplacian on the sphere. Summarizing all these, we have The simplest choice that brings the metric (15) into the form (17) is of course obtained by setting f = Q = T = Y a = 0, which corresponds to using an asymptotically flat coordinate system at infinity. This choice is generally assumed in the perturbative approach to gravitational waves in harmonic coordinates, see e.g. [51]. However, after fulfilling all the conditions, i.e. the gauge conditions (6) and the asymptotic boundary conditions (17), we are still left with the infinitesimal coordinate transformations generated by the gauge vector field ξ µ BMS ≡ ξ µ , with components The latter play the role of symmetries of the space of solutions which are parametrized by a time-independent function T (θ a ) generating super-translations, and a time-independent vector Y a (θ b ) on the sphere generating super-Lorentz transformations. These form the celebrated generalized BMS algebra [15,16,24,46,47] (i.e., the smooth version of [11][12][13][14]). The modification of the metric under the BMS group reads 3 where we remind that Under the BMS transformation, the boundary fields at infinity transform according to Eq. (20) while the bulk metric transforms into itself. The transformation law of the boundary metric agrees with Eq. (2.20) of [24]. We note that the leading uu component of the metric is given by Eq. (20a) where the divergence D a Y a only involves the determinant of the metric on the unit sphere. This is consistent with Eqs.  (28)] is the = 0 and = 1 harmonics. In order to make this explicit, we decompose the function f into STF spherical harmonics where the STF coefficients f L are linear functions of u, and find

Bondi data to linear order
Finally, we shall write the metric (15) including the bulk terms in the form The sub-dominant contributions in 1/r in the metric (23) read as Note that Q ij 1 = 0 for k = 1. Therefore, at linear order in G the next order correction term in 1/r in the metric g ab beyond the shear C ab is absent. This is just a feature of the linear metric, since at quadratic order O(G 2 ) there is a well-known term quadratic in the shear.
To leading order when r → ∞ the metric (23) is defined by the so-called Bondi mass aspect m, angular momentum aspect N a and shear C ab (see e.g. [27,32,52]). These are functions of time u and the angles θ a . The mass and angular momentum aspects are given in terms of the multipole moments to linear order in G by In the next section we will work out the mass loss and angular momentum loss formulas forṁ andṄ a to quadratic order in G. But we already note thaṫ in agreement with the Einstein equation for the angular momentum aspect.
To define the shear we introduce the usual asymptotic waveform in transverse-trace-free (TT) gauge, given in terms of the multipole moments by (see e.g. [53]) where ⊥ ijkl TT is the TT projection operator. Then the shear is given by The first term is directly related to the usual two polarization waveforms at infinity. Posing H + = lim r→∞ (rh + ) and The second term in Eq. (28) comes from the BMS transformation as

Newman-Unti metric to quadratic order
At second order in G the perturbation (2) reads as 5 In the following we will denote h 2 ≡ η µν h µν 2 and the indices are lowered by the background Minkowski metric η µν . At second order in G, the NU gauge conditions (8) imply the following equations for the functions U 2 , R 2 , Θ a 2 , respectively, 4 We adopt for the polarization vectors ε i θ = e i θ and ε i See Appendix B for a formal generalization of these equations to any PM order.
In the following, we will show how the explicit solution for the quadratic metric in harmonic coordinates, i.e., solving the Einstein field equations (1) to order G 2 for some given multipole interactions, can be used as an input in our algorithm in order to generate the corresponding Bondi-NU metric.
The main features of the quadratic metric in harmonic coordinates are [38,49,54]: (i) the presence of gravitational-wave tails, corresponding to quadratic interactions between the constant mass M and varying multipole moments M L and S L (for 2); (ii) the mass and angular momentum losses describing the corrections of the constant ADM quantities introduced in the linear metric (M and S i ) due to the GW emission; 6 (iii) the presence of the non-linear memory effect. We investigate the effects (i) and (ii) in the subsections below but postpone (iii) to future work.

Tails and the mass-quadrupole interaction
In this subsection we construct the NU metric corresponding to the monopole-quadrupole interaction M × M ij , starting from the explicit solution in harmonic coordinates given by (see Appendix B of [38], or Eq. (2.8) of [49]) The metric is composed of two types of terms: the so-called "instantaneous" ones depending on the quadrupole moment M ij and its derivatives at timeũ =t−r, and the "hereditary" tail terms depending on all times from −∞ in the past untilũ. The tail integrals are expressed in Eq. (33) by means of the Legendre function of the second kind Q (with branch cut from −∞ to 1), given by the explicit formula in terms of the Legendre polynomial P : We recall that the Legendre function Q behaves like 1/x +1 when x → +∞, and that its leading expansion when y ≡ x − 1 → 0 + reads (with H = j=1 1 j being the th harmonic number) With the known harmonic metric (33) [or see below Eq. (45)], we apply our algorithm to generate the Bondi-NU metric. We focus on the case of the mass-quadrupole interaction M × M ij , keeping track of all instantaneous terms. Plugging h µν 2 given by Eq. (33) as well as h µν 1 and U 1 given in the previous section in the right-side of Eq. (32a), and retaining only the mass-quadrupole interaction we get We remark that, as an intermediate step to obtain (36), an instantaneous term of the form −2r −1 Mñ pq M pq has been equivalently written as the last term in the second line. Such reshuffling is necessary to integrate Eq. (36). To this aim, we have to divide inside the integral the combination of Legendre functions by a factor x − 1, 7 thus we must be careful about the bound of the integral when x → 1 + . However, we have the factorization which permits to immediately integrate Eq. (36) overr (while keepingũ fixed) with result In principle this is valid up to an homogeneous solution corresponding to a linear gauge transformation starting to order G 2 . It will be of the form −ξ u 2 = −f 2 where f 2 is a function ofũ =t −r andθ a . It thus takes the same form as the linear gauge transformation already introduced to order G in Eq. (11a). Hence, we can absorb f 2 into the redefinition of f through the replacement f → f + Gf 2 , and the solution (38) is the most general in our setup. Following the same procedure outlined above to compute U 2 , we obtain Having determined U 2 , R 2 and Θ a 2 we continue our algorithm and successively obtain the contravariant components g rr , g ra and g ab of the NU metric, and then its covariant components g uu , g ua and g ab , see Sec. 2.2. We consistently keep only the terms corresponding to the mass-quadrupole M ×M ij interaction. In the end we recall that we have to express the metric components in terms of the NU coordinates x µ = (u, r, θ a ) by applying (the inverse of the) coordinate transformation (5). In order to present the result in the best way we introduce the following tail-modified quadrupole moment as defined by [53] Such definition agrees with the expression of the radiative quadrupole moment parametrizing the leading r −1 piece of the metric at future null infinity. Restoring the powers of c −1 we see that the tail provides a 1.5PN correction ∼ c −3 to the quadrupole. Generally the radiative quadrupole moment is rather defined as the second-time derivative of M rad ij , see Eq. (76a) of [53]. But here, as we not only control the leading term r −1 but also all the subleading terms r −2 , r −3 , etc. in the expansion of the metric at infinity, it will turn out to be better to define the radiative moment simply as M rad ij . We find that the NU metric g uu to quadratic order G 2 for the mass-quadrupole interaction, including all terms in the expansion at infinity, reads We recover Eq. (15a) for the linear part, and we see that to quadratic order the tails nicely enter the metric only through the replacement of the canonical moment M ij by the radiative moment M rad ij defined by Eq. (40). In fact, with this approximation (neglecting G 3 terms), we can use either M ij or M rad ij in the second line of Eq. (41). Note that the last term of Eq. (41), involving a time integral over the radiative moment, is "exact" all over the exterior region of the source. The integral is convergent under our 8 We have changed the integration variable to z = r(x − 1). In previous formulae, it is convenient to decompose ln( x+1 x−1 ) = −ln( z 2P ) + ln(1 + z 2r ) + ln( r P ), where P is the constant introduced in Eq. (11). The first term gives the tail in the quadrupole (40), the second term gives the tail in the metric (41)-(43) and the third term is cancelled after reexpressing the metric in NU coordinates. assumption of stationarity in the past. Furthermore, this term is of order O(r −4 ) at null infinity where it admits an expansion involving only powers of r −1 . We have the regular expansion when r → +∞ for u = const: where −T is the finite instant in the remote past before which the multipoles are constant. Further processing we obtain the other components of the NU metric as g ua = Ge i a n j −2 Again we find some remaining tail integrals, but which rapidly fall off when r → ∞ and admit an expansion in simple powers of r −1 . Finally we conclude that the expansion of the NU metric at infinity is regular, without the powers of ln r which plague the expansion of the metric in harmonic coordinates. In intermediate steps of the computation, however, logarithmic divergences occur in the quadratic term, but they are cancelled by the expansion of the linear term taking into accountũ = u + 2GM ln r + O(G 2 ). The fact that the NU metric admits a regular (smooth) expansion when r → +∞ to all orders, without logarithms, is nicely consistent with the earlier work [36] which proved the property of asymptotic simplicity in the sense of Geroch and Horowitz [55], i.e., with a smooth conformal boundary at null infinity, for the large class of radiative coordinate systems, containing the Bondi and NU coordinates. Indeed, a crucial assumption in the proof of [36] as well as in our work, see Eq. (42), is that the metric is stationary in the past (for u −T ).
To second order in G, as already commented, we could still add to the construction some arbitrary homogeneous solutions of the equations for U 2 , R 2 and Θ a 2 , but the corresponding terms in the metric will have exactly the same form as those found to linear order in G, see Eqs. (16), and shown to describe with appropriate boundary conditions the modification of the metric under the BMS group.
From the results (41)-(43), one can easily deduce the mass and angular momentum aspects m and N a , and the Bondi shear C ab , for the case of the mass-quadrupole interaction to order G 2 . As expected the Bondi data are entirely determined by the radiative quadrupole moment (40). Recalling the expression of the metric in the NU gauge as given by Eqs. (61) and (64) in Appendix A, where N a is defined according to the convention of [20], we get N a = 3e i a n j ε ijk S k + 2 We have added in the angular momentum aspect the linear contribution due to the total constant (ADM) angular momentum or spin S i , as read off from Eq. (25b). Notice that the difference between the Newman-Bondi and Bondi radii is a term quadratic in C ab , see Eq. (60). This term is thus quadratic in the source moment M ij , and so, for the mass-quadrupole interaction M ×M ij considered in this section, there is no difference between the NU and Bondi gauges.
We can in principle generalize the latter results to multipole interactions M × M L and M × S L (with any 2), starting from the known expressions of tail terms in the metric in harmonic coordinates: 9 Here the ellipsis refer to many non-tail contributions, in the form of instantaneous (i.e., local-in-time) terms depending on the multipole moments only at timeũ. Considering the previous results we can conjecture that the mass and angular momentum aspects will take the same form as in Eqs. (25) but with the canonical moments M L and S L replaced by the radiative moments M rad L and S rad L [56] M rad where the constants are given by (with H = j=1 More work would be needed to generalize our algorithm in order to include any multipole interactions M × M L and M × S L (especially instantaneous ones).

Mass and angular momentum losses
Taking the angular average of the mass aspect m we obtain the Bondi mass M B ≡ dΩ 4π m. At this stage, we find from Eqs. (44a) or (25a) that the Bondi mass just equals the ADM mass M ADM ≡ M . This is because we have not yet included the mass loss by GW emission which arises in this formalism from the quadratic interaction between two quadrupole moments, say M ij × M kl , as well as higher multipole moment interactions. The losses of mass and angular momentum are straightforward to include in the formalism, starting from the known results in harmonic coordinates.
The terms responsible for mass and angular momentum losses (at the lowest quadrupolequadrupole interaction level) in the harmonic-coordinate metric are (see e.g. Eq. (4.12) in [54]): where again, the ellipsis denote many instantaneous (local-in-time) terms, in contrast with the non-local time anti-derivative integrals over the multipole moments in Eq. (48). Importantly, the ellipsis in Eq. (48) also contain another type of non-local terms that are associated with the non-linear memory effect, but which we shall not discuss here. The complete quadrupole-quadrupole interaction M ij × M kl has been computed in harmonic coordinates in [54], including the description of the various GW losses and the non-linear memory effect. We thus apply our algorithm to generate the corresponding mass and angular momentum losses in the NU metric. In this calculation we only keep track of the non-local-in-time (or "hereditary") integrals, and neglect all the instantaneous terms. Furthermore, as we said we do not consider the memory effect, which is disconnected from GW losses (see e.g. [54]). Finally we are restricted to the quadrupole-quadrupole interaction, as in Eq. (48).
Looking at the second-order equations (32) we see that we are just required to solvẽ We obtain successively (changing consistently harmonic to NU coordinates) We find no such hereditary terms in R 2 . The logarithmic term in U 2 corrects the light cone deviation at linear order as given by Eq. (11a). The corresponding contributions in the NU metric follow as Combining this with previous results (44a) or (25a) we obtain the mass aspect which is now accurate enough to include the physical GW mass loss Hence the Bondi mass M B = dΩ 4π m reads (where M is the constant ADM mass) The mass loss in the right-side is characterized by the hereditary (or "semi-hereditary") 10 non-local integral, in contrast with the instantaneous contributions indicated by dots. Such instantaneous terms will be in the form of total time derivatives in the corresponding flux balance equation, and may be neglected in average over a typical orbital period for quasiperiodic systems. Thus the averaged balance equation reduces to which is of course nothing but (with this approximation) the balance equation corresponding to the standard Einstein quadrupole formula.
In a similar way we obtain the angular momentum aspect and Bondi shear as N a = 6e i a n j 1 2 The Bondi angular momentum is defined from the angular momentum aspect by As shown in [32], this quantity requires a prescription for α which is fixed to α = 1 in [14,20,24,28], α = 0 in [22,57] or α = 3 in [27]. Since the α-term gives instantaneous terms as well as higher order terms, we can simply ignore it for this computation. Hence we have Upon averaging this leads to the usual quadrupole balance equation for angular momentum 11 Note that the discussion of the GW losses in the linear momentum (or recoil) and the center-of-mass position would require the coupling between the mass quadrupole and the mass octupole moments, which is outside the scope of the present calculation.

Conclusion and perspectives
In this paper we have shown how to implement practically the transformation of the metric of an isolated matter source in the MPM (multipolar post-Minkowskian) approach from harmonic (de Donder) coordinates to Bondi-like NU (Newman-Unti) coordinates. This is of interest because the asymptotic properties of radiative space-times are generally discussed within the Bondi-Sachs-Penrose formalism, while the connection to the source's properties is done by a matching procedure to the source using the MPM expansion.
In particular we obtain explicit expressions for the NU metric valid at any order in the radial distance to the source (while staying outside the domain of the source), expressed in terms of the canonical mass and current multipole moments. Under the assumption of stationarity in the remote past, we prove that the NU metric (for particular multipole moment couplings) admits a regular expansion at future null infinity. This is consistent with the fact that the MPM expansion satisfies the property of asymptotic simplicity [36].
On the other hand the canonical moments are known in terms of the source's parameters to high PN (post-Newtonian) order. Our approach permits to rewrite explicit results derived in harmonic coordinates using the MPM approximation into the Bondi-Sachs-Penrose formalism for the asymptotic structure, including the notions of Bondi shear, and mass and angular momentum aspects. In particular, we recover from our construction the generalized BMS (Bondi-van der Burg-Metzner-Sachs) residual symmetry group leaving invariant the NU metric under appropriate boundary conditions at future null infinity. 12 To non-linear order our construction is in principle valid for any coupling between the canonical moments. In this paper we have worked out the coupling between the mass and the 11 The angular momentum aspect itself satisfies, see also Eq. quadrupole, including the contributions due to non-local (hereditary) tail effects but also all local (instantaneous) terms. Including the non-local (semi-hereditary) terms arising from the coupling between two quadrupoles, we obtain the mass and angular momentum losses due to the GW emission through the expressions of the mass and angular momentum aspects. However we ignored all the instantaneous terms in the quadrupole-quadrupole metric, as well as the contributions from the non-linear memory effect. In future work we intend to thoroughly investigate the quadrupole-quadrupole interaction in our framework, and in particular discuss the occurrence of the non-linear memory effect, thereby contrasting the perspective from approximation methods in harmonic coordinates with that from asymptotic studies in Bondi-like coordinates confined at future null infinity.
Bondi gauge and Newman-Unti gauge differ by a choice of radial coordinate [8]. They both admit identical asymptotic symmetry groups, phase spaces and physical quantities [8]. We denote in both coordinate systems the angular coordinates as θ a and the coordinate labelling the foliation of null hypersurfaces as u. Let us refer to r B as the Bondi radius and r NU as the Newman-Unti radius. The Newman-Unti radius r NU is the affine parameter along the hypersurfaces of constant u defined from g r NU u = −1 while the Bondi radius is the luminosity distance such that ∂ r B [det(g ab )/r 4 B ] = 0. The relationship between the radii is given by [8] r NU = r B + ∞ r B dr g r B u + 1 , r B = det g ab det γ ab For large radii, we have There is a modification at subleading order and at subsubleading order. We deduce that C ab and m can be read off from the metric in Newman-Unti gauge as g NU ab = r 2 NU γ ab + r NU C ab + O(r 0 NU ) , with m NU = m + 1 16 ∂ u (C ab C ab ). Instead, where r is either r B or r NU . In the convention of [20], the angular momentum aspect N a is read in Bondi gauge from We deduce from Eq. (62) that it is read in Newman-Unti gauge from

B Equations for any PM order
At any given PM order p ∈ N, the NU gauge conditions (8) To derive the equation for R p , one formally writes |g| = 1 +