Higher-order response theory in stochastic thermodynamics and optimal control

TL;DR

This work clarifies how linear-response-based methods for designing minimal-work protocols fit within a general Volterra expansion of the nonequilibrium density, and it systematically constructs higher-order (quadratic) response terms to test their usefulness. By applying the framework to an overdamped harmonic oscillator, the authors show that first-order (linear) results recover known slow-driving and weak-driving approaches, while including second-order terms yields only marginal improvements at a high computational cost and can introduce unphysical negative excess work. The study highlights fundamental limitations of near-equilibrium expansions for far-from-equilibrium driving and discusses connections to recent optimal-transport/thermodynamic-geometry perspectives as potential alternatives. The findings inform the design of nanoscale protocols and biomolecular schemes, suggesting caution with higher-order corrections and pointing toward transport-based methods for robust optimization across regimes.

Abstract

Linear response theory has found many applications in statistical physics. One of these is to compute minimal-work protocols that drive nonequilibrium systems between different thermodynamic states, which are useful for designing engineered nanoscale systems and understanding biomolecular machines. We compare and explore the relationships between linear-response-based approximations used to study optimal protocols in different driving regimes by showing that they arise as controlled truncations of a general causal response (Volterra) expansion. We then construct higher-order response terms and discuss the drawbacks and utility of their inclusion. We illustrate our results for an overdamped particle in a harmonic trap, ultimately showing that the inclusion of higher-order response in calculating optimal protocols provides marginal improvement in effectiveness despite incurring a significant computational expense, while introducing the possibility of predicting arbitrarily low and unphysical negative excess work.

Figures (4)

• Figure 1: Approximation regimes in optimal stochastic thermodynamics. Following dibyenduperscombonanca2018, we can parameterize the control protocol as $\lambda(t) = \delta\lambda(t) + \lambda_i$, and $\tau_R$ represents an appropriate relaxation time for the system. Region 1 corresponds to weak driving, region 2 corresponds to slow driving, and region 3 corresponds to protocols that drive the system further from equilibrium.
• Figure 2: Derivation of the linear response function.a. The control parameter $\lambda^a$ moves from $\lambda^a(t) - \epsilon^a$ to $\lambda^a(t)$ at time $t'$. If we assume $\lambda^a(\tau) = \alpha(\tau)$, where $\alpha(\tau)$ is the stiffness of a simple harmonic oscillator, the inverse variance of the nonequilibrium distribution, $1/\sigma^2(\tau < t')$, is equal to $\alpha(\tau < t')$ at equilibrium. Once the control parameter changes, $\delta \rho \neq 0$ and $1/\sigma^2(\tau) \neq \alpha(\tau)$. After this, over time $\delta \rho$ relaxes back to zero and $1/\sigma^2(\tau)$ relaxes to the new value, $\alpha(t)$. This picture is exact for the simple harmonic oscillator, but in general, this serves as a cartoon as the nonequilibrium density will not correspond to an equilibrium density for the same system for any control parameter value, requiring more than a single parameter to characterize their difference. b. The evolution of the nonequilibrium density $\rho(\tau)$ (dotted red line) is shown along with the equilibrium density $\rho^\mathrm{eq}(\tau)$ (solid blue line) corresponding to the current control parameter value (as specified in the upper right corner of the upper panels), with their difference indicated by the purple shaded area. In the bottom panels, assuming $\lambda^a(\tau) = \alpha(\tau)$, the stiffness of a simple harmonic oscillator, the equilibrium distribution defined by the potential abruptly changes from $\rho^\mathrm{eq}_{\lambda^a(t)-\epsilon^a}$ to $\rho^\mathrm{eq}_{\lambda^a(t)}$. Initially $\delta \rho = 0$. After $\tau = t'$, $\delta \rho$ decays back to zero as $\rho(\tau) \rightarrow \rho^\mathrm{eq}(t)$.
• Figure 3: Derivation of the quadratic response function.a. Response of the system (dotted red line) to two perturbations of a single control parameter $\lambda^a$ (solid blue line). For example, if the control parameter were the stiffness of a harmonic oscillator, the inverse variance of the nonequilibrium distribution would react to changes in this parameter. Again, this picture serves as a cartoon for systems with densities not fully characterized by a single parameter. b. We also consider perturbations to each of two different control parameters, $\lambda^a$ and $\lambda^b$.
• Figure 4: Comparison of different approximations. Solid black = true optimal protocol from seifert2007, blue dotted = Sivak-Crooks geodesic optimal protocol [Eq. \ref{['eq:geo_work']}], red dot-dash = Bonança-Deffner weak-driving optimal protocol [Eq. \ref{['eq:weak_lin_stiff_work']}], purple dashed = unconstrained linear response optimal protocol [Eq. \ref{['eq:linear_stiff_work']}], pink dotted = unconstrained quadratic response optimal protocol [Eq. \ref{['eq:quadratic_stiff_work']}]. All protocols are computed for a time discretization of 100 steps and a factor of 5 increase in the trap stiffness $\alpha_f = 5 \alpha_i$. a. Optimal protocols for $t_f=1$ according to the different approximations. b. For different protocol durations, the predicted excess work according to the approximations from Eqs. \ref{['eq:linear_stiff_work']}, \ref{['eq:weak_lin_stiff_work']}, \ref{['eq:thermo_geo_stiff_work']}, and \ref{['eq:quadratic_stiff_work']}. c. Actual excess work for optimal protocols suggested by each approximation at a given protocol durations, evaluated by numerically solving the Fokker-Planck equation. d. Deviation of the true excess work for approximated optimal protocols from exact optimal protocol, as evaluated by numerically solving the Fokker-Planck equation. Note the good performance of weak driving protocols for short times even for strong protocols (despite poor match between the approximated work as shown in panel 4.b.), the good performance of the slow-driving approximations for long times, and the modest improvement of quadratic response over linear response which peaks at a $\sim$$12\%$ decrease in the deviation from the analytically optimal excess work.