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

Quentin Henry<EMAIL_ADDRESS>Max Planck Institute for Gravitational
Physics
(Albert Einstein Institute), D-14476 Potsdam, Germany  François Larrouturou
<EMAIL_ADDRESS>Deutsches Elektronen-Synchrotron DESY, Notkestr.
85, 22607 Hamburg, Germany  Christophe Le Poncin-Lafitte
<EMAIL_ADDRESS>SYRTE, Observatoire de Paris, Université PSL,
CNRS, Sorbonne Université, LNE,
61 avenue de l’Observatoire, 75014 Paris, France

###### Abstract

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.

††preprint: DESY-23-074

##  I Introduction

High accuracy analytic predictions are a milestone of the signal analysis for
gravitational-wave (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 weak-field 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

$\left\langle\frac{\mathrm{d}E}{\mathrm{d}t}\right\rangle=-\langle\mathcal{F}\rangle\qquad\text{and}\qquad\left\langle\frac{\mathrm{d}J^{i}}{\mathrm{d}t}\right\rangle=-\langle\mathcal{G}^{i}\rangle\,,$
(1)

where $\mathcal{F}$ and $\mathcal{G}^{i}$ are respectively the energy and
angular momentum fluxes, as detected by an observer at spatial infinity and
brackets $\langle\cdot\rangle$ 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-to-leading 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
$n$PN 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.2, 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.3; (ii) taking into
account the non-linearities _via_ an iteration scheme and extracting the
_radiative_ moments, as discussed in Sec. II.4; (iii) relating those moments
to the modes and fluxes, as presented in Sec. II.5.

###  II.1 Equations of motion

In the following111The conventions employed throughout this work are as
follows: we use a mostly plus signature, the Minkowski metric being
$\eta_{\mu\nu}=(-,+,+,+)$; greek letters denote spacetime indices
$\mu,\nu,\ldots=(0,1,2,3)$ and latin ones, purely spatial indices
$i,j,\ldots=(1,2,3)$; bold font denotes three-dimensional vectors, _e.g._
$\bm{y}_{A}=y_{A}^{i}$; we use multi-index notations, _i.e._
$\mathrm{I}_{L}=\mathrm{I}_{i_{1}i_{2}\ldots i_{\ell}}$; hats and brackets
denote symmetric and trace-free (STF) operators:
$\hat{\mathrm{I}}_{L}=\mathrm{I}_{\langle
L\rangle}=\underset{L}{\text{STF}}\left[\mathrm{I}_{L}\right]$; the
d’Alembertian operator is defined with respect to the flat Minkowski metric
$\Box\equiv\eta^{\mu\nu}\partial_{\mu\nu}=\Delta-c^{-2}\partial_{t}^{2}$; dots
and numbers in parenthesis denote time differentiation
$A^{(n)}=\mathrm{d}^{n}A/\mathrm{d}t^{n}$; (anti-)symmetrizations are
weighted, _e.g._ $A_{(ij)}=(A_{ij}+A_{ji})/2$., we consider a generic matter
distribution with compact support in an Einstein-Maxwell framework. The metric
perturbation $h^{\mu\nu}\equiv\sqrt{-g}g^{\mu\nu}-\eta^{\mu\nu}$ and EM
field222Excepted for the numerical applications of Sec. IV, we work with the
dimensionless field defined in Sec. III.A of [33], namely
$A_{\mu}\equiv\sqrt{e^{2}G/c^{3}\hbar}\,A^{\text{dimfull}}_{\mu}$, where
$A^{\text{dimfull}}_{\mu}$ denotes the usual electromagnetic field. $A_{\mu}$
are governed by the coupled fields equations [33]

$\Box h^{\mu\nu}=\frac{16\pi
G}{c^{4}}\,\tau^{\mu\nu}\qquad\text{and}\qquad\Box A_{\mu}=-\frac{4\pi
G\alpha_{\text{EM}}}{c^{4}}\,\iota_{\mu}\,,$ (2)

together with the gauge conditions

$\partial_{\nu}h^{\mu\nu}=0\qquad\text{and}\qquad\sqrt{-g}\,\nabla^{\mu}A_{\mu}=\left(\eta^{\mu\nu}+h^{\mu\nu}\right)\partial_{\mu}A_{\nu}=0\,,$
(3)

where $g$ denotes the determinant of $g_{\mu\nu}$, and
$\alpha_{\text{EM}}=\frac{\mu_{0}e^{2}c}{4\pi\hbar}$ is the fine-structure
constant, _i.e._ a dimensionless number. The source terms of the wave
equations (2), $\tau^{\mu\nu}$ and $\iota_{\mu}$, play the role of the Landau-
Lifshitz pseudotensor in GR, and read

$\displaystyle\tau^{\mu\nu}\equiv|g|\,T^{\mu\nu}+\frac{c^{4}}{16\pi
G}\,\Lambda^{\mu\nu}\,,$ (4a) $\displaystyle\iota_{\mu}\equiv
J_{\mu}-\frac{c^{4}}{4\pi G\alpha_{\text{EM}}}\,\Phi_{\mu}\,.$ (4b)

We have separated the sources between the descriptions of the matter
distribution, $\\{T^{\mu\nu},J_{\mu}\\}$, and the non-linear interaction terms
$\\{\Lambda^{\mu\nu},\Phi_{\mu}\\}$. The contribution $\Lambda^{\mu\nu}$ can
be further split between the metric self-interaction term
$\Lambda^{\mu\nu}_{g}$ and the interplay between the metric and the EM field,
dubbed $\Lambda_{\text{EM}}^{\mu\nu}$. The metric self-interaction
$\Lambda^{\mu\nu}_{g}$ is the usual GR one and can found _e.g._ in Eq. (24) of
[23]. The two other contributions read explicitly [33]

$\displaystyle\Lambda^{\mu\nu}_{\text{EM}}$
$\displaystyle=\frac{4|g|}{\alpha_{\text{EM}}}\left(F^{\mu\lambda}F^{\nu}_{\
\lambda}-\frac{g^{\mu\nu}}{4}F_{\alpha\beta}F^{\alpha\beta}\right)\,,$ (5a)
$\displaystyle\Phi_{\mu}$
$\displaystyle=-h^{\alpha\beta}\partial_{\alpha\beta}A_{\mu}-\partial_{\mu}\mathfrak{g}^{\alpha\beta}\,\partial_{\alpha}A_{\beta}-\mathfrak{g}_{\mu\nu}\mathfrak{g}^{\alpha\lambda}\,F_{\alpha\beta}\partial_{\lambda}\mathfrak{g}^{\beta\nu}+\frac{1}{2}\,\mathfrak{g}^{\alpha\beta}\mathfrak{g}_{\lambda\rho}\,F_{\alpha\mu}\partial_{\beta}\mathfrak{g}^{\lambda\rho}\,,$
(5b)

where we have introduced the usual field strength
$F_{\mu\nu}\equiv\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ together with
the “gothic” metric $\mathfrak{g}^{\mu\nu}\equiv\sqrt{-g}g^{\mu\nu}$.

As for the matter sources $\\{T^{\mu\nu},J_{\mu}\\}$, 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
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

$\displaystyle T^{\mu\nu}$
$\displaystyle=\sum_{A}\sqrt{|g|_{A}}\,u_{A}^{0}v_{A}^{\mu}v_{A}^{\nu}\bigg{(}m_{A}-\frac{1}{2c^{2}}\,\left(F_{\lambda\rho}\right)_{A}\mathcal{D}^{\lambda\rho}_{A}\bigg{)}\,\delta_{A}\,,$
(6a) $\displaystyle J_{\mu}$
$\displaystyle=g_{\mu\nu}\sum_{A}\partial_{\lambda}\left(\frac{\mathcal{D}^{\nu\lambda}_{A}}{u_{A}^{0}}\delta_{A}\right)\,,$
(6b)

where $\delta_{A}\equiv\delta[x^{i}-y_{A}^{i}(t)]$ is the 3-dimensional Dirac
delta distribution, locating the matter distribution on the worldline of the
particle $A$, $\bm{y}_{A}(t)$ (thus the sources are indeed compact supported),
$u_{A}^{0}\equiv[-(g_{\mu\nu})_{A}\,v_{A}^{\mu}v_{A}^{\nu}/c^{2}]^{-1/2}$ is
the Lorentz factor, $v_{A}^{\mu}=(c,v_{A}^{i})$ (with $v_{A}^{i}$ 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

$\mathcal{D}^{\mu\nu}_{A}\equiv\begin{pmatrix}0&-c\,q^{i}_{A}\\\
c\,q^{i}_{A}&\varepsilon_{ijk}\mu^{k}_{A}\end{pmatrix}\,,$ (7)

where $\\{q_{A}^{i},\mu_{A}^{i}\\}$ are the electric- and magnetic-type
dipoles and $\varepsilon_{ijk}$ the Levi-Civita symbol.

Last, but not least, note the difference in structure between the two gauge
conditions (3): while harmonic condition on the metric is linear, the $U(1)$
gauge condition on the EM field has a non-linear sector, involving
$h^{\mu\nu}$. This fact will be crucial when explicitly implementing the MPM
iteration scheme in Sec. III.1.

###  II.2 Generic structure of the fields in the exterior zone

The equations of motion (2) are formally solved by the application of a
retarded three-dimensional Green function

$\displaystyle h^{\mu\nu}=\frac{16\pi
G}{c^{4}}\,\Box_{\text{ret}}^{-1}\tau^{\mu\nu}=-\frac{4G}{c^{4}}\int\\!\\!\mathrm{d}^{3}\bm{x}^{\prime}\,\frac{\tau^{\mu\nu}\left(\bm{x}^{\prime},t-|\bm{x}-\bm{x}^{\prime}|/c\right)}{|\bm{x}-\bm{x}^{\prime}|}\,,$
(8a) $\displaystyle A_{\mu}=-\frac{4\pi
G\alpha_{\text{EM}}}{c^{4}}\,\Box_{\text{ret}}^{-1}\iota_{\mu}=\frac{G\alpha_{\text{EM}}}{c^{4}}\int\\!\\!\mathrm{d}^{3}\bm{x}^{\prime}\,\frac{\iota_{\mu}\left(\bm{x}^{\prime},t-|\bm{x}-\bm{x}^{\prime}|/c\right)}{|\bm{x}-\bm{x}^{\prime}|}\,,$
(8b)

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^{\mu\nu}$ and $A_{\mu}$
coincide with their multipolar expansions,
$\mathcal{M}\left(h^{\mu\nu}\right)$ and $\mathcal{M}\left(A_{\mu}\right)$ (as
those are formal series). Moreover, we have by construction
$T^{\mu\nu}=J_{\mu}=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)

$\displaystyle\Box\,\mathcal{M}\left(h^{\mu\nu}\right)=\mathcal{M}\left(\Box
h^{\mu\nu}\right)=\mathcal{M}\left(\Lambda^{\mu\nu}\right)=\Lambda^{\mu\nu}\left[\mathcal{M}\left(h^{\mu\nu}\right),\mathcal{M}\left(A_{\mu}\right)\right]\,,$
(9a) $\displaystyle\Box\,\mathcal{M}\left(A_{\mu}\right)=\mathcal{M}\left(\Box
A_{\mu}\right)=\mathcal{M}\left(\Phi_{\mu}\right)=\Phi_{\mu}\left[\mathcal{M}\left(h^{\mu\nu}\right),\mathcal{M}\left(A_{\mu}\right)\right]\,.$
(9b)

Solving the vacuum equations (9) is simpler than computing the integrals (8).
Nevertheless, the $r\to\infty$ 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 [51, 50], introducing an unphysical (length) scale
$r_{0}$, thus solving (9) as333For high PN computations, it is now customary
to use a dimensional regularization scheme instead of a Hadamard one [52, 53].
Nevertheless, for the low PN order derivation of the present work, the
Hadamard one is largely sufficient and we let the question of dimensional
regularization to future studies.

$\displaystyle\mathcal{M}\left(h^{\mu\nu}\right)=\underset{B=0}{\text{PF}}\,\Box_{\text{ret}}^{-1}\left[\left(\frac{r}{r_{0}}\right)^{B}\,\mathcal{M}\left(\Lambda^{\mu\nu}\right)\right]+h_{\text{hom}}^{\mu\nu}\,,$
(10a)
$\displaystyle\mathcal{M}\left(A_{\mu}\right)=\underset{B=0}{\text{PF}}\,\Box_{\text{ret}}^{-1}\left[\left(\frac{r}{r_{0}}\right)^{B}\,\mathcal{M}\left(\Phi_{\mu}\right)\right]+A^{\text{hom}}_{\mu}\,,$
(10b)

Note that, together with the regularized propagator acting on the non-linear
sources, we have introduced homogeneous solutions, $h_{\text{hom}}^{\mu\nu}$
and $A^{\text{hom}}_{\mu}$. 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

$\displaystyle h_{\text{hom}}^{\mu\nu}=-\frac{4G}{c^{4}}\sum_{\ell\geq
0}\frac{(-)^{\ell}}{\ell!}\partial_{L}\left[\frac{1}{r}\,F_{L}^{\mu\nu}\left(t-r/c\right)\right]\,,$
(11a) $\displaystyle
A^{\text{hom}}_{\mu}=\frac{G\alpha_{\text{EM}}}{c^{4}}\sum_{\ell\geq
0}\frac{(-)^{\ell}}{\ell!}\partial_{L}\left[\frac{1}{r}\,G^{L}_{\mu}\left(t-r/c\right)\right]\,.$
(11b)

where the quantities $F^{\mu\nu}_{L}$ and $G_{\mu}^{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 $\Box
h_{\text{hom}}^{\mu\nu}=\Box A^{\text{hom}}_{\mu}=0$.

The formal separation of the solution (10) bears a physical meaning. On the
one hand, the homogeneous solutions $h_{\text{hom}}^{\mu\nu}$ and
$A^{\text{hom}}_{\mu}$ represent the linear wave emitted by the source. As
such, they encrypt information on the matter distribution, as presented in
Sec. II.3. On the other hand, the regularized propagators acting on the non-
linear sources $\Lambda^{\mu\nu}$ and $\Phi_{\mu}$ describe scattering and
diffusion processes happening during the propagation of the wave towards
spatial infinity. They are discussed in Sec. II.4.

###  II.3 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}^{\mu\nu}$ and $G_{\mu}^{L}$ entering (11) as simple, irreducible STF
_source_ (and _gauge_) moments; (ii) matching those moments to the matter
distributions $\tau^{\mu\nu}$ and $\iota_{\mu}$.

The wave equation for $h^{\mu\nu}$ in terms of $\tau^{\mu\nu}$ (2) 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
$\tau^{\mu\nu}$ in terms of the fields. Therefore, the STF reduction and
matching of the moments in terms of $\tau^{\mu\nu}$ 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^{\mu\nu}_{\text{hom}}=h^{\mu\nu}_{\text{can}}\left[\mathrm{I}_{L},\mathrm{J}_{L}\right]+\partial\varphi^{\mu\nu}\left[\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\right]\,.$
(12)

This generic homogeneous solution is split in two sectors. The _canonical_
linear metric $h^{\mu\nu}_{\text{can}}$ bears the two propagating degrees of
freedom of GR, and is written in terms of two sets of _source_ moments,
$\\{\mathrm{I}_{L},\mathrm{J}_{L}\\}$ [54, 55, 56]. Its expression can be
found _e.g._ in Eq. (36) of [23]. The second sector,
$\partial\varphi^{\mu\nu}\equiv\partial^{\mu}\varphi^{\nu}+\partial^{\nu}\varphi^{\mu}-\partial_{\lambda}\varphi^{\lambda}\,\eta^{\mu\nu}$,
is nothing but the generator of a diffeomorphism transformation, parametrized
by four sets of _gauge_ moments
$\\{\mathrm{W}_{L},\mathrm{X}_{L},\mathrm{Y}_{L},\mathrm{Z}_{L}\\}$, as
explicitly displayed _e.g._ in Eq. (37) of [23]. The matched expressions of
those moments in terms of $\tau^{\mu\nu}$ are to be found in Eqs. (123) and
(125) of [23]. Note that the monopole $\mathrm{I}$ and the dipoles
$\mathrm{I}_{i}$ and $\mathrm{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 $\mathrm{I}_{i}=0$, and denote
$\mathrm{M}=\mathrm{I}$.

####  II.3.1 Irreducible decomposition of the EM field

We need to perform the decomposition of the linear EM field (11b) in
irreducible STF tensors, while ensuring that the gauge equation (3) is
respected. Hereafter, for a quantity $Q_{L}$, we use the shortened notation

$\tilde{Q}_{L}\equiv\frac{1}{r}\,Q_{L}\left(t-r/c\right)\,.$ (13)

Following Eq. (5.3) of [38], let us decompose

$G^{L}_{0}=R_{L}\qquad\text{and}\qquad
G^{L}_{i}=T^{(+)}_{iL}+\varepsilon_{ai\langle i_{\ell}}\,T^{(0)}_{L-1\rangle
a}+\delta_{i\langle i_{\ell}}T^{(-)}_{L-1\rangle}\,,$ (14)

where

$T^{(+)}_{iL}=\underset{iL}{\text{STF}}\,G_{i}^{L}\,,\qquad
T^{(0)}_{L}=\frac{\ell}{\ell+1}\,\varepsilon_{ab\langle
i_{\ell}}\,G_{a}^{L-1\rangle b}\qquad\text{and}\qquad
T^{(-)}_{L-1}=\frac{2\ell-1}{2\ell+1}\,G_{a}^{aL-1}\,.$ (15)

Using the property
$\hat{\partial}_{iL-1}=\partial_{i}\hat{\partial}_{L-1}-\frac{\ell-1}{2\ell-1}\,\delta_{i\langle
i_{\ell-1}}\hat{\partial}_{L-2\rangle}\Delta$ together with the homogeneity
relation $\Delta\tilde{T}^{(-)}_{L}=\ddot{\tilde{T}}^{(-)}_{L}/c^{2}$, it
comes

$\displaystyle
A^{\text{hom}}_{0}=\frac{G\alpha_{\text{EM}}}{c^{4}}\sum_{\ell\geq
0}\frac{(-)^{\ell}}{\ell!}\partial_{L}\tilde{R}_{L}\,,$ (16a) $\displaystyle
A^{\text{hom}}_{i}=\frac{G\alpha_{\text{EM}}}{c^{4}}\sum_{\ell\geq
0}\frac{(-)^{\ell}}{\ell!}\Bigg{[}\frac{\ell\,\partial_{L-1}\ddot{\tilde{T}}^{(-)}_{iL-1}}{(\ell+1)(2\ell+1)c^{2}}-\ell\,\partial_{L-1}\tilde{T}^{(+)}_{iL-1}+\varepsilon_{ijk}\,\partial_{jL-1}\tilde{T}^{(0)}_{kL-1}-\frac{\partial_{iL}\tilde{T}^{(-)}_{L}}{\ell+1}\Bigg{]}\,.$
(16b)

Then, we impose the gauge condition

$\mathfrak{g}^{\mu\nu}\partial_{\mu}A_{\nu}^{\text{hom}}=\eta^{\mu\nu}\partial_{\mu}A_{\nu}^{\text{hom}}+h^{\mu\nu}\partial_{\mu}A_{\nu}^{\text{hom}}=-\frac{1}{c}\,\dot{A}_{0}^{\text{hom}}+\partial_{i}A_{i}^{\text{hom}}+\mathcal{O}(G^{2})\,.$
(17)

As we are investigating the _linear_ solution, we can drop the
$\mathcal{O}(G^{2})$ remainder for now. Nevertheless, this remainder will be
crucial during the iteration procedure, as explicitly shown in Sec. III.1.2.
The linear gauge condition thus translates into the relation

$\eta^{\mu\nu}\partial_{\mu}A_{\nu}^{\text{hom}}=-\frac{G\alpha_{\text{EM}}}{c^{4}}\sum_{\ell\geq
0}\frac{(-)^{\ell}}{\ell!}\partial_{L}\Bigg{[}\frac{\dot{\tilde{R}}_{L}}{c}+\ell\,\tilde{T}^{(+)}_{L}+\frac{1}{2\ell+1}\frac{\ddot{\tilde{T}}^{(-)}_{L}}{c^{2}}\Bigg{]}=0\,.$
(18)

Substituting $T^{(+)}_{L}$ in terms of $R_{L}$ and $T^{(-)}_{L}$ in (16), it
comes

$\displaystyle
A^{\text{hom}}_{0}=\frac{G\alpha_{\text{EM}}}{c^{4}}\sum_{\ell\geq
0}\frac{(-)^{\ell}}{\ell!}\partial_{L}\tilde{R}_{L}\,,$ (19a) $\displaystyle
A^{\text{hom}}_{i}=\frac{G\alpha_{\text{EM}}}{c^{4}}\sum_{\ell\geq
0}\frac{(-)^{\ell}}{\ell!}\Bigg{[}\frac{\partial_{L-1}\dot{\tilde{R}}_{iL-1}}{c}+\frac{\partial_{L-1}\ddot{\tilde{T}}^{(-)}_{iL}}{(\ell+1)\,c^{2}}+\varepsilon_{ijk}\,\partial_{jL-1}\tilde{T}^{(0)}_{kL-1}-\frac{\partial_{iL}\tilde{T}^{(-)}_{L}}{\ell+1}\Bigg{]}\,.$
(19b)

Next, it is common knowledge that an vector field obeying a $U(1)$ symmetry
only propagates two degrees of freedom, whereas the linear solution (19)
contains three sets of moments. Thus, in the same spirit as the splitting of
$h_{\text{hom}}^{\mu\nu}$ (12), $A_{\mu}^{\text{hom}}$ can be rearranged into
a “canonical” EM field, parametrized by two sets of STF _source_ moments,
$\\{\mathrm{E}_{L},\mathrm{B}_{L}\\}$, of respective electric and magnetic
parities, and a $U(1)$ generator, parametrized by a set of STF _gauge_
moments, $\\{\mathrm{K}_{L}\\}$

$A^{\text{hom}}_{\mu}=A_{\mu}^{\text{can}}\left[\mathrm{E}_{L},\mathrm{B}_{L}\right]+\partial_{\mu}\theta\left[\mathrm{K}_{L}\right]\,,$
(20)

where we pose

$\displaystyle
A^{\text{can}}_{0}=-\frac{G\alpha_{\text{EM}}}{c^{3}}\sum_{\ell\geq
0}\frac{(-)^{\ell}}{\ell!}\partial_{L}\tilde{\mathrm{E}}_{L}\,,$ (21a)
$\displaystyle
A^{\text{can}}_{i}=-\frac{G\alpha_{\text{EM}}}{c^{4}}\sum_{\ell\geq
1}\frac{(-)^{\ell}}{\ell!}\bigg{[}\partial_{L-1}\dot{\tilde{\mathrm{E}}}_{iL-1}+\frac{\ell}{\ell+1}\,\varepsilon_{ijk}\,\partial_{jL-1}\,\tilde{\mathrm{B}}_{kL-1}\bigg{]}\,,$
(21b) $\displaystyle\theta=-\frac{G\alpha_{\text{EM}}}{c^{4}}\sum_{\ell\geq
0}\frac{(-)^{\ell}}{\ell!}\partial_{L}\tilde{\mathrm{K}}_{L}\,.$ (21c)

The signs and the powers in $c$ have been fixed to recover the usual results
for constant charges and dipoles [57]. Comparing with Eq. (19), one can read

$\mathrm{E}_{L}=-\frac{R_{L}}{c}-\frac{1}{\ell+1}\,\frac{\dot{T}^{(-)}_{L}}{c^{2}}\,,\qquad\mathrm{B}_{L}=-\frac{\ell+1}{\ell}\,T^{(0)}_{L}\qquad\text{and}\qquad\mathrm{K}_{L}=\frac{1}{\ell+1}\,T^{(-)}_{L}\,.$
(22)

Let us note that for $\ell=0$, the linear gauge condition (17) simply yields
$\dot{\mathrm{E}}=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 $\mathrm{Q}=\mathrm{E}$.

####  II.3.2 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. (9). Therefore, one can apply the method exposed in [38, 23]
to the EM field and derive

$G^{L}_{\mu}=\underset{B=0}{\text{PF}}\int\\!\\!\mathrm{d}^{3}\bm{x}\,\left(\frac{r}{r_{0}}\right)^{B}\,\hat{x}_{L}\,\int_{-1}^{1}\mathrm{d}z\,\delta_{\ell}(z)\,\bar{\iota}_{\mu}\left(\bm{x},t+zr/c\right)\,,$
(23)

where the bar denotes a PN expansion (_i.e._ a formal $c\to\infty$ series),
and we have introduced

$\displaystyle\int_{-1}^{1}\mathrm{d}z\,\delta_{\ell}(z)\,\bar{\iota}_{\mu}\left(\bm{x},t+zr/c\right)$
$\displaystyle\equiv\frac{(2\ell+1)!!}{2^{\ell+1}\,\ell!}\int_{-1}^{1}\mathrm{d}z\,\left(1-z^{2}\right)^{\ell}\,\bar{\iota}_{\mu}\left(\bm{x},t+zr/c\right)$
(24) $\displaystyle=\sum_{k\geq
0}\frac{(2\ell+1)!!}{(2k)!!\,(2\ell+2k+1)!!}\left(\frac{r}{c}\frac{\partial}{\partial
t}\right)^{2k}\bar{\iota}_{\mu}\left(\bm{x},t\right)\,.$

Combining Eqs. (22), (15) and (23), it finally comes

$\displaystyle\mathrm{E}_{L}=-\underset{B=0}{\text{PF}}\int\\!\\!\mathrm{d}^{3}\bm{x}\left(\frac{r}{r_{0}}\right)^{B}\int_{-1}^{1}\\!\\!\mathrm{d}z\left\\{\delta_{\ell}(z)\,\hat{x}_{L}\,\frac{\bar{\iota}_{0}}{c}+\frac{2\ell+1}{(\ell+1)(2\ell+3)\,c^{2}}\,\delta_{\ell+1}(z)\,\hat{x}_{iL}\,\dot{\bar{\iota}}_{i}\right\\}\,,$
(25a) $\displaystyle\mathrm{B}_{L}=\varepsilon_{ab\langle
i_{\ell}}\,\underset{B=0}{\text{PF}}\int\\!\\!\mathrm{d}^{3}\bm{x}\left(\frac{r}{r_{0}}\right)^{B}\int_{-1}^{1}\\!\\!\mathrm{d}z\,\delta_{\ell}(z)\,\hat{x}_{L-1\rangle
a}\,\bar{\iota}_{b}\,,$ (25b)
$\displaystyle\mathrm{K}_{L}=\frac{2\ell+1}{(\ell+1)(2\ell+3)}\,\underset{B=0}{\text{PF}}\int\\!\\!\mathrm{d}^{3}\bm{x}\left(\frac{r}{r_{0}}\right)^{B}\int_{-1}^{1}\\!\\!\mathrm{d}z\,\delta_{\ell+1}(z)\,\hat{x}_{iL}\,\bar{\iota}_{i}\,.$
(25c)

Those expression are in full agreement with the results of [58], that were
derived in the framework of non-relativistic electromagnetism (corresponding
to setting $\Phi_{\mu}=0$).

###  II.4 Multipolar-post-Minkowskian iteration scheme and radiative moments

Once the linear solutions $h_{\text{hom}}^{\mu\nu}$ and $A^{\text{hom}}_{\mu}$
are decomposed in irreducible STF moments, one can turn towards the non-linear
sectors of the solutions (10). In order to do so, we proceed _via_ a post-
Minkowskian expansion

$\mathcal{M}\left(h^{\mu\nu}\right)=\sum_{n\geq
1}G^{n}h_{(n)}^{\mu\nu}\qquad\text{and}\qquad\mathcal{M}\left(A_{\mu}\right)=\sum_{n\geq
1}G^{n}A_{(n)\,\mu}\,,$ (26)

where the linear orders $h_{(1)}^{\mu\nu}$ and $A_{(1)\,\mu}$ are naturally
given by the homogeneous solutions (12) and (20). Note that we combine
multipolar and post-Minkowskian expansions, hence the name of the iteration
scheme.

####  II.4.1 Sketch of the iteration procedure

The procedure is conceptually very simple: the vacuum equations (9) 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_{\mu}$. Nevertheless, verifying the gauge condition on the $n$th PM
order $A_{(n)\,\mu}$ only requires the knowledge of lower-order solutions
$\\{h_{(m<n)}^{\mu\nu},A_{(m<n)\,\mu}\\}$, already required to compute
$A_{(n)\,\mu}$. 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 non-linearities, namely the
_tails_. An interesting fact is that, due to the form of the interaction term
$\Lambda_{\text{EM}}^{\mu\nu}$ (5a), the EM corrections to the gravitational
moments $\\{\mathrm{I}_{L},\mathrm{J}_{L}\\}$ involve at least two EM moments.
For instance, the lowest order EM interactions correcting the mass quadrupole
$\mathrm{I}_{ij}$ are given by $\mathrm{Q}\times\mathrm{E}_{ij}$ and
$\mathrm{E}_{i}\times\mathrm{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 source-to-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 [62, 63, 60]. As the structure of the linear EM field (20)
is very similar to the one of the linear metric (12), 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 $\mathrm{E}_{i}$
(resp. $\mathrm{B}_{i}$) are given by the interactions
$\mathrm{M}\times\mathrm{K}_{i}$, $\mathrm{W}\times\mathrm{E}_{i}$ and
$\mathrm{Q}\times\mathrm{W}_{i}$ (resp. $\mathrm{W}\times\mathrm{B}_{i}$ and
$\mathrm{K}\times\mathrm{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.

####  II.4.2 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\to\infty$, with $u=t-r/c$ constant. Using the relation

$\partial_{L}\bigg{[}\frac{F(t-r/c)}{r}\bigg{]}=\frac{(-)^{\ell}\,n_{L}}{c^{\ell}\,r}\,F^{(\ell)}(t-r/c)+\mathcal{O}\left(\frac{1}{r^{2}}\right)\,.$
(27)

the behavior of the linear solution (21) is easily extracted

$\displaystyle
A_{0}^{\text{lin}}=-\frac{G\alpha_{\text{EM}}}{c^{3}\,r}\sum_{\ell}\frac{1}{c^{\ell}\,\ell!}n_{L}E_{L}^{(\ell)}+\mathcal{O}\left(\frac{1}{r^{2}}\right)\,,$
(28a) $\displaystyle
A_{i}^{\text{lin}}=\frac{G\alpha_{\text{EM}}}{c^{3}\,r}\sum_{\ell}\frac{1}{c^{\ell}\,\ell!}\bigg{[}n_{L-1}E_{iL-1}^{(\ell)}-\frac{\ell}{c(\ell+1)}\,\varepsilon_{iab}n_{aL-1}B_{bL-1}^{(\ell)}\bigg{]}+\mathcal{O}\left(\frac{1}{r^{2}}\right)\,,$
(28b)

where the moments are evaluated at $u$. We recall that the gauge moments
$\\{\mathrm{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]

$A_{\mu}\sim\sum_{\ell,p}\frac{\ln^{p}r}{r}\,n_{L}\,\alpha_{\ell,p}(u)+o\left(\frac{1}{r^{2}}\right)\,,$
(29)

where $p\geq 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.

####  II.4.3 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.

$H^{\mu\nu}\sim\sum_{\ell}\frac{N_{L}}{R}\,H_{L}(U)+\mathcal{O}\left(\frac{1}{R^{2}}\right)\qquad\text{and}\qquad\mathcal{A}_{\mu}\sim\sum_{\ell}\frac{N_{L}}{R}\,A_{L}(U)+\mathcal{O}\left(\frac{1}{R^{2}}\right)\,,$
(30)

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 (53) both the metric
and the EM field takes the form (30).

In such coordinate system, one can project the metric $H^{\mu\nu}$ in a
traceless and transverse (TT) gauge [56]
$H_{0\mu}^{\text{TT}}=N_{i}H_{ij}^{\text{TT}}=H_{ii}^{\text{TT}}=0$ as

$H_{ij}^{\text{TT}}\equiv-\frac{4G}{c^{2}\,R}\,\mathcal{P}^{\text{TT}}_{ijab}\sum_{\ell}\frac{1}{c^{\ell}\,\ell!}\bigg{[}N_{L-2}\mathrm{U}_{abL-2}-\frac{2\ell}{c(\ell+1)}\,N_{cL-2}\varepsilon_{cd(a}\mathrm{V}_{b)dL-2}\bigg{]}+\mathcal{O}\left(\frac{1}{R^{2}}\right)\,,$
(31)

where we introduced the TT projector
$\mathcal{P}^{\text{TT}}_{ijab}=\perp_{i(a}\perp_{b)j}-\frac{1}{2}\perp_{ij}\perp_{ab}$
with $\perp_{ij}=\delta_{ij}-N_{ij}$. The two sets of moments
$\\{\mathrm{U}_{L},\mathrm{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_{0}^{\text{T}}=N_{i}A_{i}^{\text{T}}=0$, as

$A_{i}^{\text{T}}=\frac{G\alpha_{\text{EM}}}{c^{3}\,R}\perp_{ij}\sum_{\ell}\frac{1}{c^{\ell}\,\ell!}\bigg{[}N_{L-1}\mathrm{Q}_{jL-1}-\frac{\ell}{c(\ell+1)}\,\varepsilon_{jab}N_{aL-1}\mathrm{H}_{bL-1}\bigg{]}+\mathcal{O}\left(\frac{1}{R^{2}}\right)\,,$
(32)

where similarly, we have introduced two sets of _radiative_ electric and
magnetic moments, $\\{\mathrm{Q}_{L},\mathrm{H}_{L}\\}$. Looking back at Eq.
(28), it comes as expected

$\mathrm{Q}_{L}=\mathrm{E}_{L}^{(\ell)}+\mathcal{O}(G)\qquad\text{and}\qquad\mathrm{H}_{L}=\mathrm{B}_{L}^{(\ell)}+\mathcal{O}(G)\,,$
(33)

where the remainder encrypts all non-linearities, and notably the tails
computed in Sec. III.

###  II.5 Computation of the modes and fluxes

####  II.5.1 Gravitational sector

The detailed derivation of the energy and angular momentum fluxes associated
with the asymptotic metric (31) can be found in [56], and they read

$\displaystyle\mathcal{F}_{\text{GW}}=\sum_{\ell\geq
2}\frac{G}{c^{2\ell+1}}\frac{(\ell+1)(\ell+2)}{(\ell-1)\ell\,\ell!(2\ell+1)!!}\bigg{[}\mathrm{U}_{L}^{(1)}\mathrm{U}_{L}^{(1)}+\frac{4\ell^{2}}{c^{2}(\ell+1)^{2}}\,\mathrm{V}_{L}^{(1)}\mathrm{V}_{L}^{(1)}\bigg{]}\,,$
(34a)
$\displaystyle\mathcal{G}^{i}_{\text{GW}}=\varepsilon_{iab}\,\sum_{\ell\geq
2}\frac{G}{c^{2\ell+1}}\frac{(\ell+1)(\ell+2)}{(\ell-1)\,\ell!(2\ell+1)!!}\bigg{[}\mathrm{U}_{aL-1}\mathrm{U}_{bL-1}^{(1)}+\frac{4\ell^{2}}{c^{2}(\ell+1)^{2}}\,\mathrm{V}_{aL-1}\mathrm{V}_{bL-1}^{(1)}\bigg{]}\,.$
(34b)

From the asymptotic metric (31), one can also extract the observable
gravitational modes, $h_{\ell\mathrm{m}}$, by projecting $H_{ij}^{\text{TT}}$
onto the basis of polarization $\\{+,\times\\}$ and on spin-weighted spherical
harmonics, $Y_{-2}^{\ell\mathrm{m}}$ (following the conventions of [63, 41])
as

$h_{+}-\mathrm{i}h_{\times}=\sum_{\ell\geq
2}\sum_{\mathrm{m}=-\ell}^{\ell}h_{\ell\mathrm{m}}\,Y_{-2}^{\ell\mathrm{m}}\,.$
(35)

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]

$\displaystyle
h_{\ell\mathrm{m}}=-\frac{2G}{Rc^{\ell+2}\,\ell!}\sqrt{\frac{(\ell+1)(\ell+2)}{\ell(\ell-1)}}\,\alpha_{L}^{\ell\mathrm{m}}\,\mathrm{U}_{L}\,,$
$\displaystyle\text{if}\ \ell+\mathrm{m}\ \text{is even}\,,$ (36a)
$\displaystyle
h_{\ell\mathrm{m}}=-\frac{4G\,\mathrm{i}}{Rc^{\ell+3}\,\ell!}\sqrt{\frac{\ell(\ell+2)}{(\ell-1)(\ell+1)}}\,\alpha_{L}^{\ell\mathrm{m}}\,\mathrm{V}_{L}\,,$
$\displaystyle\text{if}\ \ell+\mathrm{m}\ \text{is odd}\,.$ (36b)

Introducing a fixed orthonormal basis $(\bm{n}_{0},\bm{\lambda}_{0},\bm{l})$
where $\bm{l}$ is normal to the orbital plane, together with
$\bm{\mathfrak{m}}_{0}=(\bm{n}_{0}+\mathrm{i}\bm{\lambda}_{0})/\sqrt{2}$, the
projector $\alpha_{L}^{\ell\mathrm{m}}$ is explicitly given by

$\alpha_{L}^{\ell\mathrm{m}}=\frac{(-)^{\mathrm{m}}\,2^{\frac{2+|\mathrm{m}|}{2}}\sqrt{\pi}\,\ell!}{\sqrt{(2\ell+1)(\ell+\mathrm{m})!\,(\ell-\mathrm{m})!}}\,\overline{\mathfrak{m}}_{0}^{\langle
M}l^{L-M\rangle}\,,$ (37)

where the overbar denotes complex conjugation.

####  II.5.2 Electromagnetic sector

As discussed in App. A.1, the energy loss can be expressed by the angular
average flux of the Landau-Lifshitz pseudotensor $\tau^{\mu\nu}$. At leading
order in $1/R$, the EM sector of this tensor reduces to the usual Poynting
vector $\Pi_{i}=c\,F^{\text{T}}_{ij}F^{\text{T}}_{j0}$. The EM energy flux
then reads

$\mathcal{F}_{\text{EM}}=\frac{c^{4}}{4\pi
G\,\alpha_{\text{EM}}}\lim_{R\to\infty}\left[\int\\!\\!\mathrm{d}^{2}\Omega\,R^{2}\,N_{a}\Pi_{a}\right]=\frac{c^{3}}{4\pi
G\,\alpha_{\text{EM}}}\lim_{R\to\infty}\left[\int\\!\\!\mathrm{d}^{2}\Omega\,R^{2}\,\dot{A}_{i}^{\text{T}}\dot{A}_{i}^{\text{T}}\right]\,,$
(38)

where we used the relation
$\partial_{j}A^{\text{T}}_{i}=-N_{j}\dot{A}^{\text{T}}_{i}/c+\mathcal{O}(1/R)$
together with the transversity. Plugging the asymptotic expression for the EM
field (32), and performing the angular integrals with the help of the formula
(83), it comes

$\mathcal{F}_{\text{EM}}=\sum_{\ell\geq
1}\frac{G\alpha_{\text{EM}}}{c^{3+2\ell}}\frac{\ell+1}{\ell\,\ell!(2\ell+1)!!}\bigg{[}\mathrm{Q}_{L}^{(1)}\mathrm{Q}_{L}^{(1)}+\frac{\ell^{2}}{c^{2}(\ell+1)^{2}}\,\mathrm{H}_{L}^{(1)}\mathrm{H}_{L}^{(1)}\bigg{]}\,.$
(39)

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 (39).

As for the angular momentum flux, it requires more care, just as in the
gravitational case (see _e.g._ [75, 26, 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 (5a) as

$\displaystyle\mathcal{G}_{\text{EM}}^{i}$ $\displaystyle=\frac{c^{4}}{16\pi
G}\,\lim_{R\to\infty}\left[\varepsilon_{iab}\int\\!\\!\mathrm{d}^{2}\Omega\,R^{3}N_{ac}\,\Lambda_{\text{EM}}^{bc}\right]$
(40) $\displaystyle=\frac{c^{4}}{4\pi
G\alpha_{\text{EM}}}\,\lim_{R\to\infty}\left[\varepsilon_{iab}\int\\!\\!\mathrm{d}^{2}\Omega\,R^{3}\,N_{ac}\,\bigg{(}F_{bd}F_{cd}-F_{0b}F_{0c}\bigg{)}\right]\,.$

But, as clear from the $R^{3}$ factor, its evaluation requires the knowledge
of the sub-dominant $\mathcal{O}(1/R^{2})$ order in the asymptotic expansion
of the field (32). 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

$\mathcal{G}^{i}_{\text{EM}}=\varepsilon_{iab}\,\sum_{\ell\geq
1}\frac{G\alpha_{\text{EM}}}{c^{3+2\ell}}\frac{\ell+1}{\ell!(2\ell+1)!!}\bigg{[}\mathrm{Q}_{aL-1}\mathrm{Q}_{bL-1}^{(1)}+\frac{\ell^{2}}{c^{2}(\ell+1)^{2}}\,\mathrm{H}_{aL-1}\mathrm{H}_{bL-1}^{(1)}\bigg{]}\,.$
(41)

Note that, in the case of a rotating electric dipole, we recover the results
of [80, 81].

Together with the gravitational fluxes (34), those fluxes enter the balance
equations for the energy and angular momentum (1), that allow to derive the
secular evolution of orbital quantities, see Sec. IV.3. In the case of quasi-
circular orbits, they are key ingredient to derive the gravitational phase
$\phi=\int\\!\\!\mathrm{d}t\,\omega$ (where $\omega$ is the gravitational
frequency), as presented in Sec. IV.4. Interestingly, the two fluxes (39) and
(41) 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 ($\mathrm{I}_{L}$, $\mathrm{J}_{L}$, $\mathrm{E}_{L}$ or
$\mathrm{B}_{L}$) onto the static structure of the space-time, described by
the constant (ADM) mass $\mathrm{M}$. The gravitational tails
$\mathrm{M}\times\mathrm{I}_{L}$ and $\mathrm{M}\times\mathrm{J}_{L}$ are
known since a long time [82, 37, 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 $\ell=1$ moments
$\mathrm{E}_{i}$ and $\mathrm{B}_{i}$ with the ADM mass $\mathrm{M}$, and will
explicitly apply the MPM iteration scheme pictured in Sec. II.4.1 to those
interactions.444Note that the interaction between the constant charge
$\mathrm{Q}$ and the gravitational dipole $\mathrm{I}_{i}$ formally enters at
the same order than the tail $\mathrm{M}\times\mathrm{E}_{i}$. Nevertheless,
we work in the center-of-mass frame, in which $\mathrm{I}_{i}=0$, so we will
discard this interaction.

###  III.1 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 does555Note that a similar departure from
gauge can also happens in the EFT framework, where the computation has to be
corrected by a removal of _anomalies_ [86]. 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.4.3).

####  III.1.1 Computation of the particular solution

In what follows, we aim to compute the quadratic interactions
$\mathrm{M}\times\mathrm{E}_{i}$ and $\mathrm{M}\times\mathrm{B}_{i}$, by
solving the PM-expanded vacuum equation (9b). The first step is naturally to
construct the non-linear source $\mathcal{M}(\Phi_{\mu})$, by injecting the
MPM expansions of the fields into the expression of $\Phi_{\mu}$. Selecting
only the $\mathrm{M}$ sector of the linear metric (12) and the appropriate
sectors of the linear EM field (21), we need to consider

$\displaystyle h^{00}_{\mathrm{M}}=-\frac{4G\mathrm{M}}{c^{2}\,r}\,,$
$\displaystyle h^{0i}_{\mathrm{M}}=0\,,\qquad\qquad\qquad\quad
h^{ij}_{\mathrm{M}}=0\,.$ (42a) $\displaystyle
A_{0}^{\mathrm{E}_{i}}=\frac{G\alpha_{\text{EM}}}{c^{3}}\,\partial_{i}\tilde{\mathrm{E}}_{i}\,,$
$\displaystyle
A_{i}^{\mathrm{E}_{i}}=\frac{G\alpha_{\text{EM}}}{c^{4}}\,\dot{\tilde{\mathrm{E}}}_{i}\,,$
(42b) $\displaystyle A_{0}^{\mathrm{B}_{i}}=0\,,$ $\displaystyle
A_{i}^{\mathrm{B}_{i}}=\frac{G\alpha_{\text{EM}}}{2c^{4}}\partial_{j}\tilde{\mathrm{B}}_{i|j}\,,$
(42c)

where we recall that we employ the shortcut (13), and have used the notation
of [87], which amounts here to
$\mathrm{B}_{i|j}\equiv\varepsilon_{ijk}\mathrm{B}_{k}$. Those MPM-expansions
are injected into the non-linear source term (5b) developed to quadratic order

$\Phi_{\mu}=-h^{\alpha\beta}\partial_{\alpha\beta}A_{\mu}-\partial_{\mu}h^{\alpha\beta}\partial_{\alpha}A_{\beta}-F_{\alpha\beta}\partial^{\alpha}h^{\beta}_{\mu}-\frac{1}{2}\,F_{\mu\alpha}\partial^{\alpha}h+{\mathcal{O}}(h^{2}A)\,,$
(43)

where indices are operated with the Minkowski metric and $h=h^{\mu}_{\mu}$ is
the trace of $h^{\mu\nu}$. It is easy to see that the equations we need to
solve are of the form

$\Box
A_{\mu}=\sum_{p,\ell}\frac{\hat{n}_{L}}{r^{p}}\,f_{p}^{L}\left(t-\frac{r}{c}\right)\,,$
(44)

and do not involve logarithmic dependencies, so we can apply the formulas
displayed in App. 1 of [59] to compute the _particular_ solutions

$a_{\mu}^{\mathrm{M}\times\mathrm{E}_{i}}=\underset{B=0}{\text{PF}}\Box_{\text{ret}}^{-1}\left[\left(\frac{r}{r_{0}}\right)^{B}\Phi^{\mathrm{M}\times\mathrm{E}_{i}}_{\mu}\right]\qquad\text{and}\qquad
a_{\mu}^{\mathrm{M}\times\mathrm{B}_{i}}=\underset{B=0}{\text{PF}}\Box_{\text{ret}}^{-1}\left[\left(\frac{r}{r_{0}}\right)^{B}\Phi^{\mathrm{M}\times\mathrm{B}_{i}}_{\mu}\right]\,,$
(45)

which explicitly read666We recall that numbers in parenthesis denote time
differentiation, namely $A^{(n)}=\mathrm{d}^{n}A/\mathrm{d}t^{n}$.

$\displaystyle a_{0}^{\mathrm{M}\times\mathrm{E}_{i}}=$
$\displaystyle\,\frac{4G^{2}\mathrm{M}\alpha_{\text{EM}}}{c^{8}}\,n^{i}\,\bigg{[}\int_{1}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}x\,Q_{1}(x)\,\mathrm{E}^{(3)}_{i}\left(t-\frac{xr}{c}\right)+\frac{3c}{4r}\,\mathrm{E}_{i}^{(2)}+\frac{c^{2}}{4r^{2}}\,\mathrm{E}_{i}^{(1)}+\frac{c^{3}}{4r^{3}}\,\mathrm{E}_{i}\bigg{]}\,,$
(46a) $\displaystyle a_{i}^{\mathrm{M}\times\mathrm{E}_{i}}=$
$\displaystyle\,-\frac{4G^{2}\mathrm{M}\alpha_{\text{EM}}}{c^{8}}\bigg{[}\int_{1}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}x\,Q_{0}(x)\,\mathrm{E}^{(3)}_{i}\left(t-\frac{xr}{c}\right)-\frac{c}{3r}\,\mathrm{E}_{i}^{(2)}-\frac{c^{2}}{3r^{2}}\,\mathrm{E}_{i}^{(1)}$
(46b) $\displaystyle\hskip
85.35826pt+\frac{c\,\hat{n}^{ij}}{8r}\,\mathrm{E}_{j}^{(2)}+\frac{c^{2}\,\hat{n}^{ij}}{8r^{2}}\,\mathrm{E}_{j}^{(1)}\bigg{]}\,,$
$\displaystyle a_{0}^{\mathrm{M}\times\mathrm{B}_{i}}=$ $\displaystyle\,0\,,$
(46c) $\displaystyle a_{i}^{\mathrm{M}\times\mathrm{B}_{i}}=$
$\displaystyle\,\frac{2G^{2}\mathrm{M}\alpha_{\text{EM}}}{c^{8}}\frac{n^{j}}{r}\bigg{[}\int_{1}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}x\,Q_{1}(x)\,\mathrm{B}_{i|j}^{(3)}\left(t-\frac{xr}{c}\right)+\frac{5}{8}\,\mathrm{B}^{(2)}_{i|j}-\frac{c}{8r}\,\mathrm{B}^{(1)}_{i|j}-\frac{c^{2}}{8r^{2}}\,\mathrm{B}_{i|j}\bigg{]}\,.$
(46d)

where the Legendre functions $Q_{\ell}$ are defined with a branch cut on
$(-\infty,1)$ as

$Q_{\ell}(x)=\frac{1}{2}\,P_{\ell}(x)\ln\bigg{(}\frac{x+1}{x-1}\bigg{)}-\sum_{j=1}^{\ell}\frac{1}{j}\,P_{\ell-j}(x)P_{j-1}(x)\,,$
(47)

with $P_{\ell}$ the Legendre polynomials.

####  III.1.2 Implementation of the gauge condition

Due to the presence of the regularization kernel $(r/r_{0})^{B}$ inside the
Green function (45), the solution $a_{\mu}$ (46) 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

$\displaystyle\partial_{i}a_{i}^{\mathrm{M}\times\mathrm{E}_{i}}-\frac{1}{c}\,\partial_{t}a_{0}^{\mathrm{M}\times\mathrm{E}_{i}}$
$\displaystyle=-\frac{4G^{2}\mathrm{M}\alpha_{\text{EM}}}{c^{7}}\frac{n^{i}}{r^{2}}\,\bigg{[}\mathrm{E}_{i}^{(2)}+\frac{c}{r}\,\mathrm{E}_{i}^{(1)}\bigg{]}\,,$
(48)
$\displaystyle\partial_{i}a_{i}^{\mathrm{M}\times\mathrm{B}_{i}}-\frac{1}{c}\,\partial_{t}a_{0}^{\mathrm{M}\times\mathrm{B}_{i}}$
$\displaystyle=0\,.$

On the other hand, the non-linear part of the gauge relation (3) yields

$h^{\mu\nu}_{\mathrm{M}}\,\partial_{\mu}A_{\nu}^{\mathrm{E}_{i}}=\frac{1}{c}\,h^{00}_{\mathrm{M}}\,\partial_{t}A_{0}^{\mathrm{E}_{i}}=\frac{4G^{2}\mathrm{M}\alpha_{\text{EM}}}{c^{7}}\frac{n^{i}}{r^{2}}\,\bigg{[}\mathrm{E}_{i}^{(2)}+\frac{c}{r}\,\mathrm{E}_{i}^{(1)}\bigg{]}\,,\qquad
h^{\mu\nu}_{\mathrm{M}}\,\partial_{\mu}A_{\nu}^{\mathrm{B}_{i}}=0\,.$ (49)

So the particular solutions (46) 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

$\displaystyle A_{0}^{\mathrm{M}\times\mathrm{E}_{i}}=$
$\displaystyle\,\frac{4G^{2}\mathrm{M}\alpha_{\text{EM}}}{c^{8}}\,n^{i}\,\bigg{[}\int_{1}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}x\,Q_{1}(x)\,\mathrm{E}^{(3)}_{i}\left(t-\frac{xr}{c}\right)+\frac{3c}{4r}\,\mathrm{E}_{i}^{(2)}+\frac{c^{2}}{4r^{2}}\,\mathrm{E}_{i}^{(1)}+\frac{c^{3}}{4r^{3}}\,\mathrm{E}_{i}\bigg{]}\,,$
(50a) $\displaystyle A_{i}^{\mathrm{M}\times\mathrm{E}_{i}}=$
$\displaystyle\,-\frac{4G^{2}\mathrm{M}\alpha_{\text{EM}}}{c^{8}}\bigg{[}\int_{1}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}x\,Q_{0}(x)\,\mathrm{E}^{(3)}_{i}\left(t-\frac{xr}{c}\right)-\frac{c}{3r}\,\mathrm{E}_{i}^{(2)}-\frac{c^{2}}{3r^{2}}\,\mathrm{E}_{i}^{(1)}$
(50b) $\displaystyle\hskip
85.35826pt+\frac{c\,\hat{n}^{ij}}{8r}\,\mathrm{E}_{j}^{(2)}+\frac{c^{2}\,\hat{n}^{ij}}{8r^{2}}\,\mathrm{E}_{j}^{(1)}\bigg{]}\,,$
$\displaystyle A_{0}^{\mathrm{M}\times\mathrm{B}_{i}}=$ $\displaystyle\,0\,,$
(50c) $\displaystyle A_{i}^{\mathrm{M}\times\mathrm{B}_{i}}=$
$\displaystyle\,\frac{2G^{2}\mathrm{M}\alpha_{\text{EM}}}{c^{8}}\frac{n^{j}}{r}\bigg{[}\int_{1}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}x\,Q_{1}(x)\,\mathrm{B}_{i|j}^{(3)}\left(t-\frac{xr}{c}\right)+\frac{5}{8}\,\mathrm{B}^{(2)}_{i|j}-\frac{c}{8r}\,\mathrm{B}^{(1)}_{i|j}-\frac{c^{2}}{8r^{2}}\,\mathrm{B}_{i|j}\bigg{]}\,.$
(50d)

###  III.2 Correction to the radiative moments

The previous computation (50) 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.4.3.

The first step is to develop the fields at future null infinity, _i.e._
$r\to\infty$ with $u=t-r/c$ constant. Using for instance Eqs. (4.5) and (4.6)
of [74], one expands the Legendre functions as

$\int_{1}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}x\,Q_{\ell}(x)F\bigg{(}t-\frac{xr}{c}\bigg{)}=-\frac{c}{2r}\int_{0}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}\tau\bigg{[}\ln\left(\frac{c\tau}{2r}\right)+2\mathcal{H}_{\ell}\bigg{]}F\big{(}u-\tau\big{)}+\mathcal{O}\left(\frac{\ln
r}{r^{2}}\right)\,,$ (51)

where $\mathcal{H}_{\ell}=\sum_{j=1}^{\ell}1/j$ is the harmonic number.777Note
that the coefficient in front of the logarithm does not depend on the value of
the index of the Legendre function. This fact will be crucial when applying
the coordinate change (53) to remove the $\ln(r)$ dependency. Therefore, the
asymptotic expansion of the previous result (50) reads

$\displaystyle
A_{0}^{\mathrm{M}\times\mathrm{E}_{i}}=-\frac{2G^{2}\mathrm{M}\alpha_{\text{EM}}\,n_{i}}{c^{7}\,r}\int_{0}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}\tau\bigg{[}\ln\left(\frac{c\tau}{2r}\right)+\frac{1}{2}\bigg{]}\mathrm{E}_{i}^{(3)}\big{(}u-\tau\big{)}+\mathcal{O}\left(\frac{\ln
r}{r^{2}}\right)\,,$ (52a) $\displaystyle
A_{i}^{\mathrm{M}\times\mathrm{E}_{i}}=\frac{2G^{2}\mathrm{M}\alpha_{\text{EM}}}{c^{7}\,r}\Bigg{\\{}\int_{0}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}\tau\bigg{[}\ln\left(\frac{c\tau}{2r}\right)+\frac{3}{4}\bigg{]}\mathrm{E}_{i}^{(3)}\big{(}u-\tau\big{)}-\frac{n^{i}}{4}n_{j}\mathrm{E}^{(2)}_{j}(u)\Bigg{\\}}+\mathcal{O}\left(\frac{\ln
r}{r^{2}}\right)\,,$ (52b) $\displaystyle
A_{i}^{\mathrm{M}\times\mathrm{B}_{i}}=-\frac{G^{2}\mathrm{M}\alpha_{\text{EM}}\,n_{j}}{c^{8}\,r}\int_{0}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}\tau\bigg{[}\ln\left(\frac{c\tau}{2r}\right)+\frac{3}{4}\bigg{]}\mathrm{B}_{i|j}^{(3)}\big{(}u-\tau\big{)}+\mathcal{O}\left(\frac{\ln
r}{r^{2}}\right)\,,$ (52c)

