  Linear perturbations of spherically symmetric spacetimes in general
relativity are described by radial wave equations, with potentials that depend
on the spin of the perturbing field. In previous work we studied the
quasinormal mode spectrum of spacetimes for which the radial potentials are
slightly modified from their general relativistic form, writing generic small
modifications as a power-series expansion in the radial coordinate. We assumed
that the perturbations in the quasinormal frequencies are linear in some
perturbative parameter, and that there is no coupling between the perturbation
equations. In general, matter fields and modifications to the gravitational
field equations lead to coupled wave equations. Here we extend our previous
analysis in two important ways: we study second-order corrections in the
perturbative parameter, and we address the more complex (and realistic) case of
coupled wave equations. We highlight the special nature of coupling-induced
corrections when two of the wave equations have degenerate spectra, and we
provide a ready-to-use recipe to compute quasinormal modes. We illustrate the
power of our parametrization by applying it to various examples, including
dynamical Chern-Simons gravity, Horndeski gravity and an effective field
theory-inspired model.
