  In this paper, we developed a basis-independent perturbative method for
calculating the non-linear optical response of arbitrary non-interacting
tight-binding models. Our method is based on the non-equilibrium Keldysh
formalism and allows an efficient numerical implementation within the framework
of the Kernel Polynomial Method for systems which are not required to be
translation-invariant. Some proof-of-concept results of the second-order
optical conductivity are presented for the special case of gapped graphene with
vacancies and an on-site Anderson disordered potential.
