Numerical computation of quasinormal modes in the first-order approach to black hole perturbations in modified gravity

We present a novel approach to the numerical computation of quasi-normal modes, based on the first-order (in radial derivative) formulation of the equations of motion and using a matrix version of the continued fraction method. This numerical method is particularly suited to the study of static black holes in modified gravity, where the traditional second-order, Schr\"odinger-like, form of the equations of motion is not always available. Our approach relies on the knowledge of the asymptotic behaviours of the perturbations near the black hole horizon and at spatial infinity, which can be obtained via the systematic algorithm that we have proposed recently. In this work, we first present our method for the perturbations of a Schwarzschild black hole and show that we recover the well-know frequencies of the QNMs to a very high precision. We then apply our method to the axial perturbations of an exact black hole solution in a particular scalar-tensor theory of gravity. We also cross-check the obtained QNM frequencies with other numerical methods.