# Final velocity and radiated energy in numerical simulations of binary black
holes

Emmanuel A. Tassone<EMAIL_ADDRESS>Carlos N. Kozameh
<EMAIL_ADDRESS>FaMAF, Universidad Nacional de Córdoba
5000, Córdoba, Argentina

###### Abstract

The evolution of global binary black holes variables such as energy or linear
momentum are mainly obtained by applying numerical methods near
coalescence,post-Newtonian expansions or a combination of both. In this paper,
we use a fully relativistic formalism presented several years ago that only
uses global variables defined at null infinity together with the gravitational
radiation emitted by the source to obtain the time evolution of such variables
for binary black holes (BBH) systems. For that, we use the Rochester catalog
composed of 776 BBHs simulations. We compute the final velocity, radiated
energy, and intrinsic angular momentum predicted by the dynamical equations in
this formalism for non-spinning, aligned/anti-aligned spins, and several
different precessing configurations. We compare obtained values with reported
values in numerical simulations. As BBHs parameter space is still not
completely covered by numerical simulations, we fit phenomenological formulas
for practical applications to the radiated energy and final velocities
obtained. Also, we compare the fits with reported values. In conclusion, we
see that our formulae and correlations for the variables described in this
work are consistent with those found in the general literature.

††preprint: APS/123-QED

##  I Introduction

The RIT Catalog of Numerical Simulations [1] is a catalog containing the
gravitational wave data for 776 numerical relativity simulations. The catalog
spans a large variety of initial conditions and provides the gravitational
strain for the evolution of the mergers during all the phases.

It is believed that spin orientation on the initial state of the binary black
hole (BBH) can constrain the evolutionary processes that lead the binary to
the merger [2]. Spin orientations can influence the kick velocities [3],
making them high enough for a merged binary to be ejected from the nucleus of
a massive host galaxy [4, 5]. The probability of these large recoils depends
on the distribution of mass ratios and spins of the progenitor binaries [6].
Spin orientations also influence the energy radiated, and particularly, the
maximum luminosity of the merger[2]. Additionally, there is also research made
on the hangup effect, which is a change in the merger duration compared to
non-spinning binaries due to the presence of spins, resulting in faster or
slower fusion of the binaries [7, 8].

Other important parameters when setting up a black hole simulation are the
masses of the black holes which are represented typically by the quotient
$q=m_{1}/m_{2}$, with $m_{1}$ being the smaller black hole progenitor. The
mass quotient, along with the spin orientation, may determine the remnant
mass, spin and recoil velocity [9].

Although the spectrum of parameters along all the 776 simulations in Rochester
catalog is wide, it is still computationally prohibitive to cover the whole
BBH parameter space exhaustively, making fitting formulas useful for practical
purposes. Many fitting formulas have been deduced that relate the final spin
of the remnant black hole with the initial spins and mass parameters [10, 11,
12, 13]. However, most of the relations established from the data analysis are
constructed based on models from the parametrized post-Newtonian formalism
(PPN). The PPN formalism allows finding approximate solutions to the Einstein
field equations for the metric tensor, and therefore, obtaining physical
quantities such as energy, velocity and angular momentum when the binaries are
far away from each other[14].

In this paper, we perform the numerical evolution of four global physical
variables defined at null infinity, namely, the Bondi energy and linear
momentum, the center of mass velocity and intrinsic angular momentum. The last
two variables are defined using a formalism developed by Kozameh and Quiroga
(KQ)[15].

The evolution of the before-mentioned variables is made using the Rochester
gravitational wave strain data catalog. We analyze the dynamics of the final
state of the resulting black hole and relate it to the most significant
variables before the coalescence. Hence, the initial spin orientations,
orbital angular momentum, initial masses, maximum energy lost, final total
velocity, final total angular momentum and its tilting angle are going to be
studied in this work. Many interesting correlations between these variables
are shown in the work.

This paper is organized as follows. In Sec. II, we describe the equations of
motion used for this work. In Sec. III, we mention some aspect of the code
employed to compute the evolution and the subsequent data analysis. In Sec IV,
we analyze the results obtained for the final states of the evolution. Then we
make a deeper analysis, classifying the results into three categories
according to the kinematics of the BBH spins: non spinning (NS), aligned (A)
and precessing (P). There are 29 simulations in first category (NS), 447 with
aligned spins (A) and 300 simulations are precessing spins BBH (P). Finally,
Sec. V closes the work with a discussion and conclusion about the most
relevant aspects found and also a discussion over the formalism employed and
potential future work.

##  II Background

The KQ formalism [15] provides a fully relativistic treatment that defines the
center of mass frame and intrinsic angular momentum using the asymptotic
symmetries represented by the Winicour linkages [16]. It is worth mentioning
that a different formulation based on null congruences with twist by Adamo,
Newman and Kozameh (ANK) [17] yields similar evolution equations at a
quadrupole level. In this sense, the KG formulation has been shown to be
consistent with PN results up to octupole terms in the gravitational radiation
[18].

The physical quantities for asymptotically flat spacetimes are given by the
Weyl scalars $\psi_{0}$, $\psi_{1}$, $\psi_{2}$, $\psi_{3}$, and $\psi_{4}$ at
null infinity. Those scalars are defined on special reference frames called
Newman-Unti coordinates [19]. Each Newman-Unti foliation is then associated to
an arbitrary timelike wordline on a fiducial flat space, called the
Observation space, where the equations of motion are given and numerically
solved. These special frames allow us to define physical variables on ”non-
inertial” frames. Using the Winicour linkages together with a gauge fixing
procedure, the KG formalism defines the center of mass and intrinsic angular
momentum. Using the Bianchi identities at null infinity, one obtains the
evolution equations for those variables on the Observation space. The reader
may refer to [15, 18, 20] for a thorough construction of the physical
quantities using the KG formalism. The relevant evolution equations for the
purpose of this work read

$\displaystyle\dot{M}$
$\displaystyle=-\frac{c}{10\sqrt{2}G}(\dot{\sigma}_{R}^{ij}\dot{\sigma}_{R}^{ij}+\dot{\sigma}_{I}^{ij}\dot{\sigma}_{I}^{ij}),$
(1) $\displaystyle\dot{P}^{i}$
$\displaystyle=\frac{2c^{2}}{15\sqrt{2}G}\dot{\sigma}_{R}^{jl}\dot{\sigma}_{I}^{kl}\epsilon^{ijk},$
(2) $\displaystyle\dot{D}^{i}$
$\displaystyle=\sqrt{2}P^{i}+\frac{3c^{2}}{7G}[\dot{\sigma}_{R}^{ijk}\sigma_{R}^{jk}+\dot{\sigma}_{I}^{ijk}\sigma_{I}^{jk}-\sigma_{R}^{ijk}\dot{\sigma}_{R}^{jk}-\sigma_{I}^{ijk}\dot{\sigma}_{I}^{jk}],$
(3) $\displaystyle\dot{J}^{i}$
$\displaystyle=-\frac{c^{3}}{5G}(\sigma_{R}^{kl}\dot{\sigma}_{R}^{jl}+\sigma_{I}^{kl}\dot{\sigma}_{I}^{jl})\epsilon^{ijk},$
(4)

