  Massive black holes (MBHs) are crucial in shaping their host galaxies. How
the MBH co-evolves with its host galaxy is a pressing problem in astrophysics
and cosmology. The valuable information carried by the binary MBH is encoded in
the gravitational waves (GWs), which will be detectable by the space-borne GW
detector LISA. In the GW data analysis, usually, only the dominant $(2,2)$ mode
of the GW signal is considered in the parameter estimation for LISA. However,
including the higher harmonics in parameter estimation can break the degeneracy
between the parameters, especially for the inclination angle and luminosity
distance. This may enable the identification of GW signals without
electromagnetic counterparts, known as ''dark sirens''. Thus, incorporating
higher harmonics will be beneficial to resolve the Hubble tension and constrain
the cosmological model. In this paper, we investigate the role of higher
harmonics in the parameter estimation for GWs emitted by binary MBHs. We
demonstrate that including $(3,3)$ mode can lead to a $10^3$-times improvement
in angular resolution and a $10^4$-times improvement in luminosity distance.
Meanwhile, our results indicate that considering higher harmonics increases the
probability of identifying over 70% host galaxies from $10^{-2}\,\rm{Gpc}^3$
cosmological volume threshold (corresponding $10^5$ host galaxies), while the
probability less than 8% for only the $(2,2)$ mode. Thus, our results
underscore the importance of including higher modes in the GW signal from
binary MBHs, for LISA at least $(3,3)$ mode.
