Electromagnetic fields in compact binaries: post-Newtonian wave generation and application to double white dwarfs systems

The aim of this work is twofold: (i) to properly define a wave-generation formalism for compact-supported sources embedded in Einstein-Maxwell theory, relying on matched post-Newtonian and multipolar-post-Minkowskian expansions; (ii) to apply this formalism (which is valid for any type of post-Newtonian sources) to the case of two stars with constant and aligned magnetic dipoles, by computing the fluxes of energy and angular momentum to the next-to-leading order, as well as the gravitational amplitude modes. Assuming eccentric orbits, we derive the evolution of orbital parameters, as well as the observables of the system, notably the gravitational phase for quasi-circular orbits. Finally, we give some numerical estimates for the contribution of the magnetic dipoles for some realistic systems.


I. INTRODUCTION
High accuracy analytic predictions are a milestone of the signal analysis for gravitationalwave (GW) detectors, notably when it comes to observing the long inspiral phase of coalescing compact objects.If such objects, typically neutron stars (NS) and white dwarfs (WD), are far from being the main source of events of the current LIGO-Virgo-KAGRA network [1], they will be a major source for the next generations of detectors, such as the Laser Interferometer Space Antenna (LISA) [2], or the Einstein Telescope (ET) [3].Indeed, those instruments will resolve tens of thousands of galactic binaries, composed primarily of WD and/or NS [4,5].It is thus crucial to provide accurate analytic templates for those sources, notably taking into account physical effects beyond the commonly used spinning point-particle approximation.Among those effects are for instance the tidal response of the stars, that have been treated to high accuracy in [6,7], or the presence of strong magnetic fields, that can be as intense as 10 9 G for white dwarfs and 10 12 G for neutron stars [8,9].Such fields can induce an orbital decay, and shift the frequency of the gravitational wave emitted by the binary [10][11][12][13][14][15][16].The LISA experiment will resolve a large number of double WD [17,18], and some of those systems will even be used for the calibration of the instrument [19][20][21].It seems therefore important to provide accurate analytic waveform templates that incorporate the electromagnetic structure of WD and NS, which is the purpose of this work.
The templates that will be used for the detection and characterization of the inspiral phase of galactic binaries are mainly based on post-Newtonian results [22], relying on both weakfield and slow-motion approximations (see e.g.[23][24][25] for reviews).In such a framework, it is customary to distinguish between a "conservative" sector, describing the motion of the two objects via a set of Noetherian quantities, and a "dissipative" one, depicting the fluxes at spatial infinity.Both sectors are linked by the famous flux-balance equations, that describe the loss of energy E and angular momentum J i of the system as dE dt = − F and where F and G i are respectively the energy and angular momentum fluxes, as detected by an observer at spatial infinity and brackets • denote the usual orbital average.We let the interested reader refer to [26] for the flux-balance equations associated with the other Noetherian quantities, i.e. linear momentum and center-of-mass position.In the case of a binary system bearing electric charges, the left-hand sides of the flux-balance equations have been widely studied, see e.g.[27][28][29][30][31][32] but, to the best of our knowledge, a proper treatment of the right-hand sides is missing.Moreover, and as previously argued, going beyond the simple electric charge case will be important for future GW detectors.In this spirit, we have achieved in [33] the construction of the conservative sector including electromagnetic (EM) effects, and computed the Noetherian energy and angular momentum at next-toleading orders for a binary system of stars bearing magnetic dipoles.The present work is the continuation of this effort, and we aim at properly defining the corresponding dissipative sector, and build a consistent wave-generation formalism, taking electromagnetic effects into account.
The computation of PN gravitational waveforms relies on matched post-Newtonian and multipolar-post-Minkowskian expansions.This "PN-MPM" framework was developed in [34][35][36][37][38] and summarized e.g. in [23].It enabled to reach high precision in the analytic determination of the gravitational waveform: the phase is currently known at 4.5PN order [39], the gravitational flux for generic orbits and the dominant gravitational modes, at 4PN order [40], and the other modes, at 3.5PN order [41][42][43] (we recall that the nPN order corresponds to a (v/c) 2n correction beyond the leading order, where v is the typical velocity of the virialized binary system and c the speed of light).Besides the PN-MPM formalism, two frameworks have also been developed to deal with the wave generation, namely the direct integration of the relaxed equations [44] and effective field theory [45,46].They both have reached the 2PN precision for point-particle [47,48].The aim of the present work is thus to construct a proper PN-MPM formalism including electromagnetic effects, and to apply it to derive the waveform emitted by a binary of stars bearing magnetic dipoles.
The plan of this paper is as follows.In Sec.II, we present the construction of the PN-MPM wave-generation formalism including electromagnetic effects, and compute the first non-linearities (the so-called tails) in Sec.III.Then, we apply our formalism to the case of two spinless stars bearing constant and aligned magnetic moments in Sec.IV, and derive the various fluxes, the evolution of orbital parameters for elliptic trajectories, as well as the gravitational phase and modes for quasi-circular orbits.We also give numerical estimates for those quantities, before concluding in Sec.V. App.A presents the construction of the electromagnetic fluxes.Lengthy expressions are displayed in App.B, and stored in an ancillary file [49].

II. WAVE-GENERATION FORMALISM
The aim of this section is to build a wave-generation formalism in the Einstein-Maxwell theory, directly inspired from the usual PN-MPM framework used in GR [34][35][36][37][38] (see [23] for a review).In such framework, the behavior of the fields at future null infinity is parametrized by some radiative multipole moments, from which the fluxes of energy and angular momentum, as well as the gravitational amplitude modes, are derived.Those radiative moments are linked to the source multipole moments, describing the matter distribution, by taking into account the non-linearities arising during the propagation of the waves towards spatial infinity.
As clear from the analysis of the generic structure of the fields in the exterior zone, performed in Sec II B, three steps are required to construct a proper PN-MPM wave-generation formalism: (i) expressing the source moments in terms of the matter distribution, as done in Sec.II C; (ii) taking into account the non-linearities via an iteration scheme and extracting the radiative moments, as discussed in Sec.II D; (iii) relating those moments to the modes and fluxes, as presented in Sec.II E.

