  We present a multipolar analysis of the gravitational recoil computed in
recent numerical simulations of binary black hole (BH) coalescence, for both
unequal masses and non-zero, non-precessing spins. We show that multipole
moments up to and including l=4 are sufficient to accurately reproduce the
final recoil velocity (within ~2%) and that only a few dominant modes
contribute significantly to it (within ~5%). We describe how the relative
amplitudes, and more importantly, the relative phases, of these few modes
control the way in which the recoil builds up throughout the inspiral, merger,
and ringdown phases. We also find that the numerical results can be reproduced
by an ``effective Newtonian'' formula for the multipole moments obtained by
replacing the radial separation in the Newtonian formulae with an effective
radius computed from the numerical data. Beyond the merger, the numerical
results are reproduced by a superposition of three Kerr quasi-normal modes
(QNMs). Analytic formulae, obtained by expressing the multipole moments in
terms of the fundamental QNMs of a Kerr BH, are able to explain the onset and
amount of ``anti-kick'' for each of the simulations. Lastly, we apply this
multipolar analysis to help explain the remarkable difference between the
amplitudes of planar and non-planar kicks for equal-mass spinning black holes.