where $M^{i}$ and $P^{i}$ are the total Bondi mass and linear momentum
respectively, $D^{i}$ is the dipole mass momentum, and $J^{i}$ the total
angular momentum. The first two equations yield the evolution of the Bondi
four-momentum whereas the last two yield the corresponding evolution of the
dipole mass momentum and total angular momentum. These last two variables are
the components of the so-called relativistic angular momentum 2-form and
behave like the electric and magnetic fields under a Lorentz transformation.
The symbols $\sigma_{R}$ and $\sigma_{I}$ denote the real and imaginary part
of the Bondi shear at null infinity. Commonly, the Bondi shear is denoted as
$\sigma_{0}$ and accounts for the deformation of the spacetime due to the
gravitational radiation that arrives at null infinity. Usually $\sigma^{0}$ is
expanded in terms of spin-weighted spherical harmonics $Y^{s}_{lm}$, very
useful for analyzing gravitational waves [21]. Yet, the KQ formalism expands
the asymptotic shear $\sigma^{0}$ in terms of tensor spherical harmonics
$Y^{s}_{i_{1}i_{2}...}$ [22]. This base is fundamental to define covariant
quantities as in eqs. 1, 2, 3 and 4. The dot on the $\sigma^{ij}$ stands for
the derivative with respect to the Bondi time coordinate at future null
infinity. It is worth mentioning that in the above equations we have omitted
quadratic octupole terms and all remaining coefficients for $l\geq 4$ since
they are negligible in BBH coalescence.

Finally we write down the equations for the center of mass $R^{i}$ and
intrinsic angular momentum $S^{i}$ as given in the KG formalism. They are
related to the previous variables by

$\displaystyle D^{i}$
$\displaystyle=MR^{i}+c^{-2}\epsilon^{ijk}\frac{P^{j}}{M}S^{k}-\frac{8}{5\sqrt{2}c}P^{j}\Delta\sigma_{R}^{ij}$
(5) $\displaystyle J^{i}$
$\displaystyle=S^{i}+\epsilon^{ijk}R^{j}P^{k}-\frac{137c^{2}}{168\sqrt{2}G}(\sigma_{R}^{ijk}\sigma_{I}^{jk}-\sigma_{I}^{ijk}\sigma_{R}^{jk})$
(6)

The center of mass $R^{i}$, is the spatial component of a covariant vector of
the formalism and represents the physical analogue of the newtonian center of
mass but also accounting for the radiated energy. One can take the derivative
with respect to the Bondi time $u$ and define the center of mass velocity
$V^{i}$ as the change of the center of mass with respect to the Bondi time.
One of the purposes of this work will be to show that the center-of-mass
velocities given by the KQ formalism agree very well with the recoil
velocities found in the RIT catalog.

The equation for $J^{i}$ resembles its newtonian counterpart except for the
extra radiative terms. In this work the last term is negligible for the BBH
catalog but it could be important in other situations where the octupole term
is not negligible. A fully relativistic version of these equations can be
found in reference [20]. Note that the orbital angular momentum can be
important when the center of mass is displaced from the origin and achieves
high velocities due to the gravitational kick.

##  III  Numerical framework

The numerical code employed to process data and run the evolution equations is
described in this section. The code can be accessed at [23].

The Rochester catalog has been used before to obtain the dynamical evolution
of several physical observables using the KG formalism [24]. However, the
latter work has a numerical error associated to the strain in the last steps
of the evolution. This error is due to a well-known effect of performing a raw
numerical integration [25]. Although the numerical error was small enough to
keep the qualitatively analysis of the physics still valid, in this work we
corrected the code so as to avoid the linear drifts coming from the
integration on sigma.

The amplitude and phase for the simulations are obtained in compressed format
from the Rochester catalog [26]. An interpolation is made with one dimensional
smoothing splines. Order k=5 splines are implemented in the interpolation.
Subsequently, the polarizations of gravitational $h^{lm}_{+},h^{lm}_{\times}$
wave strain can be obtained by the relation,

$rh^{lm}_{+}-irh^{lm}_{\times}=A_{lm}e^{i\phi_{lm}},$

where $r$ is a radial coordinate, $l$ and $m$ the modes of the gravitational
wave, $A_{lm}$ and $e^{i\phi_{lm}}$ the amplitude and phase of the
gravitational wave respectively. Asymptotically, the strain and shear are
directly related by the $\Psi_{4}$ scalar, implying

$\displaystyle h_{+}=-\sigma_{R},$ $\displaystyle h_{\times}=-\sigma_{I}.$

The code implements the calculation of the shear and then it transforms from
the spinor to the tensor harmonics basis using the transformation relations in
[27]. After that, the calculation for all the physical variables is
implemented. In particular, the evolution of eqs. 4, 1 and 2 can be carried
out for all times. Finally, the code implements functions to collect and
analyze the results obtained. Data analysis is shown in Sec. IV.

##  IV  Main results

In this section we first describe the general results obtained for the global
variables of an asymptotically flat space time containing a BH binary system.
We provide the distribution of the data and give an overview of the physics
obtained for all the simulations. Then, we dedicate a subsection for each
categories made of the Rochester repository: non-spinning Binaries, aligned-
spins binaries and precessing binaries.

We made a special analysis of the variables obtained by eqs. 1, 2, 3, 4, 5 and
6 and the correlations with initial parameters from the simulations, such as
initial total angular momentum $J_{in}$ and mass ratio $q$. It is important to
distinguish these two sets of parameters. The first set, the global variables,
are well-defined in our formalism and can be obtained at null infinity without
knowledge of the BH masses, spins and orbital angular momentum. For instance,
to obtain $J_{in}$ we only need the knowledge of the center of mass position
and the definition of intrinsic angular momentum.

The second set of parameters corresponds to all the initial conditions
contained in the Rochester repository, locally defined, containing the initial
spins, masses, and orbital angular momentum of the BH binaries. Knowledge of
this second set is very important to analyse the behaviour of BH coalescence,
and gives further understanding of this kind of isolated system.

Throughout the following subsections we will distinguish between six classes
of binaries: equal mass non-spinning (EM-NS), equal mass aligned (EM-A), equal
mass precessing (EM-P), non-equal mass non-spinning (NEM-NS), non-equal mass
aligned (NEM-A), and non-equal mass precessing (NEM-P).

###  IV.1 General results

Figure 1: Correlation between the final center of mass velocity, $V_{f}$ (in
$km/s$), and the initial total angular momentum, $J_{in}$.

The results for the final center of mass final velocity $V_{f}$ and initial
total angular momentum $J_{in}$ are presented in Fig. 1. The final velocity
$V_{f}$ is defined as,

