Until recently, all sorts of relativistic binary systems have been studied only theoretically and on the 14 September 2015 a team of LIGO and Virgo collaborators announced their first detection of a gravitational wave signal from a binary black hole system of about 36 and 29 solar masses. This announcement reaffirmed the predictions of the existance of gravitational waves as predicted by GR and most importantly the affirmation that indeed binary relativistic systems do exist in nature. LIGO, Virgo and all other promising gravitational wave detectors will thus provide with the means to be able to detect all kinds of relativistic binary systems with all sorts of physical properties  -  .
More so, the works in this direction have been in the realm of numerical relativity with a special focus of relativistic two body problems i.e. Black hole-black hole binary, black hole-neutron star binary, and neutron star-neutron star binary. The case of a three-body problem as it is the case here has never been studied before either numerically or analytically even though there is a much compelling possibility that in very dense cluster of galaxies these kinds of systems could in fact be found in the near future. As it was the case at the beginning of the research work in relativistic binary systems in the past decades, it is also very likely that there will be arising technical and computational issues for the gravitating three body problem in full numerical relativity. The analytic method used in this paper has been used by the author in  to validate other analytical results obtained by   for a binary system.
In the setting of this paper, we use the analytical results by Bishop  to study analytically the relativistic triple system consisting of two point-particles in quasi-orbit around a static Schwarzschild black hole. In practice the particles could be either two black holes, two neutron stars or in another setting a combination of both. Our objective thus is to determine the amount of the emitted gravitational radiation by the system at
in Bondi-Sachs formalism. This paper is structured as follows: Section 2 gives the background material. Section 3 defines the physical problem to be studied. Section 4 calculates the emitted gravitational radiation at
The Bondi-Sachs formalism uses coordinates
based upon a family of outgoing null hypersurfaces. We label these hypersurfaces by
, null rays by
, and the surface area coordinate by r. In this coordinates system the Bondi-Sachs metric   takes the form
being a unit sphere metric, U is the spin-weighted field given by
. For a Schwarzschild space-time,
. We define the complex quantity J by
For the Schwarzschild space-time, we have J and U being zero and thus they can be regarded as a measure of the deviation from spherical symmetry, and in addition, they contain all the dynamic content of the gravitational field in the linearized regime  . Usually we can describe this space-time by
, or by
For spherical harmonics we use
as basis functions as follows 
will be omitted in the case
are orthonormal and real. We assume the following ansatz
where r0 is the position of the matter shell, and σ the complex frequency mode which is physical damped and which further means that
. In the Bondi frame, the field equations splits into;
• the hypersurface equations and the evolution equations given by
• and the constraint equations for off the matter shell in the case of vacuum given by
Ref.  got the following second order differential equation when solving the above systems of ordinary differential equations for the Schwarzschild background;
, x is the compactification factor in this language. Bishop et al.  solved Equation (12) numerically and obtained interesting quasi-normal modes results of a Schwarzschild white hole. However in this paper, we are going to solved it for a different problem since we can apply the same physical settings in the Bondi-frame to model our problem with σ having a different physical meaning as we shall see later.
2.2. An analytic Algorithm for Calculating the Gravitational News
We shall use the following algorithm to calculate the gravitational radiation from the system.
• First we use Equation (12) and the constraints Equations (9)-(11) to get the junction conditions for the Bondi-Sachs matric variables U, ω and J at the boundary i.e. shell,
• Second we test if
, and ω are smooth across the boundary and if this is true, we then
• Calculate the News function at
3. The Problem
We consider a system consisting of two point-particles with equal mass m in quasi-orbit around a stationary Schwarzschild black hole with mass M situated at the center of mass r of the particles when l is 2. We take the orbital radius to be at r0 which means that the distance between the particles is 2r0. We take the initial position of particle 1 to be at r0 with θ and ϕ given by π/2 and νu respectively, ν is the orbital frequency and u the orbital period of the particles. We also take the initial position of a particle 2 to be at r0 with θ and ϕ given by π/2 and νu + π respectively. This imply that the rotation in the following figure is in the yz plane. The initial positions of the objects on the figure should not be confused with the actual initial positions just outline which in actual sense should be along the y axis with the particle 1 on the right and the particle 2 on the left.
The dynamics of this problem is governed by Equation (12) and for our numerical calculation purposes we shall use its Ricatti form 
where v is the orbital period of the system.
4. The Emitted Gravitational Radiation
4.1. The Linear Expansion of the Light Rays From the System to
We model the problem as follows, we start by applying Equation (5) with
where the matter density ρ in the background space-time is given by
Inside the particles orbital radius
and outside the particles orbital radius
Now integrating with respect to r we get
By multiplying Equation (18) with
and integrating over the sphere it simplifies to
From Equation (20), for
we the gravitational radiation otherwise we don’t, and that
are generally non-zero for even l and
. We now consider the case
and we note that
We note that
mode does not vary in time and hence it does not contain the emitted gravitational radiation. Thus we are only interested in
modes. We use the following normalized spherical harmonics
and the fact that
Thus from Equation (20)
and then finally we write
, Equation (32) then becomes
4.2. The Gravitational Radiation
We assume that the orbit is at the innermost stable circular orbit (ISCO), so that
. We then found the change in the Schwarzschild coordinate time
for one complete revolution of 92.3436 from which we found the orbital frequency ν of 0.0680.
To now find the numerical solutions to continue Equation (13) we make the spatial coordinate transformation of
which then imply that the ISCO is now at
. The numerical computations are done in the domains
with numerical solutions
respectively. We start the calculation with the transformed Equation (12) given by
are the Bondi metric functions, and
are the values of the expansion of the light rays β given by Equation (32) in the exterior and interior domains respectively. Bishop  has indicated that the derivatives of J should not be worked out numerically, but should be worked out analytically in terms of
and v from Equation (13) with
We define the general solutions for
outside and inside the orbital radius respectively as
where c4, c1, c2, c9, c6 and c7 are constants to be determined numerically. The functions
are analytic near
and therefore can be Taylor expand as
which then results in Equations (36) and (39) being analytic near
. We used Matlab ode45 solver to find numerical solutions of the above derivatives in Equations (38) and (39). We used stringent numerical conditions to get the results to about seven significant figures with RelTol of 10−12, AbsTol of 10−12, and the MaxStep of
and the results we found to be
We have tested for the consistency of the above results by using other Matlab solvers; ode23 and ode15s (which uses the Gears method i.e. backward differentiation formulas) and also observed the accuracy of about 15 significant figures. We went further with the test using ode23t which uses the trapezoidal rule, ode23s which is a modified Rosenbrock formula of order 2, and ode23tb which is an implicit Runge Kutta as opposed to ode45 and ode23 and found the consistency of about 8 significant figures and as opposed to 15 significant figures which is also accurate enough. This illustrate how accurate and valid the results are. These results are very crucial in obtaining the emitted gravitational radiation and hence determining the extent of their convergence is of most paramount importance.
From the hypersurface equation Equation (7) rewritten as
we are able to the Bondi metric function
. But to find the solution the integration should be done analytically where possible. We only need a solution which is valid in a neighborhood of
. Henceforth, it is convenient to make the coordinate transformation
. Equation (46) can further be rewritten as
. The constraints equations Equations (9), (10), and (11) now simplifies to
which we then apply in the domains
. Since these constraints are not completely analytic, this means that we should only evaluate them at the ISCO. We use them among others to eliminate the constants c1, c2, c6, and c7. We now assume that we end up with the solutions
Thus, from the constraints
, and the hypersurface Equation (47), we found the metric variables
. From which the expressions of the constants c9, c7, c5, and c10, were found.
We now impose the Bondi gauge conditions:
which means that for large r,
imply that the coordinate time is the same as proper time and that the regularity at
. We also impose the following junction conditions at
From the junction conditions, we were able to find the exact numerical values of the constants c1, c2, and c6 at
. The exact numerical values of the constants c9, c7, c5, and c10 were then found by substituting the values of c1, c2, and c6 back into their expressions. From here we were then able to plot the graphs of the Bondi metric functions
as observed in the following graphs.
Physically the metric functions J and U have the smooth asymptotic expansion characteristic through out the entire computational domain and this property is confirmed in Figure 1 and Figure 2. The metric function ω do not have this physical property as can be confirmed in Figure 3 but this function is crucial in the calculation procedure of the gravitation radiation in the entire domain. Physically the function J in the only one that have the time derivative and thus carries the gravitational radiation information to calculated at
and that all the other Bondi metric functions are intergrated radially from Γ to
. The above results indicate that the junction conditions at
Figure 1. The graph of
in the entire domain, and
also in the entire for the Schwarzschild space-time.
Figure 2. The graph of
for the Schwarzschild space-time.
Figure 3. The graph of
for the Schwarzschild space-time.
implemented correctly and that our numerical methods and the analytical algorithms we implemented to calculating the gravitational radiation worked properly as intended.
Then finally, since we are in the Bondi gauge, we found the gravitational news to be
which then further simplify to
with the Bondi mass loss of 0.0114 m2. The author in  has done a similar work for a single point particle in close orbit around a Schwarzschild black hole in the Bondi-frame and obtained the Bondi mass loss of 0.00089897 m2. He succeeded in validating the results by comparing it with that of the 5.5 PN formalism by Poisson  and Sasaki et al.  for the same problem. Thus the methods used in this article open up the possibilities in numerical relativity to be able to study analytically the gravitational radiation emitted by a systems consisting of one black hole and two equal orbiting black holes/neutron stars. With further improvement, the method can be develop to look at two unequal orbiting black holes or neutron stars or a combination of both with efficiency and accuracy as demonstrated in  for single orbiting black hole/neutron star.
I would like to thank Professor Nigel Bishop for suggestions and validation of the numerical methods and results and on improving the manuscript. I would also like to thank the National Research Foundation of South Africa under GUN 2053724 for financial support.
1) The constraints computed at r0
2) The Bondi metric variables computed at r0
The computed constants