  Cu$_2$ZnSnS$_4$ is an earth-abundant photovoltaic absorber material predicted
to provide a sustainable solution for commercial solar applications. One of the
main limitations restricting its commercialization is the issue of cation
disorder. Raman spectroscopy has been a sought after technique to characterize
disorder in CZTS while a clear consensus between theoretical and experimental
results is yet to be achieved. In the present study, via the virtual crystal
approximation, we take into account the progressive nature of Cu-Zn disorder in
CZTS: we obtain the phonon frequencies at zone-center within the density
functional perturbation theory formalism, and further compute the Raman spectra
for the disordered phases, achieving a consensus between theory and experiment.
These calculations confirm the presence of complete disorder in Cu-Zn 2$a$,
2$c$ and 2$d$ Wyckoff sites. They also show that the Raman intensities of two
prominent $A$ phonon modes characterized by motion of S atoms, also known to be
experimentally significant, play a key role in understanding the nature of
disorder in CZTS.
