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 (ℓ,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.