$V^{i}_{f}=\dot{R}^{i}_{f},$ (7)

where $\dot{R}^{i}_{f}=\dot{R}^{i}(t_{f})$ is the derivative with respect to
Bondi time evaluated at the final time of the simulation and the index $i$
denotes the component of the vector.

The plot in Fig. 1 shows the final velocities for the 776 simulations in the
catalog. The vertical axis $V_{f}$ is the norm of the vector velocity defined
in Eq. (7). One can see that precessing binaries exhibit higher final
velocities. This observation could be relevant both at a level of predicting
the final velocities of the residual black hole or at the level of
observational data. If the final velocity of the residual black hole can be
observed (via Doppler effect of a surrounding coalescing gas), then receding
speeds higher than $\sim 750\ km/s$ implies that the original binaries had
precessing spins. One can even speculate that if the final speeds are higher
than $\sim 3000\ km/s$, then the initial binaries were in the class of EM-P.
It is also possible to recognize in Fig. 1, that within each class, there is a
value for the initial total angular momentum $J_{in}$ with the highest final
velocity $V_{f}$. The distribution of the simulations seems to exhibit an
asymmetric peaked-like distribution, whose peak value varies among each
subclass in the classification. On account of its symmetry, the EM-NS class
has vanishing final velocities (see formula (10) below). The reader should be
aware these assertions are rather approximate due to the bias in the RIT
catalog parameter space. However, the parameter space is vast enough to
visualize the claims made in the paragraph and in this sense phenomenological
formulae are fitted in the following sections with the idea of filling the
blanks in the initial parameter space. A special exception must be made with
precessing binaries, for which the final state is known to be chaotic with
respect to the initial parameters.

(a) Distribution of angular momentum change $\Delta S$ (in code units).

(b) Distribution of final velocity $V_{f}$ (in km/s)

(c) Distribution of total radiated energy $E_{rad}$

Figure 2: Bar plots of the distribution from main variables of interest
obtained after evolution in all the simulations.

We have solved the eqs. 1, 3, 2, 4, 5 and 6 and plotted the data distribution
for the variables $E_{rad}$,$V_{f}$ and $\Delta S$ in form of histogram
(Fig.2). Again, the variables of distributions in Fig.2 are also biased by the
selection of simulations in RIT catalog. All the distributions are normalized
as a probability density function and a gaussian function is fitted to get
representative values. For the final velocity in Fig.2(b), note that most of
the simulations will end with kicks less than $\sim 1000\ km/s$. It is our
purpose to investigate in the following sections which are the factors that
cause this large variation of final velocities. Something important to remark
is that final velocities are expressed in Bondi time used in our formalism.
For comparison purposes, we should divide our final velocities by a $\sqrt{2}$
factor. This factor accounts for the transformation from the Bondi time $u$ to
the standard time $t$ due to have defined $u=\frac{t-r}{\sqrt{2}}$ and not
$u=t-r$ as found in common literature. Consequently, the kick velocities in
terms of the coordinate time $t$ are lower than those showed in Fig. 2(b).

It is more convenient to plot the data distribution for a normalized radiated
energy instead of the raw energy change, as in Fig. 2(c). We define the total
energy radiated as

$E_{rad}=1-\frac{M_{f}}{M_{i}},$ (8)

where $M_{i}$ is the initial ADM mass of the BBH and $M_{f}$ the final Bondi
mass of the remnant obtained by Eq. (1). Note the difference between the final
$M_{f}$ in this work and the commonly used Christodoulou mass in the numerical
literature. The distribution of radiated energy for our set of solutions is
found between the range $\sim 0-13\%$. The histogram is normalized and the
mean of the radiated energy is $\mu=0.047$ and the dispersion $\sigma=0.021$.

Definition (8) is directly proportional to the peak energy loss
$\dot{M}_{max}$ for binary coalescence. The bigger the peak energy loss, the
bigger will also be the total energy radiated. In Fig. 3, we fitted a linear
function of the form $E_{rad}=a_{1}\dot{M}_{max}$ to the relation and found a
proportionality constant

$a_{1}=41.348\pm 0.194$ (9)

with an error of $0.47\%$. Moreover, the proportionality relation could be
made more precisely if fitted for each mass ratio as seen in Fig. 3. Having
this result in mind, figures present in this work with $E_{rad}$ variable will
be nearly the same as figures with $\dot{M}_{max}$ variable but with a change
of scale. It must be emphasized that this proportionality relation is a
feature of the BBH coalescence and will not be valid in general. We analyze
these variables for each subclass in the following sections more specifically.

Figure 3: Correlation between total radiated energy $E_{rad}$ and the peak
radiated energy $\dot{M}_{max}$ for all simulations. Dashed line shows fit to
data. Colorbar indicates the value of mass ratio $q=m_{1}/m_{2}$. Figure 4:
Correlation between the radiated energy $E_{rad}$ and the initial angular
momentum $J_{in}$ for all the simulations

To end this section, we also give an overview of the distribution of radiated
energy with respect to the initial total angular momentum $J_{in}$ in Fig.4.
$E_{rad}$ vs $J_{in}$ clearly shows a non-linear relationship. This feature
will be analysed for each group made in the next subsection and also in the
discussion. A model to understand the behaviour of the radiated energy with
respect to the initial total angular momentum is developed in Appendix A. It
is worth mentioning that this model has a post-Newtonian derivation. Thus, we
will contrast the fully relativistic numerical evolution of variables that
have asymptotic nature with a PN model describing its dynamics. It is rather
remarkable that with few phenomenological parameters that are adjusted to fit
the numerical data, the PN formula matches the numerical data.

###  IV.2 Non-spinning binaries

In this subsection we describe the group of simulations which have no initial
spins. This group of binaries has the simplest physics to describe since the
space of initial parameters has the lowest dimension.

As a first step, we calculate the correlations of the final velocity and the
maximum gravitational luminosity of the binary system with the mass ratio $q$.
The gravitational radiation is calculated by Eq.(1) and it is directly related
to the luminosity by integrating it into a spherical surface area far away
from the source. As the final luminosity is a quantity that can be usually
measured, it is relevant to know if a direct relation exists between the
radiation per unit time with the mass ratio of the BBH. Likewise, assuming the
final velocity of the remnant lack hole can be observed, the correlation with
$q$ also yields information about this particular coalescence.

Correlations between $V_{f}$, $q$ and $\dot{M}_{max}$ are outlined in Fig. 5.
We use the Fitchett’s recoil model for circular orbits proposed in [28] whose
dependence with q is

$V_{f}(q)=a\ \frac{q^{2}(1-q)}{(1+q)^{5}}.$ (10)

with $a$ a fitting constant. We compute the parameter $a$ by a least square
method. The best fit to our data has the coefficient $a=16317.78\pm 207.93\
km/s$. The derivation of the Eq. (10) is also given in Appendix A.

