  We develop a formalism to treat higher order (nonlinear) metric perturbations
of the Kerr spacetime in a Teukolsky framework. We first show that solutions to
the linearized Einstein equation with nonvanishing stress tensor can be
decomposed into a pure gauge part plus a zero mode (infinitesimal perturbation
of the mass and spin) plus a perturbation arising from a certain scalar
("Debye-Hertz") potential, plus a so-called "corrector tensor." The scalar
potential is a solution to the spin $-2$ Teukolsky equation with a source. This
source, as well as the tetrad components of the corrector tensor, are obtained
by solving certain decoupled ordinary differential equations involving the
stress tensor. As we show, solving these ordinary differential equations
reduces simply to integrations in the coordinate $r$ in outgoing Kerr-Newman
coordinates, so in this sense, the problem is reduced to the Teukolsky equation
with source, which can be treated by a separation of variables ansatz. Since
higher order perturbations are subject to a linearized Einstein equation with a
stress tensor obtained from the lower order perturbations, our method also
applies iteratively to the higher order metric perturbations, and could thus be
used to analyze the nonlinear coupling of perturbations in the near-extremal
Kerr spacetime, where weakly turbulent behavior has been conjectured to occur.
Our method could also be applied to the study of perturbations generated by a
pointlike body traveling on a timelike geodesic in Kerr, which is relevant to
the extreme mass ratio inspiral problem.