A. Equations of motion
In the following 1 , we consider a generic matter distribution with compact support in an Einstein-Maxwell framework.The metric perturbation h µν ≡ √ −gg µν − η µν and EM field 2   A µ are governed by the coupled fields equations [33] h µν = 16πG c 4 τ µν and together with the gauge conditions where g denotes the determinant of g µν , and α EM = µ 0 e 2 c 4π is the fine-structure constant, i.e. a dimensionless number.The source terms of the wave equations (2.1), τ µν and ι µ , play the role of the Landau-Lifshitz pseudotensor in GR, and read We have separated the sources between the descriptions of the matter distribution, {T µν , J µ }, and the non-linear interaction terms {Λ µν , Φ µ }.The contribution Λ µν can be further split between the metric self-interaction term Λ µν g and the interplay between the metric and the EM field, dubbed Λ µν EM .The metric self-interaction Λ µν g is the usual GR one and can found e.g. in Eq. ( 24) of [23].The two other contributions read explicitly [33] Λ where we have introduced the usual field strength As for the matter sources {T µν , J µ }, their precise expressions do not play a role in the present section.Indeed our formalism applies to any system with compact-supported post-Newtonian matter distribution.Nevertheless, when applying the techniques developed in the 1 The conventions employed throughout this work are as follows: we use a mostly plus signature, the Minkowski metric being η µν = (−, +, +, +); greek letters denote spacetime indices µ, ν, . . .= (0, 1, 2, 3) and latin ones, purely spatial indices i, j, . . .= (1, 2, 3); bold font denotes three-dimensional vectors, e.g.present section to a concrete example in Sec.IV, we will need to specify a source modeling.Following what was done in [33], we will consider a compact binary of spinless point particles, dressed with EM dipoles.Their dynamics is free and will be specified later on.We have where ] is the 3-dimensional Dirac delta distribution, locating the matter distribution on the worldline of the particle A, y A (t) (thus the sources are indeed compact supported), A the usual 3-dimensional velocity), and the subscript A means that the quantity is regularized in the worldline of particle A, using the techniques of [50].The physical quantities describing the source are the individual masses m A and the EM dipoles, gathered into a Lorentz tensor where {q i A , µ i A } are the electric-and magnetic-type dipoles and ε ijk the Levi-Civita symbol.Last, but not least, note the difference in structure between the two gauge conditions (2.2): while harmonic condition on the metric is linear, the U(1) gauge condition on the EM field has a non-linear sector, involving h µν .This fact will be crucial when explicitly implementing the MPM iteration scheme in Sec.III A.

B. Generic structure of the fields in the exterior zone
The equations of motion (2.1) are formally solved by the application of a retarded threedimensional Green function but those integrals are too complicated to be solved directly.As we are interested with radiation towards spatial infinity, we will consider the problem in the exterior zone, i.e. the region of space outside the matter distribution.In such region of space, both fields h µν and A µ coincide with their multipolar expansions, M (h µν ) and M (A µ ) (as those are formal series).Moreover, we have by construction T µν = J µ = 0 in this region.As multipolar expansions commute with products and derivatives, we then have to solve a simpler, multipolar expanded, version of the equations of motion (2.1) (2.8b) Solving the vacuum equations (2.8) is simpler than computing the integrals (2.7).Nevertheless, the r → ∞ expansion performed induces spuriously divergent integrals, that have to be regularized.Let us emphasize that those divergences are not physical, but simply a direct consequence of the formal multipolar expansion.In order to deal with them, we employ the usual Hadamard regularization scheme [50,51], introducing an unphysical (length) scale r 0 , thus solving (2.8) as3 Note that, together with the regularized propagator acting on the non-linear sources, we have introduced homogeneous solutions, h µν hom and A hom µ .Imposing that those homogeneous solutions are regular at spatial infinity (as we recall that we are working in the exterior zone), and the usual no-incoming radiation condition (i.e. that no wave can "come from infinity"), they have to be written as where the quantities F µν L and G L µ are yet unspecified functions, STF in L. It is easy to see that those homogeneous solutions indeed satisfy the linear wave equations in vacuum h µν hom = A hom µ = 0.The formal separation of the solution (2.9) bears a physical meaning.On the one hand, the homogeneous solutions h µν hom and A hom µ represent the linear wave emitted by the source.As such, they encrypt information on the matter distribution, as presented in Sec.II C. On the other hand, the regularized propagators acting on the non-linear sources Λ µν and Φ µ describe scattering and diffusion processes happening during the propagation of the wave towards spatial infinity.They are discussed in Sec.II D.

C. Expression of the source moments
As the source moments encompass the information on the matter system, we need to express the linear (homogeneous) solution in terms of the matter distribution.The aim is twofold: (i) decomposing the quantities F µν L and G L µ entering (2.10) as simple, irreducible STF source (and gauge) moments; (ii) matching those moments to the matter distributions τ µν and ι µ .
The wave equation for h µν in terms of τ µν (2.1) is the same as the one used in pure GR.Indeed, in the metric sector, the only difference between this work and usual GR computations is the explicit expression of τ µν in terms of the fields.Therefore, the STF reduction and matching of the moments in terms of τ µν will be exactly the same as the usual one, derived in [38].The linear solution is thus parametrized by six sets STF moments as [34] h This generic homogeneous solution is split in two sectors.The canonical linear metric h µν can bears the two propagating degrees of freedom of GR, and is written in terms of two sets of source moments, {I L , J L } [54][55][56].Its expression can be found e.g. in Eq. ( 36) of [23].The second sector, ∂ϕ µν ≡ ∂ µ ϕ ν + ∂ ν ϕ µ − ∂ λ ϕ λ η µν , is nothing but the generator of a diffeomorphism transformation, parametrized by four sets of gauge moments {W L , X L , Y L , Z L }, as explicitly displayed e.g. in Eq. ( 37) of [23].The matched expressions of those moments in terms of τ µν are to be found in Eqs. ( 123) and (125) of [23].Note that the monopole I and the dipoles I i and J i are respectively the conserved Arnowitt-Deser-Misner (ADM) mass, linear momentum and angular momentum.We will operate in the center-of-mass frame, in which I i = 0, and denote M = I.

Irreducible decomposition of the EM field
We need to perform the decomposition of the linear EM field (2.10b) in irreducible STF tensors, while ensuring that the gauge equation (2.2) is respected.Hereafter, for a quantity Q L , we use the shortened notation Following Eq. (5.3) of [38], let us decompose where (2.15b) Then, we impose the gauge condition As we are investigating the linear solution, we can drop the O(G 2 ) remainder for now.Nevertheless, this remainder will be crucial during the iteration procedure, as explicitly shown in Sec.III A 2. The linear gauge condition thus translates into the relation Next, it is common knowledge that an vector field obeying a U(1) symmetry only propagates two degrees of freedom, whereas the linear solution (2.18) contains three sets of moments.Thus, in the same spirit as the splitting of h µν hom (2.11),A hom µ can be rearranged into a "canonical" EM field, parametrized by two sets of STF source moments, {E L , B L }, of respective electric and magnetic parities, and a U(1) generator, parametrized by a set of STF gauge moments, {K L } where we pose The signs and the powers in c have been fixed to recover the usual results for constant charges and dipoles [57].Comparing with Eq. (2.18), one can read Let us note that for ℓ = 0, the linear gauge condition (2.16) simply yields Ė = 0, which is nothing but the conservation of the electric monopole: the electric charge defined in our formalism is of the ADM type.In order to reinforce this fact and avoid confusion, we will denote in the following Q = E.

