Computing nonlinear response to arbitrary order in the presence of electronic interactions via time-domain integration of polarizability

We describe a theoretical approach to study the nonlinear optical response of electronic systems based on a real-time solution of the electronic dynamics in the presence of time-dependent electromagnetic fields. Using accurately parameterized tight-binding models and electronic interactions, this allows expedite and accurate calculations of the nonlinear response to arbitrary order for realistic systems. We demonstrate its capabilities by computing the excitonic spectrum and high-harmonic generation (SHG, THG, FHG) in MoS$_2$ and BN monolayers. The linear optical conductivity is seen to reproduce exactly the Kubo-formula perturbative results based on the Bethe-Salpeter equation and, moreover, reproduces very well the experimental traces. Such benchmark against experiments in linear order allows us to confidently predict the frequency dependence of the higher-harmonic susceptibilities which, within our scheme, are obtained with no additional effort, in contrast to an order-by-order perturbative approach.