We found the maximum final velocity $V_{f}=309.43$ reached at the value
$q=0.4$. Null velocities are reached for limit cases when $q=0$ or $q=1$, i.e.
mass $m_{1}=0$ or $m_{1}=m_{2}$. This is consistent with the fact that when
$q=0$, the center of mass is in the center of the mass $m_{2}$; hence, it does
not move while the mass $m_{1}$ goes to zero. Similarly, when $q=1$, the
center of mass does not move either due to the symmetry of the problem and
coincides with half the distance between the binaries.

Even though limit cases of $q$ yield vanishing velocities, Fig. 5 shows that
while $q=0$ has no energy loss, $q=1$ allows the most radiating systems. Thus,
$q=1$ non-spinning binaries systems are most likely to be detectable by its
higher luminosity, whereas intermediate values for $q$ are most likely to be
detectable by their kick velocities.

Figure 5: Correlation between the final velocity $V_{f}$ and mass ratio $q$.
Colorbar indicates the maximum energy lost $\dot{M}_{max}$. Dashed line shows
the Fitchett model fit.

Another interesting relation we have explored is the radiated energy $E_{rad}$
vs the initial total angular momentum $J_{in}$. For the NS subgroup, this plot
yields Fig 6.

Figure 6: Correlation between the total radiation emitted $E_{rad}$ and the
initial total angular momentum $J_{in}$. Colorbar indicates different levels
of the initial mass ratio $q$. Letter $k$ indicates polynomial degree of the
fit.

Fig.6 is fitted with a quadratic curve of the form

$E_{rad}=b\ J_{in}^{2},$ (11)

with $b=0.050\pm\ 0.001$. The model fitted in Eq. (11) is motivated by the
derivation made in Appendix A.

Likewise, we can also give the correlation between $\dot{M}_{max}$ and $q$.
The plot is given in Fig. 7 and the phenomenological formula containing two
parameters is derived in Appendix A. The relation is given as,

$\dot{M}=-\frac{A}{1000}\frac{q^{2}}{(1+q)^{4}}\left(1+\frac{B}{216}\left(\frac{1-q}{1+q}\right)^{2}\right),$
(12)

with $A$ and $B$ to be obtained from the numerical data. Notice that the
formula has a PN background as can be checked in Appendix A. It follows from
the data or the formula that the luminosity has a minimum at $q=0$ and a
maximum at $q=1$, i.e., equal BBH masses yield maximum luminosity. The
coefficients found with least squares method are

$\displaystyle A$ $\displaystyle=22.44\pm 0.06\ (0.30\%)\ ,$ (13)
$\displaystyle B$ $\displaystyle=-159.00\pm 4.74\ (2.98\%)\ .$ (14)

Note also that relationship (12) is injective. Thus, given the coefficients in
eqs. 13 and 14, one can invert the formula (12) and determine which mass ratio
$q$ corresponds to a particular value of the measured luminosity, i.e., global
dynamical information is used to obtain local information of the BBH.

Figure 7: Correlation between the peak energy radiated $\dot{M}_{max}$ and
mass quotient $q$. Colorbar indicates values of energy radiated which are
proportional to the peak energy loss. Dashed line shows the model fit to data.

###  IV.3 Aligned-spins binaries

In this section, we analyse the aligned-spin class. The initial spins in this
kind of configurations are likely to remain in the z-direction as aligned-
spins binaries are known to be stable configurations [29]. These kinds of
binaries are more realistic ones and hence relevant for astrophysical
applications.

In Fig.8 we outline the dependence of the final center of mass velocity
$V_{f}$ with the initial mass ratio $q$. Figures (a) and (b) highlight
different initial configurations for magnitudes and spins directions. In
figure (c) a colorbar is used to show the highest luminosity released for each
configuration. Over all the 407 aligned simulations, the norm of $S_{2}$
sweeps values in [0,0.7] interval and the absolute value of $S_{1}$ in
[0.,0.25] interval.

Fig.8(a) shows the highest final velocity $V_{f}=771.78\ km/s$ for the
quotient $q=0.6628$ in the whole class of aligned-spins binaries. In the same
plot we find that top velocities can only be achieved for BBHs with anti-
aligned spins ($\angle S_{1}S_{2}=\pi$). This property has already been
reported several times [30, 31]. On the other hand, there seems to be no clear
picture of the influence between the spin alignment of $S_{1}$ and $S_{2}$ in
lower final velocities.

Furthermore, it can be seen in Fig.8(b) that when the initial orbital angular
momentum $L_{in}$ is anti-aligned with the initial spin $S_{2}$ ($\angle
L_{in}S_{2}=\pi$), higher final velocities are reached. In short, it seems
that the higher absolute value spins, in this case $S_{2}$, should be anti-
aligned with initial orbital angular momentum $L_{in}$ and the initial lower
spin, in this case $S_{1}$, to achieve the maximum velocities for the aligned
group.

The correlation between final velocities $V_{f}$ and maximum rate of mass loss
$\dot{M}_{max}$ is shown in Fig.9. A similar pattern discussed for NS binaries
IV.2 is also present. The maximum possible velocity versus q has a local (and
absolute) maximum at $q\sim 0.66$. Final velocities vanishes at $q=0$, reach a
peak value and then decrease at $q=1$. Moreover, energy loss also vanish at
$q=0$ which indicates that both the NS and A classes are very similar at this
part of the mass ratio range.

At the other side of the range, the equal mass pattern $q=1$ is quite
interesting. Like the NS case, the maximum energy loss and vanishing final
velocity also occurs at $q=1$. However, unlike the NS case, there are now non-
vanishing final velocities in the EM-A configuration, a big difference from
the NS case.

(a) Red and blue colors denote alignment of spin $S_{2}$ with respect to
orbital angular momentum $L_{in}$. Green colors represent binaries whose
$S_{2}$ alignment with respect to $L_{in}$ can not be established or
differentiated. Different markers are used to denote alignment type between
spins $S_{1}$ and $S_{2}$.

(b) Red and blue colors denote alignment of spin $S_{2}$ with respect to
orbital angular momentum $L_{in}$. Yellow color has $S_{2}=0$ and thus no
alignment can be establish with respect to $L_{in}$. Green dots are
indifferent as $S_{1}=S_{2}$. Markers clear out the bigger value of the spin
in simulations.

Figure 8: Correlation between final velocity $V_{f}$ and mass ratio $q$ for
aligned binaries. Figure 9: Correlation between $V_{f}$ and $q$ for aligned
binaries. Colorbar shows levels of peak energy loss $\dot{M}_{max}$. Dashed
line indicates border limit to final velocities (in $km/s$).

