Numerical integration of the Teukolsky Equation in the time domain

We present a fourth order convergent (2+1) numerical code to solve the Teukolsky equation in the time domain. Our approach is to rewrite the Teukolsky equation as a system of first order differential equations. In this way we get a system that has the form of an advection equation. This is used in combination with a series expansion of the solution in powers of time. To obtain a fourth order scheme we kept terms up to fourth derivative in time and use the advection-like system of differential equations to substitute the temporal derivatives by spatial derivatives. A local stability study leads to a Courant factor of 1.5 for the nonrotating case. This scheme is used to evolve gravitational perturbations in Schwarzschild and Kerr backgrounds. Our numerical method proved to be fourth order convergent in r* and theta directions. The correct power-law tail, ~1/t^{2\ell+3}, for general initial data, and ~1/t^{2\ell+4}, for time symmetric data, was found in the simulations where the duration in time of the tail phase was long enough. We verified that it is crucial to resolve accurately the angular dependence of the mode at late times in order to obtain these values of the exponents in the power-law decay. In other cases, when the decay was too fast and round-off error was reached before a tail was developed, the quasinormal modes frequencies provided a test to determine the validity of our code.