Matching the source and gauge moments
Once the EM field has been decomposed in irreducible STF moments, one needs to link those moments to the parameter of the matter distribution, through the powerful matching procedure [38].This procedure relies on a few hypothesis, presented e.g. in Sec.4.1 of [23]: (i) the problem can be split into two zones, that partially overlap; (ii) the fields admit a well-defined PN expansion in the near-zone; (iii) the fields admit a well-defined multipolar expansion in the exterior zone.Under those hypothesis, one can match the two expansions in the overlapping zone, and derive an explicit expression for the moments in terms of the source of the wave equation.
Hopefully, the structure of our problem is identical to the pure GR one, and all the hypothesis are verified.Indeed, the first hypothesis is satisfied as long as the matter distribution has a compact support, which is one of the assumptions of this work.The second one has been verified by the study of the conservative sector, presented in [33].As for the last one, it directly follows from Eq. (2.8).Therefore, one can apply the method exposed in [23,38] to the EM field and derive where the bar denotes a PN expansion (i.e. a formal c → ∞ series), and we have introduced Those expression are in full agreement with the results of [58], that were derived in the framework of non-relativistic electromagnetism (corresponding to setting Φ µ = 0).

D. Multipolar-post-Minkowskian iteration scheme and radiative moments
Once the linear solutions h µν hom and A hom µ are decomposed in irreducible STF moments, one can turn towards the non-linear sectors of the solutions (2.9).In order to do so, we proceed via a post-Minkowskian expansion where the linear orders h µν (1) and A (1) µ are naturally given by the homogeneous solutions (2.11) and (2.19).Note that we combine multipolar and post-Minkowskian expansions, hence the name of the iteration scheme.

Sketch of the iteration procedure
The procedure is conceptually very simple: the vacuum equations (2.8) are also expanded in a post-Minkowskian fashion, and solved order by order in G. Obviously, we cannot provide here the full solution: for each practical computation, one has to determine by dimensional analysis which interactions will play a role and compute them using techniques developed notably in [37,[59][60][61].Note that the main technical difference between our case and the pure GR computations is due to the non-linear nature of the gauge condition for A µ .Nevertheless, verifying the gauge condition on the nth PM order A (n) µ only requires the knowledge of lower-order solutions {h µν (m<n) , A (m<n) µ }, already required to compute A (n) µ .Therefore this non-linear gauge fixing is not a conceptual difficulty by itself.We provide in Sec.III a practical example of such a computation, with the determination of the first nonlinearities, namely the tails.An interesting fact is that, due to the form of the interaction term Λ µν EM (2.4a), the EM corrections to the gravitational moments {I L , J L } involve at least two EM moments.For instance, the lowest order EM interactions correcting the mass quadrupole I ij are given by Q × E ij and E i × E j , which both enter at a relative 2.5PN order (as shown by a rapid dimensional analysis).
Note also that, when dealing wih pure GR, it is customary to first let apart the sourceto-source interactions and to focus on computing the interactions between gauge and source moments.Indeed, such interactions can be recast in the form of a physically equivalent description in terms of only two sets of canonical moments [60,62,63].As the structure of the linear EM field (2.19) is very similar to the one of the linear metric (2.11), it is expected that a similar physical equivalence exists.Nevertheless, the difference between source and canonical moments enter at relative 2.5PN for the gravitational moments [62], and also for the EM ones.Indeed, rapid dimensional and parity analysis show that the lowest-order gauge corrections to E i (resp.B i ) are given by the interactions M × K i , W × E i and Q × W i (resp.W × B i and K × J i ), which are all multiplied by a factor G/c 5 .Such precision is beyond the scope of this work, and the gauge moment will play no role in the practical computations of Sec.IV.We thus let the proper definition of canonical EM moments to future works.

Asymptotic structure of the fields
Once multipolar-post-Minkowskian expanded metric and EM fields have been constructed up to the required order, one has to extract the radiative moments.Again, the procedure for the metric is similar to the pure GR case, and we let the reader refer to Sec. 3 of [23] for details.We will here concentrate on the EM case, and extract its leading behavior at future null infinity r → ∞, with u = t − r/c constant.Using the relation the behavior of the linear solution (2.20) is easily extracted where the moments are evaluated at u.We recall that the gauge moments {K L } represent a U(1) symmetry and, as such, play no physical role.If the linear solution is perfectly regular at infinity, this will not be the case of the non-linear interactions, that will bear (powers of) logarithms.Indeed, the most generic structure of a MPM iterated solution at future null infinity is of the type [23] where p ≥ 0 and the "small remainder" translates the polylogarithmic behavior at lower orders.In order to properly define radiative moments, we need to eliminate those logarithms.

Radiative moments
In pure GR, the study of the asymptotic structure of the space-time has revealed that this polylogarithmic behavior is an artifact of the harmonic coordinates, and there exist "radiative" coordinate systems in which the metric is free of it [64,65].The most famous ones are the Bondi coordinates system [66,67] and the Newman-Unti one [68], but the class of radiative coordinates is quite large [69][70][71].Naturally, there exist generic constructions relating harmonic and Bondi-Newman-Unti coordinate systems exist [35,72,73] (see [74] for a practical implementation at the state-of-the art 4PN precision).
Our claim here is that there exists a radiative coordinate system (T, X i ) in which both the metric and the EM field are free of polylogarithmic structure.
with naturally U = T − R/c.Such claim is corroborated by the analysis of the Sec.III, were we show that, under the coordinate change (3.12) both the metric and the EM field takes the form (2.29).In such coordinate system, one can project the metric H µν in a traceless and transverse (TT) gauge [56] where we introduced the TT projector The two sets of moments {U L , V L } entering the asymptotic metric are naturally dubbed mass and current radiative moments.
As for the EM field, let us project it on a transverse basis, defined by A T 0 = N i A T i = 0, as where similarly, we have introduced two sets of radiative electric and magnetic moments, {Q L , H L }. Looking back at Eq. (2.27), it comes as expected where the remainder encrypts all non-linearities, and notably the tails computed in Sec.III.

Gravitational sector
The detailed derivation of the energy and angular momentum fluxes associated with the asymptotic metric (2.30) can be found in [56], and they read L U From the asymptotic metric (2.30), one can also extract the observable gravitational modes, h ℓm , by projecting H TT ij onto the basis of polarization {+, ×} and on spin-weighted spherical harmonics, Y ℓm −2 (following the conventions of [41,63]) as (2.34) In the particular case of planar orbits, as will be considered in Sec IV, the gravitational modes are linked to the radiative moments by the relations [41] ) Introducing a fixed orthonormal basis (n 0 , λ 0 , l) where l is normal to the orbital plane, together with m 0 = (n 0 + iλ 0 )/ √ 2, the projector α ℓm L is explicitly given by where the overbar denotes complex conjugation.

