  We discuss a general method to study linear perturbations of slowly rotating
black holes which is valid for any perturbation field, and particularly
advantageous when the field equations are not separable. As an illustration of
the method we investigate massive vector (Proca) perturbations in the Kerr
metric, which do not appear to be separable in the standard Teukolsky
formalism. Working in a perturbative scheme, we discuss two important effects
induced by rotation: a Zeeman-like shift of nonaxisymmetric quasinormal modes
and bound states with different azimuthal number m, and the coupling between
axial and polar modes with different multipolar index l. We explicitly compute
the perturbation equations up to second order in rotation, but in principle the
method can be extended to any order. Working at first order in rotation we show
that polar and axial Proca modes can be computed by solving two decoupled sets
of equations, and we derive a single master equation describing axial
perturbations of spin s=0 and s=+-1. By extending the calculation to second
order we can study the superradiant regime of Proca perturbations in a
self-consistent way. For the first time we show that Proca fields around Kerr
black holes exhibit a superradiant instability, which is significantly stronger
than for massive scalar fields. Because of this instability, astrophysical
observations of spinning black holes provide the tightest upper limit on the
mass of the photon: mv<4x10^-20 eV under our most conservative assumptions.
Spin measurements for the largest black holes could reduce this bound to
mv<10^-22 eV or lower.
