We formulate and calculate the second order quasi-normal modes (QNMs) of a Schwarzschild black hole (BH). Gravitational wave (GW) from a distorted BH, so called ringdown, is well understood as QNMs in general relativity. Since QNMs from binary BH mergers will be detected with high signal-to-noise ratio by GW detectors, it is also possible to detect the second perturbative order of QNMs, generated by nonlinear gravitational interaction near the BH. In the BH perturbation approach, we derive the master Zerilli equation for the metric perturbation to second order and explicitly regularize it at the horizon and spatial infinity. We numerically solve the second order Zerilli equation by implementing the modified Leaver's continued fraction method. The second order QNM frequencies are found to be twice the first order ones, and the GW amplitude is up to ~10% that of the first order for the binary BH mergers. Since the second order QNMs always exist, we can use their detections (i) to test the nonlinearity of general relativity, in particular the no-hair theorem, (ii) to remove fake events in the data analysis of QNM GWs and (iii) to measure the distance to the BH.