Lanczos-based computation of second- and third-order susceptibilities applied to extended Kitaev models

Recent advances in spectroscopic techniques such as two-dimensional coherent spectroscopy (2DCS) have highlighted the potential of higher-order response functions to reveal hidden features in quantum materials. Yet, the numerical evaluation of nonlinear susceptibilities beyond linear response remains challenging due to the steep scaling of computational cost. Building on our recently introduced method for computing second-order susceptibilities (arXiv:2502.01746), we here present its extension to third-order response. Our Lanczos-based framework enables the direct evaluation of second- and third-order susceptibilities in the multi-dimensional frequency domain. We compare several approximations and algorithmic strategies in terms of convergence behavior, numerical stability, and computational cost. As a case study, we apply the method to extended Kitaev models relevant to materials such as α-RuCl3, demonstrating how nonlinear response can distinguish between different types of excitations beyond what is accessible from linear response alone. Our results establish a practical route for studying higher-order dynamical response in arbitrary spin and Hubbard-like models.