and we recall that $A_{0}^{\mathrm{M}\times\mathrm{B}_{i}}=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

$T=t-\frac{2G\mathrm{M}}{c^{3}}\,\ln\left(\frac{r}{c\,b_{0}}\right)+\mathcal{O}\left(G^{2}\right)\,,\qquad
X^{i}=x^{i}+\mathcal{O}\left(G^{2}\right)\,.$ (53)

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 (28),
and combining it with (52), it turns out that the change of coordinate system
simply amounts to the replacement $\ln(r)\to\ln(c\,b_{0})$, just as in the GR
case.

$\displaystyle\mathcal{A}_{0}^{\mathrm{M}\times\mathrm{E}_{i}}=-\frac{2G^{2}\mathrm{M}\alpha_{\text{EM}}\,N_{i}}{c^{7}\,R}\int_{0}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}\tau\bigg{[}\ln\left(\frac{\tau}{2b_{0}}\right)+\frac{1}{2}\bigg{]}\mathrm{E}_{i}^{(3)}\big{(}U-\tau\big{)}+\mathcal{O}\left(\frac{\ln
R}{R^{2}}\right)\,,$ (54a)
$\displaystyle\mathcal{A}_{i}^{\mathrm{M}\times\mathrm{E}_{i}}=\frac{2G^{2}\mathrm{M}\alpha_{\text{EM}}}{c^{7}\,R}\Bigg{\\{}\int_{0}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}\tau\bigg{[}\ln\left(\frac{\tau}{2b_{0}}\right)+\frac{3}{4}\bigg{]}\mathrm{E}_{i}^{(3)}\big{(}U-\tau\big{)}-\frac{N^{i}}{4}\,N_{j}\mathrm{E}_{j}^{(2)}\big{(}U\big{)}\Bigg{\\}}+\mathcal{O}\left(\frac{\ln
R}{R^{2}}\right)\,,$ (54b)
$\displaystyle\mathcal{A}_{i}^{\mathrm{M}\times\mathrm{B}_{i}}=-\frac{G^{2}\mathrm{M}\alpha_{\text{EM}}\,N_{j}}{c^{8}\,R}\int_{0}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}\tau\bigg{[}\ln\left(\frac{\tau}{2b_{0}}\right)+\frac{3}{4}\bigg{]}\mathrm{B}_{i|j}^{(3)}\big{(}U-\tau\big{)}+\mathcal{O}\left(\frac{\ln
R}{R^{2}}\right)\,.$ (54c)