The relation of the energy loss is then analyzed with respect to the initial
total angular momentum $J_{in}$. In Fig. 10, the dependence between the
initial total angular momentum $J_{in}$ and total radiated energy $E_{rad}$ is
shown for five mass range ratios: $0<q\leq 0.2$, $0.2<q\leq 0.4$, $0.4<q\leq
0.6$, $0.6<q\leq 0.8$ and $0.8<q\leq 1$. We fit a quadratic polynomial model
for each range which could be useful for estimating energy on astrophysical
scenarios or also for future comparisons. We provided the coefficients of the
fit for each range of masses in Table 1.

Figure 10: Correlation between initial total angular momentum $J_{in}$ and total radiation emitted $E_{rad}$ for aligned binaries. Mass ratio intervals have been divided and selected to be studied and are illustrated with different colors $q$ | $a_{0}$ | $a_{1}$ | $a_{2}$
---|---|---|---
$0<q\leq 0.2$ | $0.037\pm 0.016$ | $-0.023\pm 0.017$ | $0.0095\pm 0.0031$
$0.2<q\leq 0.4$ | $0.032\pm 0.004$ | $-0.007\pm 0.005$ | $0.0138\pm 0.0013$
$0.4<q\leq 0.6$ | $0.093\pm 0.009$ | $-0.102\pm 0.015$ | $0.0553\pm 0.0065$
$0.6<q\leq 0.8$ | $0.12\pm 0.01$ | $-0.156\pm 0.019$ | $0.082\pm 0.009$
$0.8<q\leq 1$ | $0.149\pm 0.014$ | $-0.212\pm 0.029$ | $0.114\pm 0.014$
Table 1: Coefficients from fitted second-degree polynomials in Fig.10

To end the section of aligned-spins BBHs and for comparison purposes, we will
restrict our analysis to $q=1$ aligned BBHs and employ Reisswig et al. model
from Ref. [32]. The authors have shown that radiated energy via gravitational
waves from infinity of equal-mass binaries with aligned spins can be estimated
by a quadratic polynomial on the average initial spin
$\bar{\chi}=(\chi_{1}+\chi_{2})/2$, where $\chi_{1}$ and $\chi_{2}$ are the
projections of the initial spin in the $L_{in}$ direction. Therefore, the
model reads

$E_{rad}=a_{0}+a_{1}\bar{\chi}+a_{2}\bar{\chi}^{2}.$ (15)

Our fit to the data give us the following vector coefficients

$\vec{a}=\begin{pmatrix}a_{0}\\\ a_{1}\\\
a_{2}\end{pmatrix}=\begin{pmatrix}0.051\pm 0.001\\\ 0.040\pm 0.002\\\ 0.029\pm
0.003\end{pmatrix},$ (16)

whereas the corresponding coefficients in Ref. [32] yield

$\vec{p}=\begin{pmatrix}p_{0}\\\ p_{1}\\\
p_{2}\end{pmatrix}=\begin{pmatrix}0.036\pm 0.003\\\ 0.030\pm 0.006\\\ 0.02\pm
0.01\end{pmatrix}$ (17)

.

We see in Fig.11 there is a notable difference between the polynomials. The
difference could be attributed to the fact that the parameters for the used
simulations differ significantly and thus the final the $E_{rad}$.
Particularly, no comment is said about initial orbital angular momentum
$L_{in}$ in Ref. [32] which can generate different outcomes in the radiated
energy for the same values of the initial spins $\chi_{1}$,$\chi_{2}$.
However, it could be also that radiated energy predicted for aligned BBHs in
the framework of KQ formalism is slightly higher than radiated energy obtained
using the local quantities such as the Christodoulou mass in numerical
relativity.

The largest error is in 2nd-order coefficient and is $\sim 10\%$. We find the
maximum radiated energy evaluating at $a=1$ in the polynomial and we find the
value $E_{rad}(1)=12.2\%$. The limit is closer to $E_{rad}(1)=11.3\%$ reported
in [2] than $E_{rad}(1)=9.9\%$ reported in [32]. The maximum radiated energy
is also below the value $E_{rad}\sim 14\%$ found in Ref. [33] for head-on
collision of two highly boosted EM-NS black holes.

Figure 11: Correlation between initial total spin $\bar{\chi}$ and total
radiation emitted $E_{rad}$. Dashed grey line shows the fit to our data and
below its residuals. Dashed light blue line shows Reisswig polynomial for
comparison purposes.

###  IV.4 Precessing spins binaries

We provide in this section a deeper analysis on the parameters for the
precessing spins group. For this group of BBH, initial spins $S_{1}$ and
$S_{2}$ will, in general, not be aligned nor even along $\textbf{L}_{in}$
axis. Different behaviour with respect to the previous groups mentioned in
Secs.IV.2 and IV.3 is expected as non aligned configurations are known to
present chaotic behaviour [29].

We first study, as in previous sections, the dependence of the final velocity
$V_{f}$ with the initial mass ratio $q$. In Fig.12(b) one can appreciate that
the most notable difference in the behaviour of $V_{f}$ with respect to the
non-spinning or aligned BBHs is the achievement of highest velocities for
equal mass binaries. By contrast, $\dot{M}_{max}$ behaves equally as in other
classes, i.e., systems with $q=0$ have no energy emission as expected and
$q=1$ are the most energetic BBHs. The rate of energy loss increment as the
mass ratio does.

On the left side, Fig.12(a) shows $E_{rad}$ vs initial total angular momentum
$J_{in}$. We used the same model as with aligned and non-spinning classes
derived in Appendix A to fit the data, i.e., a 2nd degree polynomials for each
range of masses. The coefficients of the polynomials are listed in table 2.
Even though second degree polynomials may not be the best fit for the data, we
prioritize its simplicity for describing the ascending behaviour of each mass
range in the precessing group. Note also that range $0<q<0.2$ is missing in
plot 12(a) due to the lack of enough data to provide a representative curve in
that interval.

$q$ | $a_{2}$ | $a_{1}$ | $a_{0}$
---|---|---|---
$0.2<q\leq 0.4$ | $0.040\pm 0.009$ | $-0.020\pm 0.014$ | $0.017\pm 0.005$
$0.4<q\leq 0.6$ | $0.070\pm 0.013$ | $-0.071\pm 0.025$ | $0.046\pm 0.011$
$0.6<q\leq 0.8$ | $0.089\pm 0.024$ | $-0.097\pm 0.048$ | $0.059\pm 0.024$
$0.8<q\leq 1$ | $0.180\pm 0.038$ | $-0.271\pm 0.081$ | $0.143\pm 0.041$
Table 2: Coefficients of second-degree polynomials.

(a) Correlation between initial total angular momentum $J_{in}$ and total
radiation emitted $E_{rad}$.

(b) Correlation between final velocity $V_{f}$ (in $km/s$) with the initial
mass ratio $q$. Colorbar shows different levels of peak energy loss.

Figure 12: Detailed description of the many variable dependence for the
precessing binaries.

##  V Discussion and Conclusion

