  Measurements of multiple harmonic modes in the gravitational wave signals
from binary black hole events could provide an accurate test of general
relativity, however they have never been observed before. The sub-dominant
modes, other than the main (l = 2, m = 2) mode, are weak in amplitude and thus
difficult to detect in a single event at the current sensitivity of the
gravitational wave detectors. To recover sub-dominant modes, we propose an
unmodeled method for summation of signals from a population of binary black
holes. The method coherently stacks all the signal modes, effectively
increasing their signal-to-noise ratio, so the amplified signal can be
extracted from the noisy data. To test the method, we consider simulated
numerical relativity waveforms including sub-dominant modes up to (5, 5). We
inject simulated signals from a population of binary black holes into data from
the first observing run of Advanced LIGO and utilize the coherent WaveBurst
algorithm for signal detection and reconstruction. Using only the waveforms
reconstructed by coherent WaveBurst, i.e., no a priori information about the
signal model, we determine the transformation to coherently synchronize the
merger and post-merger of one signal to another. We demonstrate the
synchronization of the injected signals and show the efficient stacking of the
(2, 2) mode and the sub-dominant modes.