Therefore, the coordinate change (53) 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, $\mathrm{M}\times\mathrm{E}_{L}$ or
$\mathrm{M}\times\mathrm{B}_{L}$, the property of the asymptotic expansion of
$Q_{\ell}$ (51) leads us to think that it will be the case. Indeed, in GR, the
single coordinate change (53) 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 (54) on the transverse
basis as in Eq. (32), one can easily extract the corrections to the radiative
moments

$\displaystyle\mathrm{Q}_{i}$
$\displaystyle=\mathrm{E}_{i}^{(1)}+\frac{2G\mathrm{M}}{c^{3}}\int_{0}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}\tau\bigg{[}\ln\left(\frac{\tau}{2b_{0}}\right)+\frac{3}{4}\bigg{]}\mathrm{E}_{i}^{(3)}\big{(}U-\tau\big{)}+\mathcal{O}\left(\frac{G}{c^{5}}\right)+{\mathcal{O}}\left(G^{2}\right)\,,$
(55a) $\displaystyle\mathrm{H}_{i}$
$\displaystyle=\mathrm{B}_{i}^{(1)}+\frac{2G\mathrm{M}}{c^{3}}\int_{0}^{\infty}\\!\\!\\!\\!\\!\mathrm{d}\tau\bigg{[}\ln\left(\frac{\tau}{2b_{0}}\right)+\frac{3}{4}\bigg{]}\mathrm{B}_{i}^{(3)}\big{(}U-\tau\big{)}+\mathcal{O}\left(\frac{G}{c^{5}}\right)+{\mathcal{O}}\left(G^{2}\right)\,.$
(55b)

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 (39) and (41), the correction to $\mathrm{H}_{i}$ enters at
the same PN order than the next-to-leading corrections to $\mathrm{E}_{i}$,
namely the “memory” $\mathrm{E}_{i}\times\mathrm{I}_{jk}$, the “failed-tail”
$\mathrm{E}_{i}\times\mathrm{J}_{j}$ and the interactions with gauge moments,
not presented here. However, and as discussed in Sec. IV.2, 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 $\mu_{A}^{i}$ and vanishing electric dipoles $q_{A}^{i}=0$,
and we note $m=m_{1}+m_{2}$ the total mass, $\nu=m_{1}m_{2}/m^{2}$ the
symmetric mass ratio and $\delta=(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) :

