  We study the relativistic orbit of binary black holes in systems with small
mass ratio. The trajectory of the smaller object (another black hole or a
neutron star), represented as a particle, is determined by the geodesic
equation on the perturbed massive black hole spacetime. The particle itself
generates the gravitational perturbations leading to a problem that needs
regularization. Here we study perturbations around a Schwarzschild black hole
using Moncrief's gauge invariant formalism. We decompose the perturbations into
$\ell-$multipoles to show that all $\ell-$metric coefficients are $C^0$ at the
location of the particle. Summing over $\ell$, to reconstruct the full metric,
gives a formally divergent result. We succeed in bringing this sum to a
generalized Riemann's $\zeta-$function regularization scheme and show that this
is tantamount to subtract the $\ell\to\infty$ piece to each multipole. We
explicitly carry out this regularization and numerically compute the first
order geodesics. Application of this method to general orbits around rotating
black holes would generate accurate templates for gravitational wave laser
interferometric detectors.
