Tangent equations of motion for nonlinear response functions

Nonlinear response functions, formulated as multipoint correlation functions or Volterra kernels, encode the dynamical and spectroscopic properties of physical systems and underpin a wide range of nonlinear transport and optical phenomena. However, their evaluation rapidly becomes prohibitive at high orders because of combinatorial (often factorial) scaling or severe numerical errors. Here, we establish a systematic and efficient framework to compute nonlinear response functions directly from real-time dynamics, without explicitly constructing multipoint correlators or relying on numerically unstable finite-difference methods for order-resolved extraction. Our approach is based on the Gateaux derivative with respect to the external field in function space, which yields a closed hierarchy of tangent equations of motion (TEOM). Propagating the TEOM alongside the original dynamics isolates each perturbative order with high accuracy, providing a term-by-term decomposition of physical contributions. The computational cost scales exponentially with response order in the fully general setting and reduces to polynomial complexity when all perturbation directions are identical; both regimes avoid the factorial scaling of explicit multipoint-correlator evaluations. We demonstrate the power of TEOM by computing frequency-resolved fifth-order response functions for a solid-state electron model and by obtaining nonlinear response functions up to the 49th order with controlled accuracy in a classical Duffing oscillator. We further show that our time-evolution formulation allows optical conductivities to be evaluated directly while remaining numerically stable even near zero frequency. TEOM can be incorporated seamlessly into existing real-time evolution methods, yielding a general framework for computing nonlinear response functions in quantum and classical dynamical systems.