$\mu_{A}^{i}=\frac{\tilde{\mu}_{A}}{Gm^{2}}\,l^{i}\,,$ (56)

where $\bm{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 (56) together with $q_{A}^{i}=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}=\dot{r}n^{i}+r\dot{\phi}\lambda^{i}$,
where $\bm{\lambda}$ closes the time-dependent orthonormal basis
$(\bm{n},\bm{\lambda},\bm{l})$, and we introduce the convenient shortcuts

$\tilde{\mu}_{+}\equiv\tilde{\mu}_{1}+\tilde{\mu}_{2}\qquad\text{and}\qquad\tilde{\mu}_{-}\equiv\frac{m}{m_{2}}\,\tilde{\mu}_{2}-\frac{m}{m_{1}}\,\tilde{\mu}_{1}\,.$
(57)

In the same spirit as [33], we seek to compute the next-to-leading (NLO) PN
orders in the radiated fluxes and amplitude modes.

###  IV.1 Computation of the source moments

Let us briefly comment on the computation of the source moments.888All the
computations presented hereafter were performed with the use of the xAct
library from the Mathematica software [88]. The 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

$\displaystyle h^{00}$
$\displaystyle=-\frac{4V}{c^{2}}-\frac{4}{c^{4}}\left[2V^{2}+\frac{W_{kk}}{2}\right]+\frac{2\varphi^{2}}{c^{6}}+\mathcal{O}(c^{-8})\,,$
(58a) $\displaystyle h^{0i}$
$\displaystyle=-\frac{4V_{i}}{c^{3}}+\mathcal{O}(c^{-5})\,,$ (58b)
$\displaystyle h^{ij}$
$\displaystyle=-\frac{4}{c^{4}}\left(W_{ij}-\frac{\delta_{ij}}{2}\,W_{kk}\right)+\mathcal{O}(c^{-6})\,,$
(58c)

and the EM field, as

$A_{0}=\frac{\varphi}{c^{3}}-\frac{V\varphi}{c^{5}}+\mathcal{O}(c^{-7})\,,\qquad
A_{i}=\frac{\chi_{i}}{c^{4}}+\mathcal{O}(c^{-6})\,.$ (59)

