We develop a general framework based on the functional derivative to extract nonlinear dynamical response functions from the temporal evolution of physical quantities, without explicitly computing multipoint correlation functions. We validate our approach by calculating the second- and third-order optical responses in the Rice-Mele model and further apply it to a many-body interacting system using a tensor network method. This framework is broadly applicable to any method that can compute real-time dynamics, offering a powerful and versatile tool for investigating nonlinear responses in dynamical systems.
We develop a general framework based on the functional derivative to extract nonlinear dynamical response functions from the temporal evolution of physical quantities, without explicitly computing multipoint correlation functions. We validate our approach by calculating the second- and third-order optical responses in the Rice-Mele model and further apply it to a many-body interacting system using a tensor network method. This framework is broadly applicable to any method that can compute real-time dynamics, offering a powerful and versatile tool for investigating nonlinear responses in dynamical systems.
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.
It is well established that causal linear response functions can be found by computing the much simpler imaginary time-ordered Matsubara functions and performing an analytic continuation. This principle is the basis for much of our understanding of linear response for interacting and disordered systems, via diagrammatic perturbation theory. Similar imaginary-time approaches have recently been introduced for computing nonlinear response functions as well, for example in [Annalen der Physik 536, 2300504 (2024); Physical Review X 11, 041006 (2021)], where the authors analytically continue the Matsubara functions to obtain the Keldysh response functions. In this work, we provide a proof of this connection to all orders in perturbation theory using an equation of motion based approach. We show by induction that causal nonlinear response functions at every order can be obtained from an analytic continuation of an appropriate time-ordered Matsubara function. We demonstrate this connection explicitly for second order response functions in the Lehmann representation. As a byproduct of our approach, we derive an explicit expression for the Lehmann representation of n-th order response functions by solving the equations of motion. We also use our result to find an analytic spectral density representation for both causal response functions and Matsubara functions. As an example, we apply our method to derive the non-linear A3 term in the SU(2) spin Hall response of an insulator with spin rotation symmetry. Finally, we show how our results lead to a family of generalized sum rules, focusing explicitly on the asymptotic expression for n-th harmonic generation rate. Our work opens the door to using imaginary time approaches to study nonlinear response functions in general condensed matter systems.
The quantum geometric tensor and quantum Fisher information have recently been shown to provide a unified geometric description of the linear response of many-body systems. However, a similar geometric description of higher-order perturbative phenomena including nonlinear response in generic quantum systems is lacking. In this work, we develop a general framework for the time-dependent quantum geometry of many-body systems by treating external perturbing fields as coordinates on the space of density matrices. We use the Bures distance between the initial and time-evolved density matrix to define geometric quantities through a perturbative expansion. To lowest order, we derive a time-dependent generalization of the Bures metric related to the spectral density of linear response functions, unifying previous results for the quantum metric in various limits and providing a geometric interpretation of Fermi’s golden rule. At next order in the expansion, we define a time-dependent Bures-Levi-Civita connection for general many-body systems. We show that the connection is the sum of one contribution that is related to a second-order nonlinear response function, and a second contribution that captures the higher geometric structure of first-order perturbation theory. We show that in the quasistatic, zero-temperature limit for noninteracting fermions, this Bures connection reduces to the known expression for band-theoretic Christoffel symbols. Our work provides a systematic framework to explore many-body quantum geometry beyond the quantum metric and highlights how higher-order correlation functions can probe this geometry.
Abstract Nonlinear spectroscopy is a cornerstone of quantum science, providing unique access to multi-point correlations, quantum coherence, and couplings that are invisible to linear methods. However, classical simulation of these phenomena is fundamentally limited by the exponential growth of the Hilbert space, and practical quantum algorithms for the nonlinear regime have remained largely unexplored. Here, we present a unified quantum algorithmic framework for computing n-th order nonlinear spectroscopies. By reformulating multi-time responses as a weighted sum of expectation values at finite pump amplitudes via a generalized parameter shift rule, our approach bypasses the costly evaluation of high-order commutators and time-dependent operator expansions. This reformulation enables efficient execution via real-time evolution on current quantum hardware, ensuring inherent noise resilience. We validate the framework on IBM's superconducting quantum processors, successfully obtain higher-order response functions of a 12-qubit XXZ spin-chain. Furthermore, the versatility of our method is demonstrated by resolving quasi-particle excitation spectra in spin-liquids and identifying interaction-induced cross-peaks in atomic systems. Our results establish a practical and scalable pathway for probing complex quantum dynamics on near-term quantum devices, extending the reach of quantum simulation into the nonlinear domain.
The problem of nonlinear response is considered by employing a general time-evolution equation, and a Green's function which is the transition or conditional probability density for an unperturbed system. Expansion of the Green's function in terms of orthonormal functions enables us to express the distribution function describing the nonlinear behavior by means of matrix products whose elements are composed of correlation functions in the absence of the perturbation. In other words, it is shown how the distribution function induced by a strong perturbation may be calculated by knowing the Green's function without the perturbation. As the special case of the linear response, we have obtained Kubo's relation. The Laplace-transform technique with respect to time is found quite useful in developing the present theory in which the transient effect is also taken into account. As an application of the theory, a new relation valid in the region of the second-order perturbation connecting the transient rise and decay with the stationary alternating perturbations has been obtained.
Nonlinear spectroscopy is a cornerstone of quantum science, providing unique access to multi-point correlations, quantum coherence, and couplings that are invisible to linear methods. However, classical simulation of these phenomena is fundamentally limited by the exponential growth of the Hilbert space, and practical quantum algorithms for the nonlinear regime have remained largely unexplored. Here, we present a unified quantum algorithmic framework for computing n-th order nonlinear spectroscopies. By reformulating multi-time responses as a weighted sum of expectation values at finite pump amplitudes via a generalized parameter shift rule, our approach bypasses the costly evaluation of high-order commutators and time-dependent operator expansions. This reformulation enables efficient execution via real-time evolution on current quantum hardware, ensuring inherent noise resilience. We validate the framework on IBM's superconducting quantum processors, successfully obtain higher-order response functions of a 12-qubit XXZ spin-chain. Furthermore, the versatility of our method is demonstrated by resolving quasi-particle excitation spectra in spin-liquids and identifying interaction-induced cross-peaks in atomic systems. Our results establish a practical and scalable pathway for probing complex quantum dynamics on near-term quantum devices, extending the reach of quantum simulation into the nonlinear domain.