In this work, we have studied final velocities of the remnant black hole for
the non-spinning (NS), aligned (A) and precessing (P) groups. We have shown
that final velocities find their maximum value for different mass ratios in
each group: $V_{f}\approx 310\ km/s$ for $q=0.4$ in NS, $V_{f}=771\ km/s$ for
$q=0.66$ in A, and $V_{f}=6120.15\ km/s$ for $q=1$ in P. We have given the
final velocities in Bondi time. For comparison purposes, we should divide our
final velocities by a $\sqrt{2}$ factor, as explained in Sec.IV.1. Thus, the
maximum final velocity in the precessing group is $V_{f}=4327.59\ km/s$, which
is slightly lower than the final kick reported for this simulation in
Rochester metadata ($v_{kick}=4625.73\ km/s$).

We have found $V_{f}=219.20\ km/s$ at $q=0.4$ to be the maximum final velocity
for NS group. This value is consistent with the reported maximum value
$V_{f}=175\ km/s$ for the ratio $q\approx 0.36$ found in Ref. [3]. For the NS
group, we employed the Fitchett’s recoil model to fit the set of NS
simulations in the catalog. We obtained a fitting constant $a=16317.78\ km/s$
and associated error of $\sim 1.5\%$. The model can be used for posterior
simulations comparison or to model astrophysical BBHs with negligible
rotation. Likewise, we derived a formula to fit the correlation between
$\dot{E}_{max}$ and $q$ using a phenomenological dimensionless parameter $b$
which should be greater than $2$. The value obtained $b=2.1867$ is a
reasonable result.

For the A group, we found that the top final velocity is $V_{f}=545.73\ km/s$,
which is consistent with the maximum value in the metadata $v_{kick}=500\
km/s$ and of the same order found in Ref.[30], $V_{f}=448\ km/s$. Moreover,
our results also show that final velocities are maximal when spins are
antialigned. This result have been also reported in other works [30, 31].

We have also found that the highest attainable final speeds come from the EM-P
group. According to the Rochester repository, if the remnant black hole is
found with speeds higher than $800\ km/s$, then the BBH had precessing spins.
Moreover, if the speed of the remnant BH is higher than $3000\ km/s$ then the
inital binaries belong to the EM-P subclass.

We have also studied in depth the radiation for each class predicted by eq. 8.
For the NS group, we find that using the model derived in Appendix A, a
quadratic polynomial fits notably well the data with a coefficient $b=0.04976$
and error $\sim 1.74\%$. The simplicity of the non-spinning binaries allows
one to model the dependence of final variables with initial parameters and to
obtain formulas with smaller errors than in other classes of BBHs.

In the group of aligned spins binaries we made a fit for different ranges of
mass ratios and found higher error values than non-spinning binaries fit. The
higher errors could be attributed to the coefficients dependence of the spin.
Perhaps more complex models should be considered to reduce the coefficient
errors. Despite that, the coefficients calculated still provide a simple and
useful model to the relation between $E_{rad}$ and $J_{in}$. For the case of
$q=1$ aligned binaries we fit again the radiated energy but use the variable
$\bar{\chi}$ instead of $J_{in}$. The plot in terms of the variable
$\bar{\chi}$ allows us to compare the predicted radiated energy with other
works. We see the fit is consistent although barely higher than most reported
radiated energies. This could be a consequence of defining 8 with Bondi mass.
A more thorough analysis on the difference between radiated energy defined via
Eq.(1) and via Christodoulou mass should be made, and we left its study for
future work.

We found similar results for the precessing group. There are some ranges of
mass ratio where the fits have higher error in the linear coefficients (up to
$\sim 70\%$) and data is not enough for fitting low mass ratio binaries, i.e.,
for $0<q<0.2$.

We now make some general remarks related to this work:

  * •

Reisswig et al [32] have shown that the gravitational wave pattern for EM-NS
and EM-A are quite similar. Since the numerical evolution for EM-NS binaries
are much simpler than the EM-A subclass, they argue that one could simply use
the NS class to obtain the waveforms for both subgroups. On the other hand,
our results show that the final velocities for the EM-NS and EM-A sub classes
are completely different. Whereas the former class have vanishing recoil
velocities, EM-A cases do not. This is a distinctive difference between the
two subclasses and helps to distinguish between them. Thus, even though the
waveform templates are virtually the same for these sub classes, the numerical
evolution of the EM-A subclass yields valuable information of the binary
coalescence.

  * •

Figs.10 and 12(a) represent truly global results. The plot of $E_{rad}$ vs
$J_{in}$ is motivated by the following. Both variables can be obtained at null
infinity without knowledge of the BH masses, spins and orbital angular
momentum. To obtain $J_{in}$ we only need the knowledge of the center of mass
position and the definition of total angular momentum, whereas $E_{rad}$ is
computed with the use of the Bondi mass loss equation. This plot shows a clear
correlation between $E_{rad}$ and $J_{in}$. Using a quadratic fit we obtained
an empirical relationship for the coefficients. The results in this works show
that the initial total angular momentum $J_{in}$ is a relevant variable in the
parameter space to analyze the different numerical evolutions. It is fair to
ask why there should be a quadratic dependence on the initial total intrinsic
angular momentum of the BBH system. The answer lies in the formula for
$E_{rad}$, and it depends on quadratic gravitational radiation terms. At the
same time, the radiation depends on the total initial angular momentum and it
is conserved if we only keep up to quadratic terms in the formula for
$E_{rad}$ since $\dot{J}$ is quadratic in the shear. Thus, we should expect
this dependency in the BBH coalescence. Note the same results for $E_{rad}$
applies also for $\dot{M}_{max}$ as a linear relation between them was shown
in fig. 3.

  * •

The angular momentum loss $\Delta S$ only depends on the amount of angular
momentum radiated away, and it is computed with knowledge of the available
radiation data at null infinity. The angular momentum loss is still a main
reason of concern in our formalism. As previously shown in [24], there is a
discrepancy of the predicted intrinsic angular momentum with respect to the
predicted values in RIT catalog, which is a subject to be studied in future
works. Of course, a physically relevant definition of intrinsic angular
momentum is a difficult task but the one we provide appears to be free from
ambiguities. Despite that, we have shown throughout this paper that the
equations of motion for energy, linear momentum, center of mass position and
final velocity [eqs. 1, 2, 3 and 5] employed in the framework of KQ formalism
are consistent with the general literature.

###### Acknowledgements.

This research has been supported by grants from CONICET and the Agencia
Nacional de Promoción de la Investigación, el Desarrollo Tecnológico y la
Innovación.