The potentials $\\{V,V_{i},W_{ij},\varphi,\chi_{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 (56) for $\mu_{A}^{i}$ with
$q_{A}^{i}=0$.

With this parametrization at hand, let us focus on the computation of the EM
source moments. From their general expressions (25), we start by deriving
their PN-expanded sources $\bar{\iota}_{\mu}$ (4b). To do so, we perform
spatial and temporal indices decompositions and expand using the
parametrization of the fields (58)–(59). This procedure yields integrals with
both a compact and a non-compact supported part. The compact sectors result
from the compact sources terms (6) and (7), 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. [89, 50, 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.
(33).

###  IV.2 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 (39) and (41). A rapid
inspection of the expression for the source moments (25) in our setup reveals
that the LO of the electric-type source moments $\\{\mathrm{E}_{L}\\}$ start
at ${\mathcal{O}}(c^{-2})$. Moreover, and as expected, it is easy to realize
that $\mathrm{B}_{i}=\sum_{A}\mu_{A}^{i}+{\mathcal{O}}(c^{-2})$. As
$\mu_{A}^{i}$ is constant in our model, the radiative moment
$\mathrm{H}_{i}=\mathrm{B}_{i}^{(1)}+{\mathcal{O}}(G/c^{3})$ begins at
${\mathcal{O}}(c^{-2})$ and, due to the same fact, the associated tail (55b)
is a 2.5PN relative correction. By similar arguments, one can see that, in our
setup, the EM fluxes (39) and (41) start at ${\mathcal{O}}(c^{-9})$, and that
non-linear interactions do not enter the NLO order, _i.e._ the
${\mathcal{O}}(c^{-11})$ coefficient in the flux. Therefore we can approximate
the radiative moments by the (time derived) source moments, as in Eq. (33),
and we only require
$\\{\mathrm{E}_{i},\mathrm{E}_{ij},\mathrm{B}_{i},\mathrm{B}_{ij},\mathrm{B}_{ijk}\\}$
to
$\\{{\mathcal{O}}(c^{-4}),{\mathcal{O}}(c^{-2}),{\mathcal{O}}(c^{-2}),{\mathcal{O}}(c^{-2}),{\mathcal{O}}(c^{0})\\}$.

As for the gravitational fluxes (34), the magnetic effects enter in the
radiative moments $\\{\mathrm{U}_{L},\mathrm{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 ${\mathcal{O}}(c^{-9})$, and we
only need the knowledge of the NLO order in $\mathrm{U}_{ij}$ and the leading
orders of $\mathrm{U}_{ijk}$ and $\mathrm{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 ${\mathcal{O}}(c^{-11})$ in the fluxes, and thus one can
safely take

$\mathrm{U}_{L}=\mathrm{U}_{L}^{\text{pp}}+\frac{\alpha_{\text{EM}}}{c^{4}}\,\mathrm{U}_{L}^{\text{EM}}=\mathrm{I}_{L}^{(\ell)}+\mathcal{O}(G)\qquad\text{and}\qquad\mathrm{V}_{L}=\mathrm{V}_{L}^{\text{pp}}+\frac{\alpha_{\text{EM}}}{c^{4}}\,\mathrm{V}_{L}^{\text{EM}}=\mathrm{J}_{L}^{(\ell)}+\mathcal{O}(G)\,.$
(60)

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 (34a) and (39),
the total energy flux999The Noetherian fluxes given in Eqs. (3.14) of [33] are
vanishing in our configuration. is given by

$\mathcal{F}=\mathcal{F}_{\text{GW}}+\mathcal{F}_{\text{EM}}=\mathcal{F}_{\text{GW,pp}}+\mathcal{F}_{\text{GW,mag}}+\mathcal{F}_{\text{EM}}\,,$
(61)

where our nomenclature should be transparent. Note that the point-particle
contribution $\mathcal{F}_{\text{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 (34b) and (41), we find that the total angular momentum flux is
directed along $\bm{l}$, and read

$\mathcal{G}^{i}=\mathcal{G}^{i}_{\text{GW}}+\mathcal{G}^{i}_{\text{EM}}=\bigg{[}\mathcal{G}_{\text{GW,pp}}+\mathcal{G}_{\text{GW,mag}}+\mathcal{G}_{\text{EM}}\bigg{]}l^{i}\,.$
(62)

The explicit expressions of each contribution to NLO are displayed in App.
B.1. Then, using the formulas (36) and (37), one can extract the EM
corrections to the gravitational modes, presented in App. B.2. Those results
are collected in the ancillary file [49].

###  IV.3 The case of eccentric orbits

All the results presented in the previous section IV.2 are given for generic
(planar) orbits in terms of $(r,\dot{r},\dot{\phi})$. We will now specify them
on quasi-elliptic trajectories in order to extract physical information on the
evolution of the orbital parameters.

####  IV.3.1 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. A101010Note a missing $e_{r}$ in
front of the $\sin u$ of (A.7) of [96]. of [96]), we can establish the
relations between the orbital separation $r$, the phase $\phi$, the mean
anomaly $\ell$, the eccentric anomaly $u$ and the true anomaly $v$. They have
the same structure as the pure GR ones at 2PN, namely

$\displaystyle r=a_{r}\big{(}1-e_{r}\,\cos u\big{)}\,,$ (63a)
$\displaystyle\ell\equiv n\big{(}t-t_{0}\big{)}=u-e_{t}\,\sin
u+f_{v-u}\big{(}v-u\big{)}+f_{v}\,\sin v\,,$ (63b)
$\displaystyle\frac{\phi-\phi_{0}}{K}=v+g_{2v}\,\sin(2v)+g_{3v}\,\sin(3v)\,,$
(63c) $\displaystyle
v=2\arctan\bigg{[}\sqrt{\frac{1+e_{\phi}}{1-e_{\phi}}}\,\tan\frac{u}{2}\bigg{]}\,,$
(63d)

where $n=\frac{2\pi}{P}$, with $P$ the relativistic period, and $K=1+k$, with
$k$ the advance of periastron. The expressions of all parameters entering Eq.
(63) 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_{\star}$ and
$g_{\star}$ Eqs. (99g)–(99j), are all starting at the same order, namely
$\mathcal{O}(\alpha_{\text{EM}}/c^{6})$. The relations (63) 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.

####  IV.3.2 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

$x=\left(\frac{Gm\Omega}{c^{3}}\right)^{2/3}\,,$ (64)

where $\Omega=Kn=\frac{2\pi}{P}(1+k)$. From Eqs. (63), one can derive $r(u)$,
$\dot{r}(u)$ and $\dot{\phi}(u)$ using chain rules. We finally obtain their
expressions in terms of $(x,e_{t},u)$, which allows to plug them into the
fluxes (90) and (94) to get the exact instantaneous fluxes $\mathcal{F}$ and
$\mathcal{G}$ in terms of the eccentric anomaly. With such expressions at
hand, one can perform their orbit average as

$\left\langle\mathcal{F}\right\rangle\equiv\frac{1}{P}\int_{0}^{P}\\!\\!\mathrm{d}t\,\mathcal{F}(t)=\frac{1}{2\pi}\int_{0}^{2\pi}\\!\\!\mathrm{d}u\,\frac{\mathrm{d}\ell}{\mathrm{d}u}\,\mathcal{F}(u)\,,$
(65)

and similarly for $\left\langle\mathcal{G}\right\rangle$. Using the
integration formula (8.5) of [98]

$\frac{1}{2\pi}\int_{0}^{2\pi}\\!\\!\frac{\mathrm{d}u}{(1-e_{t}\,\cos
u)^{n}}=\frac{1}{(1-e_{t}^{2})^{n/2}}\,P_{n-1}\bigg{(}\frac{1}{\sqrt{1-e_{t}^{2}}}\bigg{)}\,,$
(66)

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 (99a), (99b), (99c) and (99) to
derive the averaged evolution of the orbital parameters $e_{r}$, $a_{r}$ and
$x$. Indeed, the flux-balance equations (1) allow us to derive the chain rule
for, _e.g._ , the time eccentricity

$\left\langle\frac{\mathrm{d}e_{t}}{\mathrm{d}t}\right\rangle=-\frac{\partial
e_{t}}{\partial E}\left\langle\mathcal{F}\right\rangle-\frac{\partial
e_{t}}{\partial J}\left\langle\mathcal{G}\right\rangle\,,$ (67)

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

####  IV.3.3 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 $\ell$) is
that it is exact. In order to express them in terms of $\ell$, one needs to
invert the relation $\ell(u)$ (63b), 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].

###  IV.4 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 quasi-circular 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
$\dot{r}=0$ as it is a higher 2.5PN correction, and set $\dot{\phi}=\omega$
where $\omega$ 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

$x=\left(\frac{Gm\omega}{c^{3}}\right)^{2/3}\,.$ (68)

This variable is naturally the circular limit of the one defined in the
eccentric case Eq. (64). The two flux-balance equations (1) are now degenerate
and the fluxes are indirectly linked by the “thermodynamic” relation

$\frac{\mathrm{d}E}{\mathrm{d}\omega}=\omega\,\frac{\mathrm{d}J}{\mathrm{d}\omega}\,.$
(69)

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

$\displaystyle\mathcal{F}=\frac{32c^{5}\nu^{2}x^{5}}{5G}\Bigg{\\{}$
$\displaystyle
1-\bigg{(}\frac{1247}{336}+\frac{35\nu}{12}\bigg{)}x-4\alpha_{\text{EM}}\Bigg{[}\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\bigg{(}\frac{1}{96}+\nu\bigg{)}\tilde{\mu}_{-}^{2}\Bigg{]}x^{2}$
(70)
$\displaystyle-\alpha_{\text{EM}}\Bigg{[}\bigg{(}\frac{755}{168}-\frac{169\nu}{6}\bigg{)}\tilde{\mu}_{+}^{2}+\bigg{(}\frac{685}{168}-\frac{169\nu}{6}\bigg{)}\delta\,\tilde{\mu}_{-}\tilde{\mu}_{+}$
$\displaystyle\qquad\qquad-\bigg{(}\frac{19}{32}+\frac{3193\nu}{1008}-\frac{169\nu^{2}}{6}\bigg{)}\tilde{\mu}_{-}^{2}\Bigg{]}x^{3}\Bigg{\\}}\,.$

We naturally recover (101a) 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

$\phi=\int\\!\\!\mathrm{d}t\,\omega=-\frac{c^{3}}{Gm}\int\\!\\!\mathrm{d}x\,\frac{x^{3/2}}{\mathcal{F}}\,\frac{\mathrm{d}E}{\mathrm{d}x}\,.$
(71)

At NLO, it comes

$\displaystyle\phi=\phi_{0}-\frac{1}{32\nu x^{5/2}}\Bigg{\\{}$ $\displaystyle
1+\bigg{(}\frac{3715}{1008}+\frac{55\nu}{12}\bigg{)}x+50\,\alpha_{\text{EM}}\Bigg{[}\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\bigg{(}\frac{1}{240}+\nu\bigg{)}\tilde{\mu}_{-}^{2}\Bigg{]}x^{2}$
(72)
$\displaystyle-\alpha_{\text{EM}}\Bigg{[}\bigg{(}\frac{5865}{14}+\frac{280\nu}{3}\bigg{)}\tilde{\mu}_{+}^{2}+\bigg{(}\frac{35015}{84}+\frac{280\nu}{3}\bigg{)}\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}$
$\displaystyle\qquad\qquad-\bigg{(}\frac{16945}{4032}+\frac{416785\nu}{1008}+\frac{280\nu^{2}}{3}\bigg{)}\tilde{\mu}_{-}^{2}\Bigg{]}x^{3}\Bigg{\\}}\,.$

where $\phi_{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_{\ell\mathrm{m}}$ (36) by the dominant contribution, and single out the
dependence on the gravitational phase $\phi$ as

$h_{\ell\mathrm{m}}=\frac{8Gm\nu
x}{c^{2}R}\,\sqrt{\frac{\pi}{5}}\,H_{\ell\mathrm{m}}\,\mathrm{e}^{-\mathrm{i}\mathrm{m}\phi}\,.$
(73)

Note that the coefficients $H_{\ell\mathrm{m}}$ for circular orbits differ
from the ones defined in generic orbits, $\mathcal{H}_{\ell\mathrm{m}}$ (96)
by a factor $c^{2}x$. The relation
$H_{\ell,-\mathrm{m}}=(-)^{\ell}\,\overline{H}_{\ell\mathrm{m}}$, where the
overbar denotes the complex conjugate, still holds, thus we only present the
coefficients with positive $\mathrm{m}$. At the relevant order, the
instantaneous modes read

$\displaystyle
H_{22}=1-\bigg{(}\frac{107}{42}-\frac{55\nu}{42}\bigg{)}x-\frac{2\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{\nu}x^{2}\bigg{[}1-\bigg{(}\frac{1}{21}-\frac{73\nu}{42}\bigg{)}x\bigg{]}\,,$
(74a) $\displaystyle
H_{21}=\frac{\mathrm{i}}{3}\,\delta\,x^{1/2}\bigg{[}1-3\,\frac{\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{\nu}\,x^{2}\bigg{]}\,,$
(74b) $\displaystyle
H_{33}=-\frac{3\mathrm{i}}{4}\sqrt{\frac{15}{14}}\,\delta\,x^{1/2}\bigg{[}1-3\,\frac{\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{\nu}\,x^{2}\bigg{]}\,,$
(74c) $\displaystyle
H_{32}=\frac{1}{3}\sqrt{\frac{5}{7}}\big{(}1-3\nu\big{)}\,x\,\bigg{[}1-4\,\frac{\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{\nu}\,x^{2}\bigg{]}\,,$
(74d) $\displaystyle
H_{31}=\frac{\mathrm{i}}{12\sqrt{14}}\,\delta\,x^{1/2}\bigg{[}1-3\,\frac{\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{\nu}\,x^{2}\bigg{]}\,,$
(74e) $\displaystyle
H_{44}=-\frac{8}{9}\sqrt{\frac{5}{7}}\big{(}1-3\nu\big{)}\,x\,\bigg{[}1-4\,\frac{\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{\nu}\,x^{2}\bigg{]}\,,$
(74f) $\displaystyle
H_{42}=\frac{\sqrt{5}}{63}\big{(}1-3\nu\big{)}\,x\,\bigg{[}1-4\,\frac{\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{\nu}\,x^{2}\bigg{]}\,,$
(74g)

and all $H_{\ell,0}=0$.111111This is due to the fact that, in this paper, we
only consider instantaneous contributions to the modes. In fact, there are
contributions of the memory to the $(\ell,0)$ modes for even $\ell$ starting
at Newtonian order. In the case of eccentric orbits, these memory
contributions also contribute to oscillating modes and could contribute to the
magnetic sector of $h_{22}$ at NLO. We neglect these contributions throughout
this paper.

###  IV.5 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. (104) and (106), 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]

$|\mu_{A}^{i}|=\frac{2\pi}{\mu_{0}}\,B_{A}R_{A}^{3}\,.$ (75)

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 1.

Parameter | Unit | Value for the first star | Value for the second star
---|---|---|---
$m_{A}$ | $M_{\odot}$ | $1.2$ | $0.3$
$R_{A}$ | km | $6.0\cdot 10^{3}$ | $15\cdot 10^{3}$
$B_{A}$ | G | $1.0\cdot 10^{9}$ | $1.0\cdot 10^{9}$
Table 1: 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 $\langle\dot{a}_{r}\rangle$ and $\langle\dot{e}_{t}\rangle$
are provided in Table 2, 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.

| $|\langle\dot{a}_{r}\rangle^{\text{EM}}/\langle\dot{a}_{r}\rangle^{\text{GR}}|$ | $|\langle\dot{e}_{t}\rangle^{\text{EM}}/\langle\dot{e}_{t}\rangle^{\text{GR}}|$
---|---|---
$e_{t}$ | Config. 1 | Config. 2 | Config. 3 | Config. 1 | Config. 2 | Config. 3
0.1 | $7.0\cdot 10^{-7}$ | $3.3\cdot 10^{-8}$ | $1.5\cdot 10^{-9}$ | $6.8\cdot 10^{-7}$ | $3.2\cdot 10^{-8}$ | $1.5\cdot 10^{-9}$
0.3 | $9.8\cdot 10^{-7}$ | $4.6\cdot 10^{-8}$ | $2.1\cdot 10^{-9}$ | $8.9\cdot 10^{-7}$ | $4.1\cdot 10^{-8}$ | $1.9\cdot 10^{-9}$
0.5 | $1.8\cdot 10^{-6}$ | $8.2\cdot 10^{-8}$ | $3.8\cdot 10^{-9}$ | $1.5\cdot 10^{-6}$ | $7.1\cdot 10^{-8}$ | $3.3\cdot 10^{-9}$
Table 2: Numerical estimation of the ratios
$|\langle\dot{a}_{r}\rangle^{\text{EM}}/\langle\dot{a}_{r}\rangle^{\text{GR}}|$
and
$|\langle\dot{e}_{t}\rangle^{\text{EM}}/\langle\dot{e}_{t}\rangle^{\text{GR}}|$
for different values of the orbital frequency and eccentricities.
Configuration 1 corresponds to $\Omega=5\cdot 10^{-2}$Hz, configuration 2, to
$\Omega=5\cdot 10^{-3}$Hz and configuration 3, to $\Omega=5\cdot 10^{-4}$Hz,
all of those based on the values of Table 1.

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 (104) and
(106) are formally exact for any bound orbit (_i.e._ $0\leq 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 $\sim 10^{3}$ shorter than those of white dwarfs,
see Eq. (75). 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 (25), 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, (31) and (32). This asymptotic parametrization allows us to extract
the fluxes of energy and angular momentum (34), (39) and (41), 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 (55). 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
non-ultrarelativistic 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) (104), (106) and (108), 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 $\sim 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 non-
aligned dipoles. In such configurations, we could have quantitative changes,
as the leading order in the EM fluxes (39) and (41) 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.

###### Acknowledgements.

It is a pleasure to thank C. Dlapa, G. Faye, M. Riva and L. Speri for
interesting discussions. F.L. received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
program (grant agreement No 817791). C.LPL. acknowledges the financial support
of Centre National d’Études Spatiales (CNES) for LISA.

##  Appendix A Electro-magnetic fluxes

This appendix discusses the construction of the EM fluxes (39) and (41) in
terms of the radiative moments.

###  A.1 Generic expressions

Flux balance equations within the PN-MPM formalism have been treated in, e.g.
[26]. In this reference, the radiated fluxes are defined using the angular
integral of the Landau-Lifshitz pseudotensor in the radiative zone. In the
present paper, this tensor is now composed of the pure GR and the EM parts.
However, as discussed in Sec. II.4, both fields bear the same asymptotic
multipolar structure. Thus, the same defintion holds when generalizing to EM
effects. This allows us to directly use Eqs. (4.4) of [26] by defining the
total energy and angular momentum fluxes as

$\displaystyle\mathcal{F}\equiv\lim_{R\to\infty}c\,R^{2}\int\\!\\!\mathrm{d}^{2}\Omega\,N_{i}\tau^{0i}\,,$
(76a)
$\displaystyle\mathcal{G}^{i}\equiv\lim_{R\to\infty}R^{3}\,\varepsilon_{iab}\int\\!\\!\mathrm{d}^{2}\Omega\,N_{ac}\tau^{bc}\,.$
(76b)

The compact sector of $\tau^{\mu\nu}$ (4a) cannot contribute to the surface
integrals, so we are left with the contribution of
$\Lambda^{\mu\nu}=\Lambda_{g}^{\mu\nu}+\Lambda_{\text{EM}}^{\mu\nu}$. The
contribution of $\Lambda_{g}^{\mu\nu}$ has been treated in [56] and gives the
usual gravitational fluxes (34). Focusing on $\Lambda_{\text{EM}}^{\mu\nu}$,
we define the “EM” part of the fluxes as

$\displaystyle\mathcal{F}_{\text{EM}}\equiv\lim_{R\to\infty}\left[\frac{c^{5}\,R^{2}}{16\pi
G}\int\\!\\!\mathrm{d}^{2}\Omega\,N_{i}\Lambda_{\text{EM}}^{0i}\right]\mathrel{\mathop{=}\limits_{R\to\infty}}\frac{c^{4}\,R^{2}}{4\pi
G\alpha_{\text{EM}}}\int\\!\\!\mathrm{d}^{2}\Omega\,N_{i}\,\Pi_{i}\,,$ (77a)
$\displaystyle\mathcal{G}_{\text{EM}}^{i}\equiv\lim_{R\to\infty}\left[\frac{c^{4}\,R^{3}}{16\pi
G}\,\varepsilon_{iab}\int\\!\\!\mathrm{d}^{2}\Omega\,N_{ac}\,\Lambda_{\text{EM}}^{bc}\right]\mathrel{\mathop{=}\limits_{R\to\infty}}\frac{c^{4}\,R^{3}}{4\pi
G\alpha_{\text{EM}}}\,\varepsilon_{iab}\int\\!\\!\mathrm{d}^{2}\Omega\,N_{ac}\,\bigg{(}F_{bd}F_{cd}-F_{0b}F_{0c}\bigg{)}\,,$
(77b)

where we recall $\Pi_{i}=c\,F_{ij}F_{j0}$ is the Poynting vector.

###  A.2 Gauge-independent derivation of the energy flux

In Sec. II.5.2, the energy flux has been derived from the asymptotic expansion
of $A_{\mu}$ in a transverse gauge. In order to emphasize that our result (39)
is indeed gauge-independent, let us re-derive it using a generic (_i.e._ not
transverse) form of the EM field

$\displaystyle A_{0}=-\frac{G\alpha_{\text{EM}}}{c^{3}\,R}\sum_{\ell\geq
1}\frac{1}{c^{\ell}\,\ell!}N_{L}\mathrm{Q}_{L}+\mathcal{O}\left(\frac{1}{R^{2}}\right)\,,$
(78a) $\displaystyle A_{i}=\frac{G\alpha_{\text{EM}}}{c^{3}\,R}\sum_{\ell\geq
1}\frac{1}{c^{\ell}\,\ell!}\bigg{[}N_{L-1}\mathrm{Q}_{iL-1}-\frac{\ell}{c(\ell+1)}\,\varepsilon_{iab}N_{aL-1}\mathrm{H}_{bL-1}\bigg{]}+\mathcal{O}\left(\frac{1}{R^{2}}\right)\,,$
(78b)

where the moments are naturally evaluated in the (radiative) retarded time
$T-R/c$. After some massaging, using notably the relations

$\displaystyle\partial_{i}\Bigg{(}\frac{N_{L}}{R}\,\mathrm{Q}_{L}\Bigg{)}=-\frac{N_{iL}}{R\,c}\,\dot{\mathrm{Q}}_{L}+\mathcal{O}\left(\frac{1}{R^{2}}\right)\,,$
(79a)
$\displaystyle\partial_{i}\Bigg{(}\varepsilon_{jab}\,\frac{N_{aL-1}}{R}\,\mathrm{H}_{bL-1}\Bigg{)}=-\varepsilon_{jab}\,\frac{N_{iaL-1}}{R\,c}\,\dot{\mathrm{H}}_{bL-1}+\mathcal{O}\left(\frac{1}{R^{2}}\right)\,,$
(79b)

the field strength reads at leading order

$\displaystyle F_{0i}=\frac{G\alpha_{\text{EM}}}{c^{4}\,R}\sum_{\ell\geq
1}\frac{1}{c^{\ell}\,\ell!}\bigg{[}\perp_{ij}N_{L-1}\dot{\mathrm{Q}}_{jL-1}-\frac{\ell}{c(\ell+1)}\,\varepsilon_{iab}N_{aL-1}\dot{\mathrm{H}}_{bL-1}\bigg{]}\,,$
(80a) $\displaystyle F_{ij}=\frac{G\alpha_{\text{EM}}}{c^{4}\,R}\sum_{\ell\geq
1}\frac{1}{c^{\ell}\,\ell!}\bigg{[}N_{jL-1}\dot{\mathrm{Q}}_{iL-1}-N_{iL-1}\dot{\mathrm{Q}}_{jL-1}-\frac{\ell\,N_{aL-1}}{c(\ell+1)}\,\big{(}\varepsilon_{iab}N_{j}-\varepsilon_{jab}N_{i}\big{)}\dot{\mathrm{H}}_{bL-1}\bigg{]}\,,$
(80b)

where we recall that the transverse projector is
$\perp_{ij}=\delta_{ij}-N_{ij}$. At the leading order, the Poynting vector
$\Pi_{i}=c\,F_{ij}F_{j0}$ becomes

$\Pi_{i}=\frac{G^{2}\alpha_{\text{EM}}^{2}}{c^{7}\,R^{2}}\sum_{\ell,k}\frac{N_{i}N_{L-1}N_{K-1}}{c^{\ell+k}\,\ell!\,k!}\,\perp_{ab}\,\Bigg{[}\dot{\mathrm{Q}}_{aL-1}\dot{\mathrm{Q}}_{bK-1}+\frac{\ell
k}{c^{2}(\ell+1)(k+1)}\,\dot{\mathrm{H}}_{aL-1}\dot{\mathrm{H}}_{bK-1}\bigg{]}\,,$
(81)

where we have not written the cross terms
$\propto\mathrm{Q}_{L}\mathrm{H}_{K}$, as their odd parity will eliminate them
when performing the angular integrals. The energy flux (77a) then reads

$\displaystyle\mathcal{F}_{\text{EM}}$
$\displaystyle=\frac{G\alpha_{\text{EM}}}{4\pi
c^{3}}\sum_{\ell,k}\frac{1}{c^{\ell+k}\,\ell!\,k!}\,\Bigg{[}\dot{\mathrm{Q}}_{aL-1}\dot{\mathrm{Q}}_{bK-1}+\frac{\ell
k}{c^{2}(\ell+1)(k+1)}\,\dot{\mathrm{H}}_{aL-1}\dot{\mathrm{H}}_{bK-1}\bigg{]}\int\\!\\!\mathrm{d}^{2}\Omega\perp_{ab}N_{L-1}N_{K-1}$
(82) $\displaystyle=\sum_{\ell\geq
1}\frac{G\alpha_{\text{EM}}}{c^{3+2\ell}}\frac{\ell+1}{\ell\,\ell!(2\ell+1)!!}\bigg{[}\dot{\mathrm{Q}}_{L}\dot{\mathrm{Q}}_{L}+\frac{\ell^{2}}{c^{2}(\ell+1)^{2}}\,\dot{\mathrm{H}}_{L}\dot{\mathrm{H}}_{L}\bigg{]}$

where we have used the fact that the moments
$\\{\mathrm{Q}_{L},\mathrm{H}_{L}\\}$ are STF and the integration formula

$\dot{\mathrm{Q}}_{aL-1}\dot{\mathrm{Q}}_{bK-1}\int\\!\\!\mathrm{d}^{2}\Omega\,\perp_{ab}N_{L-1}N_{K-1}=\frac{4\pi\,(\ell+1)!}{\ell\,(2\ell+1)!!}\,\delta_{k\ell}\,\dot{\mathrm{Q}}_{L}\dot{\mathrm{Q}}_{L}\,.$
(83)

This generic-gauge derivation agrees with the result (39), performed in the
transverse gauge, which was expected as the energy flux should be a gauge-
invariant quantity.

###  A.3 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 (77b), and as discussed in
lengths in _e.g._ [56, 75, 26], its computation requires the knowledge of the
sub-leading terms $\mathcal{O}(1/R^{2})$ of the radiative field (54). 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 (21). This partial flux will be
written in terms of (derivatives) of the source moments $\mathrm{E}_{L}$ and
$\mathrm{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
$\big{(}\mathrm{E}_{L}^{(\ell)},\mathrm{B}_{L}^{(\ell)})\to(\mathrm{Q}_{L},\mathrm{H}_{L})$,
see Eq. (33). This method, although not a formal proof, has the advantage of
the simplicity: the linear metric (21) is easily expanded at sub-leading
order, using

$\hat{\partial}_{L}\bigg{[}\frac{F(t-r/c)}{r}\bigg{]}=\frac{(-)^{\ell}\hat{n}_{L}}{r\,c^{\ell}}\bigg{[}F^{(\ell)}(t-r/c)+\frac{\ell(\ell+1)}{2}\,\frac{c}{r}\,F^{(\ell-1)}(t-r/c)\bigg{]}+\mathcal{O}\left(\frac{1}{r^{3}}\right)\,.$
(84)

Injecting it into the field strength, it comes

$\displaystyle F_{0i}=\frac{G\alpha_{\text{EM}}}{c^{4}\,r}\sum_{\ell\geq
1}\frac{1}{c^{\ell}\,\ell!}\bigg{[}\perp_{ij}n_{L-1}\mathrm{E}^{(\ell+1)}_{jL-1}-\frac{\ell}{c(\ell+1)}\,\varepsilon_{iab}n_{aL-1}\mathrm{B}^{(\ell+1)}_{bL-1}\bigg{]}$
$\displaystyle\qquad\quad+\frac{G\alpha_{\text{EM}}}{c^{3}\,r^{2}}\sum_{\ell\geq
1}\frac{1}{c^{\ell}\,\ell!}\bigg{[}\frac{\ell(\ell+1)}{2}\,n_{L-1}\mathrm{E}^{(\ell)}_{iL-1}-\frac{(\ell+1)(\ell+2)}{2}\,n_{iL}\mathrm{E}^{(\ell)}_{L}-\frac{\ell^{2}}{2c}\,\varepsilon_{iab}n_{aL-1}\mathrm{B}^{(\ell)}_{bL-1}\bigg{]}\,,$
(85a) $\displaystyle F_{ij}=\frac{G\alpha_{\text{EM}}}{c^{4}\,r}\sum_{\ell\geq
1}\frac{1}{c^{\ell}\,\ell!}\bigg{[}n_{jL-1}\mathrm{E}^{(\ell+1)}_{iL-1}-n_{iL-1}\mathrm{E}^{(\ell+1)}_{jL-1}-\frac{\ell\,n_{aL-1}}{c(\ell+1)}\,\big{(}\varepsilon_{iab}n_{j}-\varepsilon_{jab}n_{i}\big{)}\mathrm{B}^{(\ell+1)}_{bL-1}\bigg{]}$
$\displaystyle\qquad\quad+\frac{G\alpha_{\text{EM}}}{c^{3}\,r^{2}}\sum_{\ell\geq
1}\frac{1}{c^{\ell}\,\ell!}\bigg{[}\frac{\ell(\ell+1)}{2}\bigg{(}n_{jL-1}\mathrm{E}^{(\ell)}_{iL-1}-n_{iL-1}\mathrm{E}^{(\ell)}_{jL-1}\bigg{)}$
$\displaystyle\hskip
142.26378pt+\bigg{(}\frac{2\ell}{(\ell+1)c}\,\varepsilon_{ija}+\frac{\ell(\ell+2)}{2c}(\varepsilon_{iab}n_{jb}-\varepsilon_{jab}n_{ib})\bigg{)}n_{L-1}\mathrm{B}^{(\ell)}_{aL-1}$
$\displaystyle\hskip
142.26378pt+\frac{\ell(\ell-1)}{(\ell+1)c}\bigg{(}\varepsilon_{iab}\mathrm{B}^{(\ell)}_{jbL-2}-\varepsilon_{jab}\mathrm{B}^{(\ell)}_{ibL-2}\bigg{)}n_{aL-2}\bigg{]}\,.$
(85b)

The flux (77b) thus becomes

$\mathcal{G}_{\text{EM}}^{i}=\frac{G\alpha_{\text{EM}}}{4\pi
c^{3}}\sum_{\ell,k}\frac{\ell+1}{c^{\ell+k}\,\ell!\,k!}\varepsilon_{iab}\Bigg{[}\mathrm{E}_{L}^{(\ell)}\mathrm{E}_{bK-1}^{(k+1)}+\frac{k\,\ell}{c^{2}(k+1)(\ell+1)}\mathrm{B}_{L}^{(\ell)}\mathrm{B}_{bK-1}^{(k+1)}\Bigg{]}\int\\!\\!\mathrm{d}^{2}\Omega\,n_{aK-1}n_{L}\,.$
(86)

Using the integration formula

$\varepsilon_{iab}\,\mathrm{E}^{(\ell)}_{L}\mathrm{E}^{(k+1)}_{bK-1}\,\int\\!\\!\mathrm{d}^{2}\Omega\,n_{aK-1}n_{L}=\frac{4\pi\,\ell!}{(2\ell+1)!!}\,\delta_{k\ell}\,\varepsilon_{iab}\,\mathrm{E}^{(\ell)}_{aL-1}\mathrm{E}^{(\ell+1)}_{bL-1}\,,$
(87)

it comes

$\mathcal{G}_{\text{EM}}^{i}=\varepsilon_{iab}\,\sum_{\ell\geq
1}\frac{G\alpha_{\text{EM}}}{c^{3+2\ell}}\frac{\ell+1}{\ell!\,(2\ell+1)!!}\Bigg{[}\mathrm{E}^{(\ell)}_{aL-1}\mathrm{E}^{(\ell+1)}_{bL-1}+\frac{\ell^{2}}{c^{2}(\ell+1)^{2}}\,\mathrm{B}^{(\ell)}_{aL-1}\mathrm{B}^{(\ell+1)}_{bL-1}\bigg{]}\,.$
(88)

Under the replacement
$\big{(}\mathrm{E}_{L}^{(\ell)},\mathrm{B}_{L}^{(\ell)})\to(\mathrm{Q}_{L},\mathrm{H}_{L})$,
this result gives Eqs. (41), 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.4, the balance equations (1) are degenerate in the
case of quasi-circular orbits. Thus, the phase (71) can also be obtained in
terms of the balance equation for the angular momentum, i.e.

$\phi=-\frac{c^{3}}{Gm}\int\\!\\!\mathrm{d}x\,\frac{x^{3/2}}{\mathcal{G}}\,\frac{\mathrm{d}J}{\mathrm{d}x}\,,$
(89)

where $\mathcal{G}=|\mathcal{G}^{i}|$ and $J=|J^{i}|$. Implementing this
formula, we recover exactly the previous result (72). 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 (34b) 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 $\mathcal{F}_{\text{EM}}$
and $\mathcal{G}_{\text{EM}}$ is the same as the proportionality coefficients
of those between $\mathcal{F}_{\text{GW}}$ and $\mathcal{G}_{\text{GW}}$.

##  Appendix B Lengthy expressions

All the lengthy results presented hereafter are collected in the ancillary
file [49].

###  B.1 Expressions of the fluxes for generic orbits

This appendix presents the explicit expressions of the fluxes computed in Sec
IV.2. We recall that the total energy flux for arbitrary planar motion has
been decomposed as

$\mathcal{F}=\mathcal{F}_{\text{GW,pp}}+\mathcal{F}_{\text{GW,mag}}+\mathcal{F}_{\text{EM}}\,.$
(90)

The (instantaneous) point particle sector, $\mathcal{F}_{\text{GW,pp}}$, can
be found up to 3PN in Eq. (5.2) of [98], and is given at 1PN by

$\displaystyle\mathcal{F}_{\text{GW,pp}}$
$\displaystyle=\frac{32G^{3}m^{4}\nu^{2}}{5c^{5}r^{4}}\Bigg{\\{}r^{2}\dot{\phi}^{2}+\frac{\dot{r}^{2}}{12}$
(91) $\displaystyle\hskip
79.6678pt+\frac{1}{c^{2}}\Bigg{[}\bigg{(}\frac{785}{336}-\frac{71\nu}{28}\bigg{)}r^{4}\dot{\phi}^{4}-\bigg{(}\frac{117}{28}-\frac{45\nu}{14}\bigg{)}r^{2}\dot{r}^{2}\dot{\phi}^{2}-\bigg{(}\frac{8}{21}-\frac{3\nu}{14}\bigg{)}\dot{r}^{4}$
$\displaystyle\hskip
113.81102pt-\bigg{(}\frac{170}{21}-\frac{10\nu}{21}\bigg{)}Gmr\dot{\phi}^{2}+\bigg{(}\frac{9}{14}+\frac{5\nu}{42}\bigg{)}\frac{Gm\dot{r}^{2}}{r}+\frac{1-4\nu}{21}\frac{G^{2}m^{2}}{r^{2}}\Bigg{]}\Bigg{\\}}\,.$

As for the EM contributions, they read at NLO

$\displaystyle\mathcal{F}_{\text{GW,mag}}$
$\displaystyle=-\frac{192G^{5}m^{6}\nu\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{5c^{9}r^{6}}\Bigg{\\{}r^{2}\dot{\phi}^{2}-\frac{\dot{r}^{2}}{12}$
$\displaystyle\hskip
142.26378pt+\frac{1}{c^{2}}\Bigg{[}\bigg{(}\frac{277}{112}-\frac{569\nu}{84}\bigg{)}r^{4}\dot{\phi}^{4}-\bigg{(}\frac{1159}{168}-\frac{1373\nu}{84}\bigg{)}r^{2}\dot{r}^{2}\dot{\phi}^{2}$
$\displaystyle\hskip
170.71652pt+\bigg{(}\frac{17}{42}-\frac{37\nu}{42}\bigg{)}\dot{r}^{4}-\bigg{(}\frac{1331}{168}-\frac{27\nu}{7}\bigg{)}Gmr\dot{\phi}^{2}$
$\displaystyle\hskip
170.71652pt+\bigg{(}\frac{17}{18}-\frac{109\nu}{84}\bigg{)}\frac{Gm\dot{r}^{2}}{r}+\frac{1-4\nu}{14}\frac{G^{2}m^{2}}{r^{2}}\Bigg{]}\Bigg{\\}}\,,$
(92a) $\displaystyle\mathcal{F}_{\text{EM}}$
$\displaystyle=\frac{4G^{5}m^{6}\nu^{2}\alpha_{\text{EM}}\tilde{\mu}_{-}^{2}}{15c^{9}r^{6}}\Bigg{\\{}r^{2}\dot{\phi}^{2}+4\dot{r}^{2}$
$\displaystyle\hskip
113.81102pt+\frac{1}{c^{2}}\Bigg{[}\bigg{(}\frac{397}{28}-\frac{465\nu}{14}\bigg{)}r^{4}\dot{\phi}^{4}+\bigg{(}\frac{1177}{14}-\frac{906\nu}{7}\bigg{)}r^{2}\dot{r}^{2}\dot{\phi}^{2}$
$\displaystyle\hskip
142.26378pt-\bigg{(}\frac{128}{7}-\frac{324\nu}{7}\bigg{)}\dot{r}^{4}-\bigg{(}\frac{71}{14}+\frac{45\nu}{7}\bigg{)}Gmr\dot{\phi}^{2}$
$\displaystyle\hskip
142.26378pt-\bigg{(}\frac{488}{7}-\frac{648\nu}{7}\bigg{)}\frac{Gm\dot{r}^{2}}{r}+8\frac{1-4\nu}{7}\frac{G^{2}m^{2}}{r^{2}}\Bigg{]}\Bigg{\\}}$
$\displaystyle\quad+\frac{64G^{5}m^{6}\nu^{2}\delta\alpha_{\text{EM}}\tilde{\mu}_{+}\tilde{\mu}_{-}}{105c^{11}r^{6}}\Bigg{[}\frac{33}{4}r^{4}\dot{\phi}^{4}+\frac{471}{16}r^{2}\dot{r}^{2}\dot{\phi}^{2}-6\dot{r}^{4}+\frac{9}{8}Gmr\dot{\phi}^{2}-\frac{123}{8}\frac{Gm\dot{r}^{2}}{r}+\frac{G^{2}m^{2}}{r^{2}}\Bigg{]}$
$\displaystyle\quad+\frac{32G^{5}m^{6}\nu^{2}\alpha_{\text{EM}}\tilde{\mu}_{+}^{2}}{105c^{11}r^{6}}\Bigg{[}10r^{4}\dot{\phi}^{4}+\frac{83}{4}r^{2}\dot{r}^{2}\dot{\phi}^{2}+4\dot{r}^{4}+Gmr\dot{\phi}^{2}+\frac{4Gm\dot{r}^{2}}{r}+\frac{G^{2}m^{2}}{r^{2}}\Bigg{]}\,.$
(92b)

Note that we have written $\mathcal{F}_{\text{GW,mag}}$ in terms of
$\tilde{\mu}_{1,2}$ instead of $\tilde{\mu}_{\pm}$, for the sake of
compactness. Indeed, it comes

$\tilde{\mu}_{1}\tilde{\mu}_{2}=\nu\left(\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}\right)\,.$
(93)

Similarly to the energy flux, the total angular momentum flux is split
according to

$\mathcal{G}^{i}=\bigg{[}\mathcal{G}_{\text{GW,pp}}+\mathcal{G}_{\text{GW,mag}}+\mathcal{G}_{\text{EM}}\bigg{]}l^{i}\,,$
(94)

where the (instantaneous) point-particle sector, $\mathcal{G}_{\text{GW,pp}}$,
can be found at 3PN _e.g._ in Eq. (3.2) of [105]. At the required order, the
contributions read

$\displaystyle\mathcal{G}_{\text{GW,pp}}$
$\displaystyle=\frac{16G^{2}m^{3}\nu^{2}\dot{\phi}}{5c^{5}r}\Bigg{\\{}r^{2}\dot{\phi}^{2}-\frac{\dot{r}^{2}}{2}+\frac{Gm}{r}$
$\displaystyle\hskip
79.6678pt+\frac{1}{c^{2}}\Bigg{[}\bigg{(}\frac{307}{168}-\frac{137\nu}{42}\bigg{)}r^{4}\dot{\phi}^{4}+\bigg{(}\frac{85}{84}+\frac{283\nu}{84}\bigg{)}r^{2}\dot{r}^{2}\dot{\phi}^{2}+\bigg{(}\frac{37}{42}+\frac{17\nu}{84}\bigg{)}\dot{r}^{4}$
$\displaystyle\hskip
113.81102pt-\bigg{(}\frac{29}{21}+\frac{95\nu}{42}\bigg{)}Gmr\dot{\phi}^{2}+\bigg{(}\frac{64}{21}+\frac{\nu}{12}\bigg{)}\frac{Gm\dot{r}^{2}}{r}-\bigg{(}\frac{745}{84}-\frac{\nu}{42}\bigg{)}\frac{G^{2}m^{2}}{r^{2}}\Bigg{]}\Bigg{\\}}\,,$
(95a) $\displaystyle\mathcal{G}_{\text{GW,mag}}$
$\displaystyle=-\frac{48G^{4}m^{5}\nu\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}\,\dot{\phi}}{5c^{9}r^{3}}\Bigg{\\{}r^{2}\dot{\phi}^{2}-\frac{3\dot{r}^{2}}{2}+\frac{2Gm}{r}$
$\displaystyle\hskip
142.26378pt+\frac{1}{c^{2}}\Bigg{[}\bigg{(}\frac{323}{168}-\frac{22\nu}{3}\bigg{)}r^{4}\dot{\phi}^{4}-\bigg{(}\frac{167}{14}-\frac{1221\nu}{28}\bigg{)}r^{2}\dot{r}^{2}\dot{\phi}^{2}$
$\displaystyle\hskip
170.71652pt+\bigg{(}\frac{27}{7}-\frac{1181\nu}{84}\bigg{)}\dot{r}^{4}+\bigg{(}\frac{125}{84}-\frac{229\nu}{21}\bigg{)}Gmr\dot{\phi}^{2}$
$\displaystyle\hskip
170.71652pt+\bigg{(}\frac{215}{84}+\frac{227\nu}{42}\bigg{)}\frac{Gm\dot{r}^{2}}{r}-\bigg{(}\frac{356}{21}-\frac{116\nu}{21}\bigg{)}\frac{G^{2}m^{2}}{r^{2}}\Bigg{]}\Bigg{\\}}\,,$
(95b) $\displaystyle\mathcal{G}_{\text{EM}}$
$\displaystyle=\frac{4G^{5}m^{6}\nu^{2}\alpha_{\text{EM}}\tilde{\mu}_{-}^{2}\,\dot{\phi}}{15c^{9}r^{4}}\Bigg{\\{}1+\frac{1}{c^{2}}\Bigg{[}\bigg{(}\frac{223}{14}-\frac{282\nu}{7}\bigg{)}r^{2}\dot{\phi}^{2}-\bigg{(}\frac{151}{28}+\frac{207\nu}{14}\bigg{)}\dot{r}^{2}$
$\displaystyle\hskip
170.71652pt-\bigg{(}\frac{159}{28}+\frac{55\nu}{14}\bigg{)}\frac{Gm}{r}\Bigg{]}\Bigg{\\}}$
$\displaystyle\quad+\frac{614G^{5}m^{6}\nu^{2}\delta\alpha_{\text{EM}}\tilde{\mu}_{+}\tilde{\mu}_{-}\,\dot{\phi}}{105c^{11}r^{4}}\Bigg{[}r^{2}\dot{\phi}^{2}-\frac{77\dot{r}^{2}}{307}+\frac{25Gm}{307r}\Bigg{]}$
$\displaystyle\quad+\frac{16G^{5}m^{6}\nu^{2}\alpha_{\text{EM}}\tilde{\mu}_{+}^{2}\,\dot{\phi}}{5c^{11}r^{4}}\Bigg{[}r^{2}\dot{\phi}^{2}-\frac{\dot{r}^{2}}{14}+\frac{Gm}{7r}\Bigg{]}\,.$
(95c)

###  B.2 Expressions of the gravitational modes for generic orbits

This appendix presents the instantaneous gravitational modes computed in Sec
IV.2. For the sake of simplicity, we single out the dependence on the
gravitational phase $\phi$ and rescale the modes $h_{\ell\mathrm{m}}$ (36) as

$h_{\ell\mathrm{m}}=\frac{8Gm\nu}{c^{4}R}\,\sqrt{\frac{\pi}{5}}\,\mathcal{H}_{\ell\mathrm{m}}\,\mathrm{e}^{-\mathrm{i}\mathrm{m}\phi}\,,$
(96)

and further split
$\mathcal{H}_{\ell\mathrm{m}}=\mathcal{H}_{\ell\mathrm{m}}^{\text{pp}}+\mathcal{H}_{\ell\mathrm{m}}^{\text{EM}}$.
Moreover, those coefficients obey
$\mathcal{H}_{\ell,-\mathrm{m}}=(-)^{\ell}\,\overline{\mathcal{H}}_{\ell\mathrm{m}}$,
so we only present the ones with a positive $\mathrm{m}$. The point-particle
sector, $\mathcal{H}_{\ell\mathrm{m}}^{\text{pp}}$ 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

$\displaystyle\mathcal{H}_{22}^{\text{EM}}=-\frac{3}{2}\frac{G^{3}m^{3}\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{c^{4}r^{3}\,\nu}\Bigg{\\{}1+\frac{1}{c^{2}}\Bigg{[}\bigg{(}\frac{1}{21}-\frac{8\nu}{7}\bigg{)}r^{2}\dot{\phi}^{2}-\mathrm{i}\bigg{(}\frac{73}{21}-\frac{73\nu}{7}\bigg{)}r\,\dot{r}\dot{\phi}$
$\displaystyle\hskip
199.16928pt-\bigg{(}\frac{9}{14}-\frac{24\nu}{7}\bigg{)}\dot{r}^{2}-\bigg{(}\frac{187}{42}-\frac{20\nu}{7}\bigg{)}\frac{Gm}{r}\Bigg{]}\Bigg{\\}}\,,$
(97a)
$\displaystyle\mathcal{H}_{21}^{\text{EM}}=-\mathrm{i}\,\frac{G^{3}m^{3}\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}\,\delta}{c^{5}r^{2}\,\nu}\,\dot{\phi}\,,$
(97b)
$\displaystyle\mathcal{H}_{20}^{\text{EM}}=\sqrt{\frac{3}{2}}\frac{G^{3}m^{3}\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{c^{4}r^{3}\,\nu}\Bigg{\\{}1+\frac{1}{c^{2}}\Bigg{[}\bigg{(}\frac{3}{7}-\frac{16\nu}{7}\bigg{)}r^{2}\dot{\phi}^{2}-\bigg{(}\frac{9}{14}-\frac{24\nu}{7}\bigg{)}\dot{r}^{2}-\bigg{(}\frac{187}{42}-\frac{20\nu}{7}\bigg{)}\frac{Gm}{r}\Bigg{]}\Bigg{\\}}\,,$
(97c)
$\displaystyle\mathcal{H}_{33}^{\text{EM}}=\frac{\mathrm{i}}{4}\sqrt{\frac{105}{2}}\frac{G^{3}m^{3}\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}\,\delta}{c^{5}r^{3}\,\nu}\bigg{(}r\dot{\phi}+\frac{2}{7}\,\mathrm{i}\,\dot{r}\bigg{)}\,,$
(97d)
$\displaystyle\mathcal{H}_{32}^{\text{EM}}=-\sqrt{\frac{5}{7}}\frac{G^{3}m^{3}\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{c^{6}r^{2}\,\nu}\big{(}1-3\nu\big{)}\,\bigg{(}r\dot{\phi}-\frac{\mathrm{i}\,\dot{r}}{4}\bigg{)}\,\dot{\phi}\,,$
(97e)
$\displaystyle\mathcal{H}_{31}^{\text{EM}}=-\frac{\mathrm{i}}{4}\sqrt{\frac{7}{2}}\frac{G^{3}m^{3}\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}\,\delta}{c^{5}r^{3}\,\nu}\bigg{(}r\dot{\phi}+\frac{6}{7}\,\mathrm{i}\,\dot{r}\bigg{)}\,,$
(97f)
$\displaystyle\mathcal{H}_{30}^{\text{EM}}=-\frac{\mathrm{i}}{2}\sqrt{\frac{3}{14}}\frac{G^{3}m^{3}\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{c^{6}r^{2}\,\nu}\big{(}1-3\nu\big{)}\,\dot{r}\,\dot{\phi}\,,$
(97g)
$\displaystyle\mathcal{H}_{44}^{\text{EM}}=\frac{53}{24}\sqrt{\frac{5}{7}}\frac{G^{3}m^{3}\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{c^{6}r^{3}\,\nu}\big{(}1-3\nu\big{)}\,\bigg{[}r^{2}\dot{\phi}^{2}+\frac{26\,\mathrm{i}}{53}r\,\dot{r}\dot{\phi}-\frac{8}{53}\dot{r}^{2}+\frac{12}{53}\frac{Gm}{r}\bigg{]}\,,$
(97h)
$\displaystyle\mathcal{H}_{42}^{\text{EM}}=-\frac{5\sqrt{5}}{84}\frac{G^{3}m^{3}\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{c^{6}r^{3}\,\nu}\big{(}1-3\nu\big{)}\,\bigg{[}r^{2}\dot{\phi}^{2}+\frac{13\,\mathrm{i}}{5}r\,\dot{r}\dot{\phi}-\frac{8}{5}\dot{r}^{2}+\frac{12}{5}\frac{Gm}{r}\bigg{]}\,,$
(97i)
$\displaystyle\mathcal{H}_{40}^{\text{EM}}=-\frac{11}{28\sqrt{2}}\frac{G^{3}m^{3}\alpha_{\text{EM}}\tilde{\mu}_{1}\tilde{\mu}_{2}}{c^{6}r^{3}\,\nu}\big{(}1-3\nu\big{)}\,\bigg{[}r^{2}\dot{\phi}^{2}+\frac{8}{11}\dot{r}^{2}-\frac{12}{11}\frac{Gm}{r}\bigg{]}\,.$
(97j)

###  B.3 Quasi-Keplerian parameters

The parameters entering the quasi-Keplerian representation of the motion (63)
are given herebelow in terms of the two dimensionless invariants linked to the
(unreduced) energy $E$ and norm of the angular momentum $J=|J^{i}|$

$\varepsilon\equiv\frac{-2E}{m\,\nu\,c^{2}}\qquad\text{and}\qquad
j\equiv\frac{-2E\,J^{2}}{G^{2}m^{5}\,\nu^{3}}\,,$ (98)

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_{\star}=\sqrt{1-j}+\mathcal{O}(c^{-2})$. It relates to the usual definition
$h=J/(Gm^{2}\nu)$ by $j=-2h^{2}E/(m\nu)$. At the required order, we obtain

$\displaystyle
a_{r}=\frac{Gm}{c^{2}\,\varepsilon}\Bigg{\\{}1-\frac{7-\nu}{4}\varepsilon+\frac{\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{j\,\nu}\,\varepsilon^{2}\,\bigg{[}1+\frac{12+\nu}{j}\,\varepsilon\bigg{]}\Bigg{\\}}\,,$
(99a) $\displaystyle
n=\frac{c^{3}}{Gm}\,\varepsilon^{3/2}\,\Bigg{\\{}1-\frac{15-\nu}{8}\varepsilon+\big{(}7-2\nu\big{)}\frac{\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{j^{3/2}\nu}\,\varepsilon^{3}\Bigg{\\}}\,,$
(99b) $\displaystyle
K=1+\frac{3}{j}\,\varepsilon-\frac{3\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{j^{2}\,\nu}\,\varepsilon^{2}\bigg{[}1+\bigg{(}\frac{45}{2j}-\frac{7}{2}+\nu\bigg{)}\,\varepsilon\bigg{]}\,,$
(99c) $\displaystyle
e_{r}^{2}=1-j+\bigg{(}6-\nu-\frac{15-5\nu}{4}\,j\bigg{)}\varepsilon$
$\displaystyle\hskip
56.9055pt-\frac{2\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{j\,\nu}\,\varepsilon^{2}\,\bigg{[}2-j+\bigg{(}\frac{24+2\nu}{j}-\frac{11+7\nu}{2}-\frac{11-3\nu}{2}\,j\bigg{)}\varepsilon\bigg{]}\,,$
(99d) $\displaystyle
e_{t}^{2}=1-j-\bigg{(}2-2\nu-\frac{17-7\nu}{4}\,j\bigg{)}\varepsilon$
$\displaystyle\hskip
56.9055pt-\frac{2\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{j\,\nu}\,\varepsilon^{2}\,\bigg{[}1+\bigg{(}\frac{12+\nu}{j}-\frac{7-2\nu}{j^{1/2}}-\frac{37-7\nu}{4}+(7-2\nu)j^{1/2}\bigg{)}\varepsilon\bigg{]}\,,$
(99e) $\displaystyle
e_{\phi}^{2}=1-j+\bigg{(}6-\frac{15-\nu}{4}\,j\bigg{)}\varepsilon$
$\displaystyle\hskip
56.9055pt-\frac{2\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{j\,\nu}\,\varepsilon^{2}\,\bigg{[}3-2j+\bigg{(}\frac{336-13\nu}{8\,j}-\frac{63-16\nu}{4}-\frac{88+9\nu}{8}\,j\bigg{)}\varepsilon\bigg{]}\,,$
(99f) $\displaystyle
f_{v-u}=-\big{(}7-2\nu)\frac{\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{j^{3/2}\,\nu}\,\varepsilon^{3}\,,$
(99g) $\displaystyle
f_{v}=\frac{3\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{2}\,\frac{\sqrt{1-j}}{j^{3/2}}\,\varepsilon^{3}\,,$
(99h) $\displaystyle
g_{2v}=-\frac{3-4\nu}{4}\frac{\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{\nu}\frac{1-j}{j^{3}}\,\varepsilon^{3}\,,$
(99i) $\displaystyle
g_{3v}=\frac{\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{8}\frac{(1-j)^{3/2}}{j^{3}}\,\varepsilon^{3}\,.$
(99j)

The 1PN point-particle sector of those expressions naturally agrees with [92].
For the sake of completeness, we also display the expression of $x$ (64) at
the sought order

$x=\varepsilon\,\Bigg{\\{}1+\bigg{[}\frac{2}{j}-\frac{15-\nu}{12}\bigg{]}\,\varepsilon-\frac{2\alpha_{\text{EM}}\,\tilde{\mu}_{1}\tilde{\mu}_{2}}{j^{2}\,\nu}\,\varepsilon^{2}\bigg{[}1+\bigg{(}\frac{43}{2j}-\frac{57-13\nu}{12}-\frac{7-2\nu}{3}j^{1/2}\bigg{)}\varepsilon\bigg{]}\Bigg{\\}}\,.$
(100)

###  B.4 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

$\displaystyle\left\langle\mathcal{F}\right\rangle=\frac{32\,\nu^{2}c^{5}}{5G}x^{5}\Bigg{[}\mathcal{F}_{\text{N}}+x\,\mathcal{F}_{1\text{PN}}+\alpha_{\text{EM}}\,x^{2}\bigg{(}\mathcal{F}_{\text{mag,LO}}+x\,\mathcal{F}_{\text{mag,NLO}}\bigg{)}\Bigg{]}\,,$
(101a)
$\displaystyle\left\langle\mathcal{G}\right\rangle=\frac{4}{5}\,m\nu^{2}c^{2}\,x^{7/2}\Bigg{[}\mathcal{G}_{\text{N}}+x\,\mathcal{G}_{1\text{PN}}+\alpha_{\text{EM}}\,x^{2}\bigg{(}\mathcal{G}_{\text{mag,LO}}+x\,\mathcal{G}_{\text{mag,NLO}}\bigg{)}\Bigg{]}\,,\,.$
(101b)

with (let us emphasize that the following expressions are exact, in the sense
that they have not been truncated in powers of $e_{t}$)

$\displaystyle\mathcal{F}_{\text{N}}$
$\displaystyle=\frac{1}{(1-e_{t}^{2})^{7/2}}\Bigg{[}1+\frac{73\,e_{t}^{2}}{24}+\frac{37\,e_{t}^{4}}{96}\Bigg{]}\,,$
(102a) $\displaystyle\mathcal{F}_{1\text{PN}}$
$\displaystyle=\frac{1}{(1-e_{t}^{2})^{9/2}}\Bigg{[}-\frac{1247}{336}-\frac{35\nu}{12}+\bigg{(}\frac{10475}{672}-\frac{1081\nu}{36}\bigg{)}e_{t}^{2}$
$\displaystyle\hskip
85.35826pt+\bigg{(}\frac{10043}{384}-\frac{311\nu}{12}\bigg{)}e_{t}^{4}+\bigg{(}\frac{2179}{1792}-\frac{851\nu}{576}\bigg{)}e_{t}^{6}\Bigg{]}\,,$
(102b) $\displaystyle\mathcal{F}_{\text{mag,LO}}$
$\displaystyle=\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{11/2}}\Bigg{[}-4-\frac{101\,e_{t}^{2}}{2}-\frac{1999\,e_{t}^{4}}{48}-\frac{59\,e_{t}^{6}}{32}\Bigg{]}$
$\displaystyle\quad+\frac{\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{11/2}}\bigg{[}\frac{1}{24}+\frac{19\,e_{t}^{2}}{48}+\frac{23\,e_{t}^{4}}{64}+\frac{3\,e_{t}^{6}}{128}\Bigg{]}\,,$
(102c) $\displaystyle\mathcal{F}_{\text{mag,NLO}}$
$\displaystyle=\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{13/2}}\Bigg{[}\frac{663}{56}+\frac{47\nu}{2}-\bigg{(}\frac{28579}{84}-\frac{4903\nu}{9}\bigg{)}e_{t}^{2}-\bigg{(}\frac{138987}{112}-\frac{50419\nu}{36}\bigg{)}e_{t}^{4}$
$\displaystyle\hskip
142.26378pt-\bigg{(}\frac{318767}{672}-\frac{157811\nu}{288}\bigg{)}e_{t}^{6}-\bigg{(}\frac{69525}{7168}-\frac{1537\nu}{96}\bigg{)}e_{t}^{8}\Bigg{]}$
$\displaystyle\quad+\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{6}}\Bigg{[}-\frac{49}{3}+\frac{14\nu}{3}-\bigg{(}\frac{8995}{72}-\frac{1285\nu}{36}\bigg{)}e_{t}^{2}$
$\displaystyle\hskip
142.26378pt-\bigg{(}\frac{7091}{96}-\frac{1013\nu}{48}\bigg{)}e_{t}^{4}-\bigg{(}\frac{259}{144}-\frac{37\nu}{72}\bigg{)}e_{t}^{6}\Bigg{]}$
$\displaystyle\quad+\frac{\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}}{(1-e_{t}^{2})^{13/2}}\Bigg{[}\frac{5}{12}+\frac{259\,e_{t}^{2}}{48}+\frac{1963\,e_{t}^{4}}{192}+\frac{1441\,e_{t}^{6}}{384}+\frac{141\,e_{t}^{8}}{1024}\Bigg{]}$
$\displaystyle\quad+\frac{\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{13/2}}\Bigg{[}\frac{19}{32}-\frac{191\nu}{144}+\bigg{(}\frac{93}{8}-\frac{2627\nu}{144}\bigg{)}e_{t}^{2}+\bigg{(}\frac{10915}{384}-\frac{10333\nu}{288}\bigg{)}e_{t}^{4}$
$\displaystyle\hskip
113.81102pt+\bigg{(}\frac{1079}{96}-\frac{5099\nu}{384}\bigg{)}e_{t}^{6}+\bigg{(}\frac{1533}{4096}-\frac{983\nu}{2048}\bigg{)}e_{t}^{8}\Bigg{]}\,,$
(102d)

and

$\displaystyle\mathcal{G}_{\text{N}}$
$\displaystyle=\frac{8+7e_{t}^{2}}{(1-e_{t}^{2})^{2}}\,,$ (103a)
$\displaystyle\mathcal{G}_{1\text{PN}}$
$\displaystyle=\frac{1}{(1-e_{t}^{2})^{3}}\Bigg{[}-\frac{1247}{42}-\frac{70\nu}{3}+\bigg{(}\frac{3019}{42}-\frac{335\nu}{3}\bigg{)}e_{t}^{2}+\bigg{(}\frac{8399}{336}-\frac{275\nu}{12}\bigg{)}e_{t}^{4}\Bigg{]}\,,$
(103b) $\displaystyle\mathcal{G}_{\text{mag,LO}}$
$\displaystyle=\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{4}}\Bigg{[}-32-174\,e_{t}^{2}-\frac{63\,e_{t}^{4}}{2}\Bigg{]}+\frac{\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{4}}\bigg{[}\frac{1}{3}+e_{t}^{2}+\frac{e_{t}^{4}}{8}\Bigg{]}\,,$
(103c) $\displaystyle\mathcal{G}_{\text{mag,NLO}}$
$\displaystyle=\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{5}}\Bigg{[}\frac{271}{7}+204\nu-\bigg{(}\frac{63443}{42}-\frac{7192\nu}{3}\bigg{)}e_{t}^{2}$
$\displaystyle\hskip
142.26378pt-\bigg{(}\frac{363301}{168}-\frac{15529\nu}{6}\bigg{)}e_{t}^{4}-\bigg{(}\frac{17091}{112}-206\nu\bigg{)}e_{t}^{6}\Bigg{]}$
$\displaystyle\quad-\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{11/2}}\,\frac{2(7-2\nu)(16+41\,e_{t}^{2}-60\,e_{t}^{4}-7\,e_{t}^{6})}{3}$
$\displaystyle\quad+\frac{\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}}{(1-e_{t}^{2})^{5}}\Bigg{[}\frac{10}{3}+\frac{145\,e_{t}^{2}}{6}+\frac{35\,e_{t}^{4}}{2}+\frac{15\,e_{t}^{6}}{16}\Bigg{]}$
$\displaystyle\quad+\frac{\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{5}}\Bigg{[}\frac{19}{4}-\frac{191\nu}{18}+\bigg{(}\frac{385}{8}-\frac{2867\nu}{36}\bigg{)}e_{t}^{2}$
$\displaystyle\hskip
85.35826pt+\bigg{(}\frac{1263}{32}-\frac{2845\nu}{48}\bigg{)}e_{t}^{4}+\bigg{(}\frac{131}{64}-\frac{311\nu}{96}\bigg{)}e_{t}^{6}\Bigg{]}\,.$
(103d)

Naturally, the point-particle sectors at 1PN agree with Eqs (8.9) of [98] and
(4.11) of [105].

###  B.5 Time evolution of the orbital parameters

Using the chain rule (67) together with Eq. (99) and the averaged fluxes
(101), we find the evolution of the time-eccentricity as

$\left\langle\frac{\mathrm{d}e_{t}}{\mathrm{d}t}\right\rangle=-\frac{c^{3}\,e_{t}\,\nu}{Gm}\,x^{4}\,\Bigg{[}\mathcal{E}_{\text{N}}+x\,\mathcal{E}_{1\text{PN}}+\alpha_{\text{EM}}\,x^{2}\bigg{(}\mathcal{E}_{\text{mag,LO}}+x\,\mathcal{E}_{\text{mag,NLO}}\bigg{)}\Bigg{]}\,,$
(104)

where

$\displaystyle\mathcal{E}_{\text{N}}$
$\displaystyle=\frac{1}{(1-e_{t}^{2})^{5/2}}\Bigg{[}\frac{304}{15}+\frac{121\,e_{t}^{2}}{15}\Bigg{]}\,,$
(105a) $\displaystyle\mathcal{E}_{1\text{PN}}$
$\displaystyle=\frac{1}{(1-e_{t}^{2})^{7/2}}\Bigg{[}-\frac{939}{35}-\frac{4084\nu}{45}+\bigg{(}\frac{29917}{105}-\frac{7753\nu}{30}\bigg{)}e_{t}^{2}+\bigg{(}\frac{13929}{280}-\frac{1664\nu}{45}\bigg{)}e_{t}^{4}\Bigg{]}\,,$
(105b) $\displaystyle\mathcal{E}_{\text{mag,LO}}$
$\displaystyle=\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{9/2}}\Bigg{[}-\frac{3752}{15}-\frac{1190\,e_{t}^{2}}{3}-37\,e_{t}^{4}\Bigg{]}+\frac{\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{9/2}}\bigg{[}2+3\,e_{t}^{2}+\frac{e_{t}^{4}}{4}\Bigg{]}\,,$
(105c) $\displaystyle\mathcal{E}_{\text{mag,NLO}}$
$\displaystyle=\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{11/2}}\Bigg{[}-\frac{16220}{21}+\frac{72916\nu}{45}-\bigg{(}\frac{1778741}{210}-\frac{27724\nu}{3}\bigg{)}e_{t}^{2}$
$\displaystyle\hskip
142.26378pt-\bigg{(}\frac{267145}{42}-\frac{543343\nu}{90}\bigg{)}e_{t}^{4}-\bigg{(}\frac{321277}{1120}-\frac{3099\nu}{10}\bigg{)}e_{t}^{6}\Bigg{]}$
$\displaystyle\quad-\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{5}}\,\frac{(7-2\nu)(3248+6891\,e_{t}^{2}+61\,e_{t}^{4})}{45}$
$\displaystyle\quad+\frac{\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}}{(1-e_{t}^{2})^{11/2}}\Bigg{[}\frac{268}{15}+\frac{2123\,e_{t}^{2}}{30}+\frac{559\,e_{t}^{4}}{15}+\frac{261\,e_{t}^{6}}{160}\Bigg{]}$
$\displaystyle\quad+\frac{\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{11/2}}\Bigg{[}\frac{1013}{30}-\frac{872\nu}{15}+\bigg{(}\frac{7659}{40}-\frac{14677\nu}{60}\bigg{)}e_{t}^{2}$
$\displaystyle\hskip
113.81102pt+\bigg{(}\frac{27527}{240}-\frac{2683\nu}{20}\bigg{)}e_{t}^{4}+\bigg{(}\frac{641}{128}-\frac{1927\nu}{320}\bigg{)}e_{t}^{6}\Bigg{]}\,.$
(105d)

Let us 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 as

$\left\langle\frac{\mathrm{d}a_{r}}{\mathrm{d}t}\right\rangle=\nu\,c\,x^{3}\,\Bigg{[}\mathcal{A}_{\text{N}}+x\,\mathcal{A}_{1\text{PN}}+\alpha_{\text{EM}}\,x^{2}\bigg{(}\mathcal{A}_{\text{mag,LO}}+x\,\mathcal{A}_{\text{mag,NLO}}\bigg{)}\Bigg{]}\,,$
(106)

where

$\displaystyle\mathcal{A}_{\text{N}}$
$\displaystyle=\frac{1}{(1-e_{t}^{2})^{7/2}}\Bigg{[}-\frac{64}{5}-\frac{584\,e_{t}^{2}}{15}-\frac{74\,e_{t}^{4}}{15}\Bigg{]}\,,$
(107a) $\displaystyle\mathcal{A}_{1\text{PN}}$
$\displaystyle=\frac{1}{(1-e_{t}^{2})^{9/2}}\Bigg{[}\frac{2972}{105}+\frac{176\nu}{5}-\bigg{(}\frac{30442}{105}-380\nu\bigg{)}e_{t}^{2}$
$\displaystyle\hskip
85.35826pt-\bigg{(}\frac{879}{2}-\frac{1687\nu}{5}\bigg{)}e_{t}^{4}-\bigg{(}\frac{11717}{420}-\frac{296\nu}{15}\bigg{)}e_{t}^{6}\Bigg{]}\,,$
(107b) $\displaystyle\mathcal{A}_{\text{mag,LO}}$
$\displaystyle=\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{11/2}}\Bigg{[}\frac{576}{5}+\frac{11816\,e_{t}^{2}}{15}+\frac{2716\,e_{t}^{4}}{5}+\frac{174\,e_{t}^{6}}{5}\Bigg{]}$
$\displaystyle\quad-\frac{\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{11/2}}\bigg{[}\frac{8}{15}+\frac{76\,e_{t}^{2}}{15}+\frac{23\,e_{t}^{4}}{5}+\frac{3\,e_{t}^{6}}{10}\Bigg{]}\,,$
(107c) $\displaystyle\mathcal{A}_{\text{mag,NLO}}$
$\displaystyle=\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{13/2}}\Bigg{[}\frac{1156}{3}-\frac{2192\nu}{5}+\bigg{(}\frac{171160}{21}-\frac{78572\nu}{9}\bigg{)}e_{t}^{2}+\bigg{(}\frac{125977}{6}-\frac{884126\nu}{45}\bigg{)}e_{t}^{4}$
$\displaystyle\hskip
142.26378pt+\bigg{(}\frac{55471}{7}-\frac{109891\nu}{15}\bigg{)}e_{t}^{6}+\bigg{(}\frac{459559}{1680}-\frac{4103\nu}{15}\bigg{)}e_{t}^{8}\Bigg{]}$
$\displaystyle\quad+\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{6}}\,\frac{2(7-2\nu)(96+1452\,e_{t}^{2}+1353\,e_{t}^{4}+74\,e_{t}^{6})}{15}$
$\displaystyle\quad-\frac{\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}}{(1-e_{t}^{2})^{13/2}}\Bigg{[}\frac{16}{3}+\frac{1036\,e_{t}^{2}}{15}+\frac{1963\,e_{t}^{4}}{15}+\frac{1441\,e_{t}^{6}}{30}+\frac{141\,e_{t}^{8}}{80}\Bigg{]}$
$\displaystyle\quad-\frac{\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{13/2}}\Bigg{[}\frac{42}{5}-\frac{152\nu}{9}+\bigg{(}\frac{2366}{15}-\frac{10474\nu}{45}\bigg{)}e_{t}^{2}+\bigg{(}\frac{1917}{5}-\frac{41339\nu}{90}\bigg{)}e_{t}^{4}$
$\displaystyle\hskip
113.81102pt+\bigg{(}\frac{9349}{60}-\frac{10241\nu}{60}\bigg{)}e_{t}^{6}+\bigg{(}\frac{1773}{320}-\frac{991\nu}{160}\bigg{)}e_{t}^{8}\Bigg{]}\,.$
(107d)

Quite naturally, the Newtonian factors of Eqs. (104) and (106) are nothing but
the well-known 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

$\left\langle\frac{\mathrm{d}x}{\mathrm{d}t}\right\rangle=\frac{c^{3}\,\nu}{Gm}\,x^{5}\,\Bigg{[}\mathcal{X}_{\text{N}}+x\,\mathcal{X}_{1\text{PN}}+\alpha_{\text{EM}}\,x^{2}\bigg{(}\mathcal{X}_{\text{mag,LO}}+x\,\mathcal{X}_{\text{mag,NLO}}\bigg{)}\Bigg{]}\,,$
(108)

with

$\displaystyle\mathcal{X}_{\text{N}}$
$\displaystyle=\frac{1}{(1-e_{t}^{2})^{7/2}}\Bigg{[}\frac{64}{5}+\frac{584\,e_{t}^{2}}{15}+\frac{74\,e_{t}^{4}}{15}\Bigg{]}\,,$
(109a) $\displaystyle\mathcal{X}_{1\text{PN}}$
$\displaystyle=\frac{1}{(1-e_{t}^{2})^{9/2}}\Bigg{[}-\frac{2972}{105}-\frac{176\nu}{5}+\bigg{(}\frac{1462}{7}-380\nu\bigg{)}e_{t}^{2}$
$\displaystyle\hskip
85.35826pt+\bigg{(}\frac{12217}{30}-\frac{1687\nu}{5}\bigg{)}e_{t}^{4}+\bigg{(}\frac{11717}{420}-\frac{296\nu}{15}\bigg{)}e_{t}^{6}\Bigg{]}\,,$
(109b) $\displaystyle\mathcal{X}_{\text{mag,LO}}$
$\displaystyle=\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{11/2}}\Bigg{[}-128-\frac{10768\,e_{t}^{2}}{15}-\frac{7472\,e_{t}^{4}}{15}-\frac{118\,e_{t}^{6}}{5}\Bigg{]}$
$\displaystyle\quad+\frac{\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{11/2}}\bigg{[}\frac{8}{15}+\frac{76\,e_{t}^{2}}{15}+\frac{23\,e_{t}^{4}}{5}+\frac{3\,e_{t}^{6}}{10}\Bigg{]}\,,$
(109c) $\displaystyle\mathcal{X}_{\text{mag,NLO}}$
$\displaystyle=\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{13/2}}\Bigg{[}-\frac{18768}{35}+\frac{2368\nu}{5}-\bigg{(}6812-\frac{382784\nu}{45}\bigg{)}e_{t}^{2}-\bigg{(}\frac{366455}{21}-\frac{828794\nu}{45}\bigg{)}e_{t}^{4}$
$\displaystyle\hskip
142.26378pt-\bigg{(}\frac{203953}{30}-\frac{306497\nu}{45}\bigg{)}e_{t}^{6}-\bigg{(}\frac{20513}{112}-\frac{3133\nu}{15}\bigg{)}e_{t}^{8}\Bigg{]}$
$\displaystyle\quad+\frac{\tilde{\mu}_{+}^{2}+\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}-\nu\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{6}}\,\frac{2(7-2\nu)(96-4484\,e_{t}^{2}-4530\,e_{t}^{4}-7\,e_{t}^{6})}{45}$
$\displaystyle\quad+\frac{\delta\,\tilde{\mu}_{+}\tilde{\mu}_{-}}{(1-e_{t}^{2})^{13/2}}\Bigg{[}\frac{16}{3}+\frac{1036\,e_{t}^{2}}{15}+\frac{1963\,e_{t}^{4}}{15}+\frac{1441\,e_{t}^{6}}{30}+\frac{141\,e_{t}^{8}}{80}\Bigg{]}$
$\displaystyle\quad+\frac{\tilde{\mu}_{-}^{2}}{(1-e_{t}^{2})^{13/2}}\Bigg{[}\frac{42}{5}-\frac{152\nu}{9}+\bigg{(}\frac{2246}{15}-\frac{10474\nu}{45}\bigg{)}e_{t}^{2}+\bigg{(}\frac{1857}{5}-\frac{41339\nu}{90}\bigg{)}e_{t}^{4}$
$\displaystyle\hskip
113.81102pt+\bigg{(}\frac{9289}{60}-\frac{10241\nu}{60}\bigg{)}e_{t}^{6}+\bigg{(}\frac{1773}{320}-\frac{991\nu}{160}\bigg{)}e_{t}^{8}\Bigg{]}\,,$
(109d)

and we recognize the Newtonian sector as the “enhancement” factor due to
Peters and Mathews [99].

## References

  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run (2021), arXiv:2111.03606 [gr-qc] .