Tangent equations of motion for nonlinear response functions

Abstract

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.

Figures (13)

• Figure 1: Schematic illustration of the Gateaux derivative of the state functional $\rho[f](t)$ (e.g., a density matrix in quantum dynamics, or a state variable in classical dynamics). The two horizontal axes represent a simplified two-dimensional slice of the function space $\mathcal{F}$ of external fields, and the vertical axis denotes the state space $\mathcal{M}$ of density matrices. For each time $t$, the set $\{\rho[f](t) \mid f \in \mathcal{F} \}$ forms a surface in $\mathcal{F} \times \mathcal{M}$ (orange at time $t$ and blue at time $t + \mathrm{d}t$). The blue and yellow curves are the trajectories $\rho[f](t)$ and $\rho[f{+}\epsilon g](t)$, respectively, where $g$ is the perturbation field. The cyan arrows on the gray tangent planes correspond to the Gateaux derivatives $\mathcal{D}_{g} \rho[f](t)$ and $\mathcal{D}_{g} \rho[f](t{+}\mathrm{d}t)$, which are tangent vectors induced by an infinitesimal variation of the external field in the direction of $g$. The time evolution of $\mathcal{D}_{g} \rho[f](t)$ is governed by the first-order tangent equations of motion (TEOM) [Eq. \ref{['eq:teom_1st']}]. More generally, TEOM extend the dynamical state (quantum or classical) from $\rho[f](t)$ to a finite-order collection of field-functional derivatives, which are propagated alongside $\rho[f](t)$ and used to reconstruct the nonlinear response functions.
• Figure 2: Temporal profiles of the electric current and its Gateaux derivatives with $g_{\alpha_1\!} = g_{\text{cos}_1\!}$ and $\omega_1 = 1.5$. (a) Electric current. The inset shows an enlarged view for $\vert t \vert < 1$, where the crosses indicate the data points calculated on the adaptive time grid. (b) First Gateaux derivative of the symmetrized electric current. (c) Second Gateaux derivative of the antisymmetrized electric current. The inset spanning panels (b) and (c) presents the absolute values of the currents at $t = 0$ as functions of the field amplitude, and the dashed lines indicate power-law scalings $\propto F_0$ (red), $\propto F_0^2$ (blue), and $\propto F_0^3$ (green).
• Figure 3: Convergence behavior of the second- and third-order response functions, for $\omega_1 = 1.5$. [(a) and (b)] Response functions obtained with Gaussian, exponential, and Hann windows, compared with perturbation-theory results. The inset in (a) is an enlarged view from $\omega = -1.1$ to $-0.95$. [(c) and (d)] Dependence on the field amplitude. [(e) and (f)] Dependence on the number of $k$ points. The inset in (e) shows an enlarged view from $\omega = -1.6$ to $-1.4$, where the red curve corresponds to $N = 800$.
• Figure 4: Dependence of spectral resolution on the temporal widths of the Gaussian perturbation field $\tau_{g}$ and the Gaussian window function $\tau_{w}$, for (a) $\tau_{g} = \tau_{w} = 100$, (b) $\tau_{g} = 100$ and $\tau_{w} = 10$, (c) $\tau_{g} = 10$ and $\tau_{w} = 100$, and (d) $\tau_{g} = \tau_{w} = 10$.
• Figure 5: Same data as Fig. \ref{['fig:resolution']}, but replotted as a function of $(\omega_1, \omega_2)$. (a) $\tau_{g} = 10$ and $\tau_{w} = 100$; (b) $\tau_{g} = 100$ and $\tau_{w} = 10$.