Electromagnetic sector
As discussed in App.A 1, the energy loss can be expressed by the angular average flux of the Landau-Lifshitz pseudotensor τ µν .At leading order in 1/R, the EM sector of this tensor reduces to the usual Poynting vector Π i = c F T ij F T j0 .The EM energy flux then reads where we used the relation ) together with the transversity.Plugging the asymptotic expression for the EM field (2.31), and performing the angular integrals with the help of the formula (A8), it comes This result is consistent with the computation of [46], and the canonical Larmor formula [57] is recovered in the case of a simple electric dipole.Moreover, to emphasize the gauge invariance of the flux, a derivation relaxing the transverse gauge condition is presented in App.A 2 and fully agree with our result (2.38).
As for the angular momentum flux, it requires more care, just as in the gravitational case (see e.g.[26,[75][76][77][78][79] for subtleties linked to the definition of an angular momentum flux in GR).As discussed in App.A 1, it is also defined with respect to the EM sector of the Landau-Lifshitz pseudotensor (2.4a) as (2.39) But, as clear from the R 3 factor, its evaluation requires the knowledge of the sub-dominant O(1/R 2 ) order in the asymptotic expansion of the field (2.31).This is quite similar to what is happening in pure GR, see e.g.[56], and calls for a full BMS-type analysis of the asymptotic structure of the fields.Instead of performing such tedious computation, we follow the spirit of [26] and provide in App.A 3 a derivation leading to the result (2.40) Note that, in the case of a rotating electric dipole, we recover the results of [80,81].
Together with the gravitational fluxes (2.33), those fluxes enter the balance equations for the energy and angular momentum (1.1), that allow to derive the secular evolution of orbital quantities, see Sec.IV C. In the case of quasi-circular orbits, they are key ingredient to derive the gravitational phase φ = dt ω (where ω is the gravitational frequency), as presented in Sec.IV D. Interestingly, the two fluxes (2.38) and (2.40) start formally at the same order than the gravitational ones, namely 2.5PN.This is a strong difference with the results of the conservative sector, where the EM effects enter formally as 1PN corrections to the Noetherian quantities [33].

III. COMPUTATION OF THE LEADING TAIL EFFECTS
The leading corrections to the fluxes derived in the previous section are the so-called tail interactions.Those interactions depict the scattering of a propagating moment (I L , J L , E L or B L ) onto the static structure of the space-time, described by the constant (ADM) mass M. The gravitational tails M × I L and M × J L are known since a long time [37,[82][83][84] (see also [85] for a derivation using Schwarzschild coordinates, and [74] for a computation in a radiative coordinate system), and we will not redo it here.Instead we will focus on the leading EM tails, namely the interaction of the EM ℓ = 1 moments E i and B i with the ADM mass M, and will explicitly apply the MPM iteration scheme pictured in Sec.II D 1 to those interactions. 4.The electro-magnetic tail effect The MPM iteration method [37] is a two-steps procedure, described in details e.g. in Sec. 2 of [23].First, one computes a particular solution of the sourced wave equation, by applying the Hadamard-regularized propagator.But, due to the regularization procedure, this particular solution may not obey the gauge condition, even if the source does 5 and one needs to correct for it.Therefore, the second step is to add a proper homogeneous solution that corrects for the violation of the gauge condition.The sum of those two solutions is the correct value of the field, from which the radiative moments can be extracted (after the suitable coordinate change, as explained in Sec.II D 3).

Computation of the particular solution
In what follows, we aim to compute the quadratic interactions M × E i and M × B i , by solving the PM-expanded vacuum equation (2.8b).The first step is naturally to construct the non-linear source M(Φ µ ), by injecting the MPM expansions of the fields into the expression of Φ µ .Selecting only the M sector of the linear metric (2.11) and the appropriate sectors of the linear EM field (2.20), we need to consider ) where we recall that we employ the shortcut (2.12), and have used the notation of [87], which amounts here to B i|j ≡ ε ijk B k .Those MPM-expansions are injected into the nonlinear source term (2.4b) developed to quadratic order where indices are operated with the Minkowski metric and h = h µ µ is the trace of h µν .It is easy to see that the equations we need to solve are of the form and do not involve logarithmic dependencies, so we can apply the formulas displayed in App. 1 of [59] to compute the particular solutions which explicitly read6 where the Legendre functions Q ℓ are defined with a branch cut on (−∞, 1) as with P ℓ the Legendre polynomials.

Implementation of the gauge condition
Due to the presence of the regularization kernel (r/r 0 ) B inside the Green function (3.4), the solution a µ (3.5) has a priori no reason to verify the gauge condition, as discussed in Sec. 2 of [23].We thus need to verify it and to correct it, if needed, by adding a suitable homogeneous solution.
Using relations among the Legendre functions, it comes On the other hand, the non-linear part of the gauge relation (2.2) yields So the particular solutions (3.5) already satisfy the gauge condition, no further work is needed.Note that taking into account the non-linear part of the gauge relation was crucial here.But, importantly, we cannot generalize this fact, and a gauge-fixing procedure, such as the one described in Sec. 2 of [23] is expected to be vital for higher-order computations.
The first PM iteration of the EM field for the tail interaction thus reads

Correction to the radiative moments
The previous computation (3.9) gives the first PM correction to the MPM expanded EM field.Let us now extract the corresponding corrections to the radiative moments, following the procedure described in Sec.II D 3.
The first step is to develop the fields at future null infinity, i.e. r → ∞ with u = t − r/c constant.Using for instance Eqs.(4.5) and (4.6) of [74], one expands the Legendre functions as where H ℓ = ℓ j=1 1/j is the harmonic number. 7Therefore, the asymptotic expansion of the previous result (3.9) reads .11a) (3.11b) and we recall that A M×B i 0 = 0. Let us now perform a change of coordinate system, to move from the current harmonic one to a radiative one, where the asymptotic expansion of the metric is free of polylogarithmic structure.Such transformation is well known (see e.g. the detailed discussion in [74]) and reads at lowest order Here, b 0 is an (unphysical) constant linked to the time origin in the new coordinate system, and will disappear from all observable quantities.Applying this coordinate change to the asymptotic expansion of the linear metric (2.27), and combining it with (3.11), it turns out that the change of coordinate system simply amounts to the replacement ln(r) → ln(c b 0 ), just as in the GR case.
Therefore, the coordinate change (3.12) removes the logarithmic structure present in the asymptotic expansions of both the metric and the EM field.Can we expect this simultaneous cancellation to hold at higher orders?For a tail with generic multipolarity, M × E L or M × B L , the property of the asymptotic expansion of Q ℓ (3.10) leads us to think that it will be the case.Indeed, in GR, the single coordinate change (3.12) removes the logarithms of all the tails.As for the other interactions, there is good hope that one can perform this simultaneous cancellation by adapting the procedure described in [74] with EM moments.This is obviously more a wishful thinking than a formal proof, and it will have to be verified order by order.
Projecting the asymptotic expansion of the EM field (3.13) on the transverse basis as in Eq. (2.31), one can easily extract the corrections to the radiative moments This result is only the leading order correction to the two lowest-order radiative moments entering the flux.Note that, because of the powers in c in the formulas (2.38) and (2.40), the correction to H i enters at the same PN order than the next-to-leading corrections to E i , namely the "memory" E i × I jk , the "failed-tail" E i × J j and the interactions with gauge moments, not presented here.However, and as discussed in Sec.IV B, none of those corrections enter at the precision required for our numerical applications.
Let us stress that this formalism is valid for any type of compact-supported and post-Newtonian sources, i.e. enjoying weak field and small velocities.This is a first step towards a proper description of environmental and/or physical effects.It can thus be applied to multiple star systems embedded with EM fields, but also supernovae, "mountainous" neutron stars or even systems containing non-ultrarelativistic accretion disks.In Sec.IV, we apply this formalism to an astrophysically-motivated model, encompassing notably the case of binaries made of white dwarfs.