## References

  * Healy and Lousto [2020] J. Healy and C. O. Lousto, Third rit binary black hole simulations catalog, Physical Review D 102, 104018 (2020).
  * Hemberger  _et al._ [2013] D. A. Hemberger, G. Lovelace, T. J. Loredo, L. E. Kidder, M. A. Scheel, B. Szilágyi, N. W. Taylor, and S. A. Teukolsky, Final spin and radiated energy in numerical simulations of binary black holes with equal masses and equal, aligned or antialigned spins, Physical Review D 88, 064014 (2013).
  * González  _et al._ [2007] J. A. González, M. Hannam, U. Sperhake, B. Bruegmann, and S. Husa, Supermassive recoil velocities for binary black-hole mergers with antialigned spins, Physical Review Letters 98, 231101 (2007).
  * Campanelli  _et al._ [2007a] M. Campanelli, C. Lousto, Y. Zlochower, and D. Merritt, Large merger recoils and spin flips from generic black hole binaries, The Astrophysical Journal Letters 659, L5 (2007a).
  * Campanelli  _et al._ [2007b] M. Campanelli, C. O. Lousto, Y. Zlochower, and D. Merritt, Maximum gravitational recoil, Physical Review Letters 98, 231102 (2007b).
  * Healy  _et al._ [2014] J. Healy, C. O. Lousto, and Y. Zlochower, Remnant mass, spin, and recoil from spin aligned black-hole binaries, Phys. Rev. D90, 104004 (2014), arXiv:1406.7295 [gr-qc] .
  * Healy and Lousto [2018] J. Healy and C. O. Lousto, Hangup effect in unequal mass binary black hole mergers and further studies of their gravitational radiation and remnant properties, Physical Review D 97, 084002 (2018).
  * Campanelli  _et al._ [2006] M. Campanelli, C. O. Lousto, and Y. Zlochower, Spinning-black-hole binaries: The orbital hang-up, Physical Review D 74, 041501 (2006).
  * Zlochower and Lousto [2015] Y. Zlochower and C. O. Lousto, Modeling the remnant mass, spin, and recoil from unequal-mass, precessing black-hole binaries: The intermediate mass ratio regime, Phys. Rev. D 92, 024022 (2015).
  * Rezzolla  _et al._ [2008] L. Rezzolla, E. Barausse, E. N. Dorband, D. Pollney, C. Reisswig, J. Seiler, and S. Husa, Final spin from the coalescence of two black holes, Physical Review D 78, 044002 (2008).
  * Hofmann  _et al._ [2016] F. Hofmann, E. Barausse, and L. Rezzolla, The final spin from binary black holes in quasi-circular orbits, The Astrophysical Journal Letters 825, L19 (2016).
  * Lousto and Zlochower [2014] C. O. Lousto and Y. Zlochower, Black hole binary remnant mass and spin: A new phenomenological formula, Physical Review D 89, 104052 (2014).
  * Lousto  _et al._ [2010] C. O. Lousto, M. Campanelli, Y. Zlochower, and H. Nakano, Remnant masses, spins and recoils from the merger of generic black hole binaries, Classical and Quantum Gravity 27, 114006 (2010).
  * Blanchet [2014] L. Blanchet, Gravitational radiation from post-newtonian sources and inspiralling compact binaries, Living reviews in relativity 17, 1 (2014).
  * Kozameh and Quiroga [2016] C. N. Kozameh and G. D. Quiroga, Center of mass and spin for isolated sources of gravitational radiation, Physical Review D 93, 064050 (2016).
  * Winicour [1980] J. Winicour, Angular momentum in general relativity, in  _General relativity and gravitation_ , Vol. 2, edited by A. Held (Plenum Press, New York, 1980).
  * Adamo  _et al._ [2012] T. M. Adamo, E. T. Newman, and C. Kozameh, Null geodesic congruences, asymptotically-flat spacetimes and their physical interpretation, Living reviews in relativity 15, 1 (2012).
  * Kozameh  _et al._ [2018] C. N. Kozameh, J. I. Nieva, and G. D. Quiroga, Spin and center of mass comparison between the post-newtonian approach and the asymptotic formulation, Physical Review D 98, 064005 (2018).
  * Newman and Unti [1963] E. T. Newman and T. Unti, A class of null flat-space coordinate systems, Journal of Mathematical Physics 4, 1467 (1963).
  * Kozameh  _et al._ [2020] C. Kozameh, J. Nieva, and G. Quiroga, Relativistic center of mass in general relativity, Physical Review D 101, 024028 (2020).
  * Thorne [1980] K. S. Thorne, Multipole expansions of gravitational radiation, Reviews of Modern Physics 52, 299 (1980).
  * Newman and Silva-Ortigoza [2005] E. T. Newman and G. Silva-Ortigoza, Tensorial spin-s harmonics, Classical and Quantum Gravity 23, 497 (2005).
  * [23] E. A. Tassone, pykqevolution, https://github.com/Emmatassone/pyKQEvolution.
  * Tassone  _et al._ [2021] E. A. Tassone, P. A. Mandrilli, C. N. Kozameh, G. D. Quiroga, and J. I. Nieva, Numerical evolution of the center of mass and angular momentum for binary black holes, Phys. Rev. D 104, 084038 (2021).
  * Reisswig and Pollney [2011] C. Reisswig and D. Pollney, Notes on the integration of numerical relativity waveforms, Classical and Quantum Gravity 28, 195015 (2011).
  * [26] R. numerical realtivity group, Catalog of numerical simulations, https://ccrg.rit.edu/numerical-simulations.
  * Mandrilli  _et al._ [2020] P. Mandrilli, J. Nieva, and G. Quiroga, Correspondence between tensorial spin-s and spin-weighted spherical harmonics, General Relativity and Gravitation 52, 1 (2020).
  * Fitchett [1983] M. Fitchett, The influence of gravitational wave momentum losses on the centre of mass motion of a newtonian binary system, Monthly Notices of the Royal Astronomical Society 203, 1049 (1983).
  * Kozameh and Quiroga [2020] C. Kozameh and G. Quiroga, Spin stability of a binary black hole coalescence, Physical Review D 102, 064015 (2020).
  * Pollney  _et al._ [2007] D. Pollney, C. Reisswig, L. Rezzolla, B. Szilágyi, M. Ansorg, B. Deris, P. Diener, E. N. Dorband, M. Koppitz, A. Nagar, _et al._ , Recoil velocities from equal-mass binary black-hole mergers: a systematic investigation of spin-orbit aligned configurations, Physical Review D 76, 124002 (2007).
  * Koppitz  _et al._ [2007] M. Koppitz, D. Pollney, C. Reisswig, L. Rezzolla, J. Thornburg, P. Diener, and E. Schnetter, Recoil velocities from equal-mass binary-black-hole mergers, Physical Review Letters 99, 041102 (2007).
  * Reisswig  _et al._ [2009] C. Reisswig, S. Husa, L. Rezzolla, E. N. Dorband, D. Pollney, and J. Seiler, Gravitational-wave detectability of equal-mass black-hole binaries with aligned spins, Physical Review D 80, 124026 (2009).
  * Sperhake  _et al._ [2008] U. Sperhake, V. Cardoso, F. Pretorius, E. Berti, and J. A. Gonzalez, High-energy collision of two black holes, Physical review letters 101, 161101 (2008).
  * Kidder [1995] L. E. Kidder, Coalescing binary systems of compact objects to (post) 5/2-newtonian order. v. spin effects, Physical Review D 52, 821 (1995).

