  The gravitational wave strain emitted by a perturbed black hole (BH) ringing
down is typically modeled analytically using first-order BH perturbation
theory. In this Letter we show that second-order effects are necessary for
modeling ringdowns from BH merger simulations. Focusing on the strain's
$(\ell,m)=(4,4)$ angular harmonic, we show the presence of a quadratic effect
across a range of binary BH mass ratios that agrees with theoretical
expectations. We find that the quadratic $(4,4)$ mode's amplitude exhibits
quadratic scaling with the fundamental $(2,2)$ mode -- its parent mode. The
nonlinear mode's amplitude is comparable to or even larger than that of the
linear $(4,4)$ mode. Therefore, correctly modeling the ringdown of higher
harmonics -- improving mode mismatches by up to 2 orders of magnitude --
requires the inclusion of nonlinear effects.