IV. APPLICATION TO REALISTIC SYSTEMS
Let us now apply the wave-generation formalism constructed in the previous sections to a simple realistic model, consisting of two stars bearing magnetic dipoles.The stars are modeled by point particles of mass m A with magnetic dipoles µ i A and vanishing electric dipoles q i A = 0, and we note m = m 1 + m 2 the total mass, ν = m 1 m 2 /m 2 the symmetric mass ratio and δ = (m 1 − m 2 )/m.We will use here the same setup as was done in the study of the conservative sector [33], namely two constant magnetic dipoles, aligned and normal to the orbital plane (as discussed in Sec.IV.A of [33], the motion is indeed planar under those hypothesis) : where l is the normal to the orbital plane and the normalization has been chosen for dimensional purposes.This configuration has been shown to be the state of equilibrium for a binary system of stars with magnetic dipoles [16].Note that, contrary to the study of the conservative sector, we can work "on-shell" and directly enforce (4.1) together with q i A = 0 in the intermediate computations.
As the motion is planar, we can decompose the center-of-mass (CoM) relative velocity of the two bodies as v i = ṙn i + r φλ i , where λ closes the time-dependent orthonormal basis (n, λ, l), and we introduce the convenient shortcuts In the same spirit as [33], we seek to compute the next-to-leading (NLO) PN orders in the radiated fluxes and amplitude modes.

