  Sorting and assigning phonon branches (e.g., longitudinal acoustic) of phonon
modes is important for characterizing the phonon bands of a crystal and the
determination of phonon properties such as the Gr\"uneisan parameter and group
velocity. To do this, the phonon band indices (including the longitudinal and
transverse acoustic) have to be assigned correctly to all phonon modes across a
path or paths in the Brillouin zone. As our solution to this challenging
problem, we propose a computationally efficient and robust two-stage hybrid
method that combines two approaches with their own merits. The first is the
perturbative approach in which we connect the modes using degenerate
perturbation theory. In the second approach, we use numerical fitting based on
least-squares fits to circumvent local connectivity errors at or near exact
degenerate modes. The method can be easily generalized to other condensed
matter problems involving Hermitian matrix operators such as electronic bands
in tight-binding Hamiltonians or in a standard density-functional calculation,
and photonic bands in photonic crystals.