##  Appendix A Relation between $\dot{E}$, $J_{in}$ and $q$

Starting from the definition of the mass quadrupole moments, [34] or [14], we
have the following expressions:

$\displaystyle I^{ij}_{N}$ $\displaystyle=\eta mx^{<i}x^{j>}$ (18)
$\displaystyle J^{ij}_{N}$ $\displaystyle=-\frac{\delta
m}{m}x^{<i}L^{j>}+\frac{3}{2}\eta mx^{<i}(m_{2}\chi_{2}-m_{1}\chi_{1})^{j>},$
(19)

with the mass parameters given as $m=m_{1}+m_{2}$, $\delta
m=m_{1}-m_{2}$,$L^{i}$ the orbital angular momentum, and
$\eta=m_{1}m_{2}/m^{2}$ the symmetric mass ratio. For simplicity we consider
the non spinning, equal mass case. Thus, $\delta m=0$, $m=1$, and
$\eta=\frac{1}{4}$.With these mass parameters, the quadrupole radiative
moments are

$\displaystyle I^{ij}_{N}$ $\displaystyle=\eta
m\big{[}x^{i}x^{j}-\frac{1}{3}\delta_{ij}x^{2}\big{]}$ (20) $\displaystyle
J^{ij}_{N}$ $\displaystyle=0.$ (21)

Explicitly, the non-zero radiative moments in the N approximation are given
by,

$\displaystyle I_{N}^{zz}$ $\displaystyle=$ $\displaystyle-\frac{1}{12}mr^{2}$
(22) $\displaystyle I_{N}^{xx}$ $\displaystyle=$
$\displaystyle\frac{1}{24}mr^{2}[1+3\cos(2\omega t)]$ (23) $\displaystyle
I_{N}^{yy}$ $\displaystyle=$ $\displaystyle\frac{1}{24}mr^{2}[1-3\cos(2\omega
t)]$ (24) $\displaystyle I_{N}^{xy}$ $\displaystyle=$ $\displaystyle
I_{N}^{yx}=\frac{1}{8}mr^{2}[\sin(2\omega t)]$ (25)

Since the orbital angular momentum $L=\eta mr^{2}\omega$ is conserved in the N
approximation we obtain

$\displaystyle I_{N}^{(3)zz}$ $\displaystyle=$ $\displaystyle 0$ (26)
$\displaystyle I_{N}^{(3)xx}$ $\displaystyle=$ $\displaystyle
4L_{in}\omega^{2}[\sin(2\omega t)]$ (27) $\displaystyle I_{N}^{(3)yy}$
$\displaystyle=$ $\displaystyle-4L_{in}\omega^{2}[\sin(2\omega t)]$ (28)
$\displaystyle I_{N}^{(3)xy}$ $\displaystyle=$ $\displaystyle
I_{N}^{(3)yx}=-4L_{in}\omega^{2}[\cos(2\omega t)],$ (29)

where the symbol (3) denotes three time derivatives. Using Eq. (1) and taking
an average over a period, we have

$\displaystyle\dot{M}$
$\displaystyle=-\frac{1}{10\sqrt{2}}\dot{\sigma}_{R}^{ij}\dot{\sigma}_{R}^{ij}$
(30) $\displaystyle=-\frac{32}{10\sqrt{2}}L_{in}^{2}\omega^{4}.$ (31)

Thus,

$\dot{M}_{max}=-\frac{32}{10\sqrt{2}}\omega_{max}^{4}L_{in}^{2}.$ (32)

If the masses are not equal, the resulting formula for the non spinning
subclass is given by

$\dot{M}_{max}=-\frac{16}{10\sqrt{2}}\omega_{max}^{4}L_{in}^{2}\left(1+\frac{1}{36}(\frac{1-q}{1+q})^{2}(M_{in}\omega_{max})^{\frac{2}{3}}\right).$
(33)

The highest allowed value for $\omega$ in the PN approximation is $0.05$.
However, The resulting numerical factor of $4\times 10^{-5}$ in the above
equation should not be compared with the one obtained from the numerical
evolution since we are in the PN approximation and the binaries are not even
close to the coalescing time. However, it is important to see the quadratic
dependence of the maximum luminosity on the initial value of the total angular
momentum of the system. If we want to apply this model to the other spinning
classes (A or P), we use $L^{i}=J^{i}-(S_{1}+S_{2})^{i}$, which gives us a
second order polynomial with non vanishing zeroth and first order
coefficients. Although the coefficients will depend on the initial spins
($S_{1},S_{2}$), the model will be still useful to give a representative value
of the dependence.

One can also write the maximum luminosity in terms of the ratio $q$ by writing
the orbital angular momentum as

$L=\eta mr^{2}\omega=\frac{q}{(1+q)^{2}}mr^{2}\omega.$ (34)

Inserting this equation into (33) and using the circular orbit formula

$\omega^{2}r^{2}=\frac{m}{r}$

one obtains

$\dot{M}=-\frac{16}{10\sqrt{2}}\frac{q^{2}}{(1+q)^{4}}\left(\frac{m}{r}\right)^{5}\left(1+\frac{1}{36}\frac{m}{r}\left(\frac{1-q}{1+q}\right)^{2}\right).$
(35)

The last stable circular orbit before coalescence is given by

$r_{stable}=6m,$

Thus, the luminosity at that particular orbit is given by

$\dot{M}=-\frac{16}{10\sqrt{2}}\frac{q^{2}}{(1+q)^{4}}\left(\frac{1}{6}\right)^{5}\left(1+\frac{1}{216}\left(\frac{1-q}{1+q}\right)^{2}\right).$
(36)

We now assume that the final details of the coalescence do not change the
parameter dependence given by eq. (36). However, the huge amount of energy
must be taken into account by giving two phenomenological constants, one for
each contribution of the mass and magnetic quadrupole terms. We thus assume
that

$\dot{M}=-\frac{A}{1000}\frac{q^{2}}{(1+q)^{4}}\left(1+\frac{B}{216}\left(\frac{1-q}{1+q}\right)^{2}\right),$
(37)

and the values of $A$ and $B$ are found by fitting the formula with the
numerical data.

One can also make a phenomenological correlation between the final value of
the velocity $V_{f}$ o momentum $P_{f}$ and $q$. Note that Eq. (2) depends on
the product of the real and imaginary part of $\sigma$. Thus, an analogous
derivation to the one performed for $\dot{M}$ yields

$V_{f}=a\eta^{2}\frac{\delta m}{m}=a\frac{q^{2}}{(1+q)^{4}}\frac{1-q}{1+q},$
(38)

where the parameter $a$ can be obtained from the numerical plot. This last
result has been previously derived in Ref. [28].