A. Computation of the source moments
Let us briefly comment on the computation of the source moments. 8The techniques have been widely developed in the literature and we let the interested reader refer to the references given hereafter for more details.The first step is to parametrize the gravitational and EM fields in the near-zone in terms of so-called PN potentials.This parametrization was crucial for the derivation of the conservative sector, and thus has been presented in [33].However, in the present problem, we restrict ourselves to the computation of the NLO, and do not take the electric dipoles into account.This implies that we require those parametrizations at a lower PN order than in the previous work [33].The gravitational field is thus parametrized as ) and the EM field, as The potentials {V, V i , W ij , ϕ, χ i } satisfy the wave equations presented in Eqs.(3.7) of [33], and have been computed off-shell to solve the conservative problem.In the present work, we can derive them on-shell, i.e. we can replace the accelerations by the already derived equations of motions, and use directly the ansatz (4.1) for µ i A with q i A = 0.With this parametrization at hand, let us focus on the computation of the EM source moments.From their general expressions (2.24), we start by deriving their PN-expanded sources ῑµ (2.3b).To do so, we perform spatial and temporal indices decompositions and expand using the parametrization of the fields (4.3)-(4.4).This procedure yields integrals with both a compact and a non-compact supported part.The compact sectors result from the compact sources terms (2.5) and (2.6), while the non-compact sector consists of products of the potentials.The techniques used to compute the resulting integrals are based on the Hadamard partie finie regularization and are described in, e.g.[50,[89][90][91].After integration, we end up with expressions of the source moments in terms of the positions and velocities of the particles, expressed in a general frame.The last step is thus to express them in the CoM frame, which is done by applying the method depicted in Sec.III.E. of [33].Even in the CoM frame, the expressions of the source moments are too long to be displayed here, but they constitute only an intermediate result.The relevant (i.e.observable) quantities are indeed the radiative moments, which, at the NLO, are essentially their time derivatives, see Eq. (2.32).

B. Computation of the fluxes and modes
Before using the method described above to compute the source moments, let us investiguate the required precision needed to reach the NLO order for magnetic effects.Let us first concentrate on the EM fluxes (2.38) and (2.40).A rapid inspection of the expression for the source moments (2.24) in our setup reveals that the LO of the electric-type source moments {E L } start at O(c −2 ).Moreover, and as expected, it is easy to realize that B i = A µ i A + O(c −2 ).As µ i A is constant in our model, the radiative moment H i = B (1) i + O(G/c 3 ) begins at O(c −2 ) and, due to the same fact, the associated tail (3.14b) is a 2.5PN relative correction.By similar arguments, one can see that, in our setup, the EM fluxes (2.38) and (2.40) start at O(c −9 ), and that non-linear interactions do not enter the NLO order, i.e. the O(c −11 ) coefficient in the flux.Therefore we can approximate the radiative moments by the (time derived) source moments, as in Eq. (2.32), and we only require As for the gravitational fluxes (2.33), the magnetic effects enter in the radiative moments {U L , V L } at a relative 2PN order, as was already the case in the conservative sector [33].Therefore, the EM sector of those fluxes are also of order O(c −9 ), and we only need the knowledge of the NLO order in U ij and the leading orders of U ijk and V ij .Similarly to the case of the EM fluxes, one can show that there is no EM corrections to the non-linear interactions up to O(c −11 ) in the fluxes, and thus one can safely take All the required (gravitational and EM) radiative moments were derived using the time derivatives of the source moments computed following the procedure depicted previously.Gathering them and applying the formulas (2.33a) and (2.38), the total energy flux 9 is given by where our nomenclature should be transparent.Note that the point-particle contribution F GW,pp is currently known at a relative 4PN order [40], but is only required at NLO in this work.Similarly, working out the formulas (2.33b) and (2.40), we find that the total angular momentum flux is directed along l, and read The explicit expressions of each contribution to NLO are displayed in App.B 1.Then, using the formulas (2.35) and (2.36), one can extract the EM corrections to the gravitational modes, presented in App.B 2. Those results are collected in the ancillary file [49].

C. The case of eccentric orbits
All the results presented in the previous section IV B are given for generic (planar) orbits in terms of (r, ṙ, φ).We will now specify them on quasi-elliptic trajectories in order to extract physical information on the evolution of the orbital parameters.

Quasi-Keplerian parametrization of the orbits
In [33], only quasi-circular orbits were studied.Thus, we first need to derive the quasi-Keplerian parametrization of the trajectories at the required order before proceeding any further.In a post-Newtonian framework, celestial bodies are naturally expected to follow nearly Keplerian orbits.For spinless point particles, the first corrections to the Keplerian orbits were computed in [92], then extended at 2PN [93], 3PN [94], and finally up to the current 4PN precision [95].The aim of this section is to derive a quasi-Keplerian representation of the motion at NLO, by taking the EM corrections into account.
Following the standard procedure, depicted e.g. in [96,97] (notably, all required integrals are displayed in App.A10 of [96]), we can establish the relations between the orbital separation r, the phase φ, the mean anomaly ℓ, the eccentric anomaly u and the true anomaly v.They have the same structure as the pure GR ones at 2PN, namely r = a r 1 − e r cos u , (4.8a) where n = 2π P , with P the relativistic period, and K = 1+k, with k the advance of periastron.The expressions of all parameters entering Eq. (4.8) are given in terms of the invariant energy and angular momentum in App.B 3, and collected in the ancillary file [49].Note that the coefficients translating a departure from a Kepler-like structure, f ⋆ and g ⋆ Eqs.(B10g)-(B10j), are all starting at the same order, namely O(α EM /c 6 ).The relations (4.8) depict entirely the conservative motion of the binary on a quasi-elliptic orbit, and as such, contains all the information we need in the following section.

Evolution of the orbital parameters
Instead of using the gauge-invariant energy and angular momentum, we will describe the evolution of the orbital parameters in terms of the time-eccentricity e t and the orbital frequency-related parameter where Ω = Kn = 2π P (1 + k).From Eqs. (4.8), one can derive r(u), ṙ(u) and φ(u) using chain rules.We finally obtain their expressions in terms of (x, e t , u), which allows to plug them into the fluxes (B1) and (B5) to get the exact instantaneous fluxes F and G in terms of the eccentric anomaly.With such expressions at hand, one can perform their orbit average as and similarly for G .Using the integration formula (8.5) of [98] 1 2π with P n the usual Legendre polynomial, we derived the expressions for the orbital averaged fluxes, presented in App.B 4. Finally, one can use the fluxes together with the explicit expressions (B10a), (B10b), ( B10c) and (B10e) to derive the averaged evolution of the orbital parameters e r , a r and x.Indeed, the flux-balance equations (1.1) allow us to derive the chain rule for, e.g., the time eccentricity where we have denoted the norm of the conserved angular momentum as J = |J i |.Those generalization of the Peters & Mathews formulas [99,100] can be found in App.B 5 as well as in the ancillary file [49].

Gravitational modes
In order to obtain a full waveform, we have also computed the gravitational amplitude modes for quasi-elliptic orbits in terms of the eccentric anomaly u, using the same techniques as those employed to derive the instantaneous expressions of the fluxes.The benefit of such parametrization in place of an expression in terms of the time (or equivalently the mean anomaly ℓ) is that it is exact.In order to express them in terms of ℓ, one needs to invert the relation ℓ(u) (4.8b), which involves an expansion in the time eccentricity e t , see e.g.[101].
The expressions of the instantaneous part of the amplitude modes are too lengthy to be displayed, even in an Appendix, but can be found in the ancillary file [49], and we find full agreement between our point-particle 1PN sector and the results of [97].

D. Quasi-circular orbits : gravitational phase and modes
Finally, let us study the system on quasi-circular orbits.For most of the relevant quantities, one would only need to set the eccentricity to vanish to obtain the results for quasicircular orbits.However, the explicit expression of the phase evolution of the system in terms of the orbital frequency is not provided within the quasi-Keplerian framework.To express the phase, we take ṙ = 0 as it is a higher 2.5PN correction, and set φ = ω where ω is the orbital frequency for a quasi-circular motion.In such a configuration, the fluxes, amplitudes and phase only depend on the gauge invariant PN parameter This variable is naturally the circular limit of the one defined in the eccentric case Eq. ( 4 This relation is a consequence of the first law of binary black hole mechanics [102,103] for systems with constant individual masses, and applies for the conserved quantities of the system.Making use of the flux-balance equations, we checked that it holds for the radiative sector.Thus we can focus solely on the flux-balance equation for the energy to derive the phase.Using the generalized Kepler's law, Eq. (4.5) of [33], on can express the flux in terms of x at NLO We naturally recover (B12a) for e t = 0. From this quantity and the expression of the energy in terms of x, Eq. (4.6a) of [33], one can extract the gravitational phase as At NLO, it comes (4.17)where φ 0 is a constant of integration and the phase has been collected in the ancillary file [49].
Then, we can also express the gravitational modes in circular orbits.For the sake of simplicity, we rescale the observable gravitational modes h ℓm (2.35) by the dominant contribution, and single out the dependence on the gravitational phase φ as Note that the coefficients H ℓm for circular orbits differ from the ones defined in generic orbits, H ℓm (B7) by a factor c 2 x.The relation H ℓ,−m = (−) ℓ H ℓm , where the overbar denotes the complex conjugate, still holds, thus we only present the coefficients with positive m.At the relevant order, the instantaneous modes read and all H ℓ,0 = 0.11

E. Numerical applications
One of the motivations for this work is to estimate the significance of the magnetic contribution to the dissipative sector.For instance, it has been shown in [15] that EM effects on eccentric orbits were not negligible for LISA data analysis, when accumulating a sufficient observation time.Indeed, eccentricity affects each harmonics differently, breaking the degeneracy and making the presence of magnetic dipoles possibly detectable.However, the analysis done in [15] relies solely on the leading order of the Peters & Mathews formulas for point-particles [99,100], neglecting EM effects in the dissipative sector.This approach is not formally consistent in a PN sense.In the present work, we have generalized the Peters & Mathews formulas by including the magnetic dipoles to NLO, see Eqs. (B15) and (B17), which completes the model used in [15].
We can estimate the order of magnitude of the newly computed EM terms for realistic systems, by comparing them to the pure point-particle contribution.To this end, we choose an astrophysically motivated model, that relates the magnitude of the magnetic dipole of a star to its magnetic field at its surface and its radius as [104] As for the numerical values of the masses and magnetic fields of the stars, we took those already used in [11,33], which are typical values for white dwarfs, see Table I.

Parameter Unit
Value for the first star Value for the second star Numerical values taken for the physical parameters of each star, reproduced from [11].m A is the mass of the star A, R A its equatorial radius, and B A the magnitude of the magnetic field at its surface.
Finally, we selected three different configurations for the fundamental observed frequency (which is half the orbital frequency), 10 −1 , 10 −2 and 10 −3 Hz.Those values are motivated by the bandwidth of the LISA detector [2].Numerical estimates of the magnitude of the magnetic contributions to ȧr and ėt are provided in Table II, and it is clear that these corrections, of order at most 10 −6 , are quite small, even in the extreme case e t = 0.5.This suggests that the model used in [15] would yield qualitatively similar conclusions if the correct evolution equations (taking into account the magnetic corrections) were considered.I.
As was the case in the conservative sector [33], the leading EM corrections are smaller than the NLO point-particle ones, but larger than the NNLO point-particle ones, and this for each of the different chosen configurations.This makes the EM contributions roughly comparable to a 1.5PN point-particle term.On the other hand, the EM contributions tend to dominate over the point-particle ones for eccentricities close to 1. But, if our formulas (B15) and (B17) are formally exact for any bound orbit (i.e.0 ≤ e t < 1) within the PN regime, it is a well-known fact that the PN approximation breaks down for high eccentricities.This is due to the fact that, at the periastron of a highly eccentric orbit, the weak-field condition is naturally invalid.Thus, the enhancement of magnetic effects for high eccentricities is spurious.Furthermore, when the two stars are close enough, depicting the magnetic interaction with dipoles is no more possible, as other effects (such mass exchange or common magnetic envelop) enter the picture.
Finally, let us comment on the application to binaries of neutron stars.As discussed in Sec.I, those have much stronger magnetic fields than white dwarfs.However, their magnetic dipoles are significantly smaller, due to the fact that their radius are ∼ 10 3 shorter than those of white dwarfs, see Eq. (4.20).Thus, in our setup of constant and aligned magnetic dipoles, the EM contributions to the decay of orbital parameters for binaries of neutron stars systems are expected to be negligible.

V. SUMMARY AND CONCLUSION
Including electromagnetic effects in the templates used by future gravitational-wave detectors, such as LISA or ET, will be important, notably for an accurate characterization of white dwarfs binaries [10][11][12][13][14][15][16].In this spirit, a proper post-Newtonian treatment of the motion of two stars with EM dipoles has been done in [33], and the usual Noetherian quantities (energy, angular momentum, etc.) have been derived.Pushing further in this direction, the present work tackles the problem of the generation of gravitational-wave by such binaries.
Relying on matched multipolar-post-Minkowskian and post-Newtonian expansions, we first concentrate on the near-zone behavior of the fields, and derive source and gauge moments (2.24), that entirely describe the physics of the matter distribution.From those moments, we iteratively construct non-linear radiative moments that parametrize the values of the fields at future null infinity, (2.30) and (2.31).This asymptotic parametrization allows us to extract the fluxes of energy and angular momentum (2.33), (2.38) and (2.40), as well as the gravitational modes, parametrizing the waveform as detected by an observer at infinity.As a proof of concept, we apply the iteration scheme to explicitly derive the EM tail, i.e. the leading order non-linear correction, of the two lowest-order EM radiative moments (3.14).We emphasize that this formalism is valid for any type of post-Newtonian, compact-supported source (i.e.enjoying weak field and small velocities).It can thus be applied to other astrophysical systems containing EM effects within the PN regime, for instance non-spherical rotating neutron stars, supernovae or even systems containing nonultrarelativistic accretion disks.
Next, we apply this formalism to the concrete case of a binary composed of stars bearing magnetic dipoles.Under the assumption of constant and aligned magnetic dipoles, we derive the fluxes and modes for generic orbits.Enforcing eccentric orbits, we derive the evolution of orbital parameters (eccentricity, semi-major axis and frequency) (B15), (B17) and (B19), and provide the EM corrections to the gravitational modes to the next-to-leading order.For the specific case of quasi-circular orbits, we derive the correction to the gravitational phase.All those results are collected in the ancillary file [49].
Specifying realistic astrophysical configurations, we perform numerical estimates that show the magnitude of the magnetic contribution to the radiated gravitational waves, and dissipative dynamics.Notably, we find that the ratio of the magnetic effects vs. the pure point-particle terms in the secular evolution of the orbital elements are at most ∼ 10 −6 .This implies that the model used in the study [15] would yield the qualitatively similar conclusions if it were to be consistent in a post-Newtonian sense.However, and similarly to [33], we find that the leading magnetic corrections are smaller than a 1PN point-particle effect, but bigger than a 2PN one, making the magnetic effects roughly comparable to a 1.5PN quantity.
If the present work formally tackles the problem of post-Newtonian gravitational wave generation in an Einstein-Maxwell framework, the efforts in building realistic waveform models including EM effects can naturally be improved.In fact, our model assumes constant and aligned magnetic dipoles.If such configuration is indeed an equilibrium state [16], it would be quite interesting to apply our formalism to systems with non-constant and/or nonaligned dipoles.In such configurations, we could have quantitative changes, as the leading order in the EM fluxes (2.38) and (2.40) is vanishing due to the constant dipole hypothesis.This requires to fix the dynamics of the dipoles, which could be done for instance by coupling them to the spins of the stars.This is left for future work.
the field strength reads at leading order where we recall that the transverse projector is ⊥ ij = δ ij − N ij .At the leading order, the Poynting vector where we have not written the cross terms ∝ Q L H K , as their odd parity will eliminate them when performing the angular integrals.The energy flux (A2a) then reads where we have used the fact that the moments {Q L , H L } are STF and the integration formula This generic-gauge derivation agrees with the result (2.38), performed in the transverse gauge, which was expected as the energy flux should be a gauge-invariant quantity.

Angular momentum flux
Let us now turn to the more subtle derivation of the angular momentum flux.As clear from the R 3 appearing in its expression (A2b), and as discussed in lengths in e.g.[26,56,75], its computation requires the knowledge of the sub-leading terms O(1/R 2 ) of the radiative field (3.13).The canonical way of proceeding involves a full analysis of the asymptotic structure at future null infinity, and is plagued by subtleties, see e.g.[75].Instead, we adopt a more pragmatic approach, following the spirit of [26].We start by computing the "source" flux, i.e. the dominant PM order of the flux, by using the linearized metric (2.20).This partial flux will be written in terms of (derivatives) of the source moments E L and B L , hence the name.In order to reconstruct the "full" flux from this source flux, one should in principle add all higher PM orders via tedious computations.Nevertheless, using the argument of [26], we claim that this source flux has the same structure that the "radiative" flux, under the naive replacements E (2.32).This method, although not a formal proof, has the advantage of the simplicity: the linear metric (2.20) is easily expanded at sub-leading order, using Injecting it into the field strength, it comes The flux (A2b) thus becomes (A11) Using the integration formula it comes L ) → (Q L , H L ), this result gives Eqs.(2.40), that we have used in this work.Note however that when applying our formalism to a concrete case in Sec.IV, we have only pushed the accuracy at NLO and the non-linearities did not contribute.In this precise case, the "source" and "full" fluxes are naturally identical.
As discussed in Sec.IV D, the balance equations (1.1) are degenerate in the case of quasicircular orbits.Thus, the phase (4.16) can also be obtained in terms of the balance equation for the angular momentum, i.e.
where G = |G i | and J = |J i |.Implementing this formula, we recover exactly the previous result (4.17).This derivation can be seen as a check regarding the definition of the radiated EM angular momentum.Indeed, the general definition of the GW angular momentum (2.33b) in terms of the radiative moments was, of course, already known and its explicit value contains also magnetic contributions which combine to those of the EM flux to give the expected result.The fact that both derivations agree thus strongly suggests that the adopted definition is correct.As a sidenote, it is also remarkable that the link between the coefficients of F EM and G EM is the same as the proportionality coefficients of those between F GW and G GW .

Appendix B: Lengthy expressions
All the lengthy results presented hereafter are collected in the ancillary file [49].

Expressions of the fluxes for generic orbits
This appendix presents the explicit expressions of the fluxes computed in Sec IV B. We recall that the total energy flux for arbitrary planar motion has been decomposed as The (instantaneous) point particle sector, F GW,pp , can be found up to 3PN in Eq. (5.2) of [98], and is given at 1PN by (B3b) Note that we have written F GW,mag in terms of μ1,2 instead of μ± , for the sake of compactness.Indeed, it comes μ1 μ2 = ν μ2 Similarly to the energy flux, the total angular momentum flux is split according to where the (instantaneous) point-particle sector, G GW,pp , can be found at 3PN e.g. in Eq. (3.2) of [105].At the required order, the contributions read This appendix presents the instantaneous gravitational modes computed in Sec IV B. For the sake of simplicity, we single out the dependence on the gravitational phase φ and rescale the modes h ℓm (2.35) as and further split H ℓm = H pp ℓm + H EM ℓm .Moreover, those coefficients obey H ℓ,−m = (−) ℓ H ℓm , so we only present the ones with a positive m.The point-particle sector, H pp ℓm can be found e.g. in [97], and we don't replicate it here.As for the EM corrections, the instantaneous modes read to the relevant order H EM 44 = 53 24 5 7

Quasi-Keplerian parameters
The parameters entering the quasi-Keplerian representation of the motion (4.8) are given herebelow in terms of the two dimensionless invariants linked to the (unreduced) energy E and norm of the angular momentum where we recall that we deal with bounded trajectories, thus E < 0. While the energy parameter is a 1PN quantity, the second parameter is Newtonian and is such that the eccentricities are simply e ⋆ = √ 1 − j + O(c −2 ).It relates to the usual definition h = J/(Gm 2 ν) by j = −2h 2 E/(mν).At the required order, we obtain The 1PN point-particle sector of those expressions naturally agrees with [92].For the sake of completeness, we also display the expression of x (4.9) at the sought order

Expressions of the fluxes for eccentric orbits
In the spirit of Eqs.(8.8) of [98] and (4.10) of [105], respectively presenting the averaged instantaneous energy and angular momentum fluxes in a system of modified harmonic coordinates, we split at the NLO F = 32 ν 2 c 5 5G x 5 F N + x F 1PN + α EM x 2 F mag,LO + x F mag,NLO , (B12a) with (let us emphasize that the following expressions are exact, in the sense that they have not been truncated in powers of e t ) Naturally, the point-particle sectors at 1PN agree with Eqs (8.9) of [98] and (4.11) of [105].

Time evolution of the orbital parameters
Using the chain rule (4.12) together with Eq. (B10e) and the averaged fluxes (B12), we find the evolution of the time-eccentricity as where and we recognize the Newtonian sector as the "enhancement" factor due to Peters and Mathews [99].
y A = y i A ; we use multi-index notations, i.e.I L = I i1i2...i ℓ ; hats and brackets denote symmetric and tracefree (STF) operators: ÎL = I L = STF L [I L ]; the d'Alembertian operator is defined with respect to the flat Minkowski metric ≡ η µν ∂ µν = ∆ − c −2 ∂ 2 t ; dots and numbers in parenthesis denote time differentiation A (n) = d n A/dt n ; (anti-)symmetrizations are weighted, e.g.A (ij) = (A ij + A ji )/2. 2 Excepted for the numerical applications of Sec.IV, we work with the dimensionless field defined in Sec.III.A of [33], namely A µ ≡ e 2 G/c 3 A dimfull µ , where A dimfull µ denotes the usual electromagnetic field.

TABLE II .
Numerical estimation of the ratios | ȧr EM / ȧr GR | and | ėt EM / ėt GR | for different values of the orbital frequency and [105]s emphasize again that the expressions presented in this section are exact, in the sense that they have not been truncated in powers of e t .Similarly, the semi-major axis decays asda r dt = ν c x 3 A N + x A 1PN + α EM x 2 A mag,LO + x A mag,NLOQuite naturally, the Newtonian factors of Eqs.(B15) and (B17) are nothing but the wellknown Peters' formulas[100].As for the 1PN point-particle sectors, we have a full agreement with Eqs.(6.17) and (6.19) of[105].As for the evolution of the orbital period, it comes X N + x X 1PN + α EM x 2 X mag,LO + x X mag,NLO , + + δ μ+ μ− − ν μ2 + + δ μ+ μ− − ν μ2 − (1 − e 2 t ) 11/2 − 128 − 10768 e 2