Thermodynamic constraints on kinetic perturbations of homogeneous driven diffusions

We analyze the static response to kinetic perturbations of nonequilibrium steady states that can be modeled as diffusions. We demonstrate that kinetic response is purely a nonequilibirum effect, measuring the degree to which the Fluctuation-Dissipation Theorem is violated out of equilibrium. For driven diffusions in a flat landscape, we further demonstrate that such response is constrained by the strength of the nonequilibrium driving via quantitative inequalities.

Introduction.-The question addressed by linear response theory is how a system reacts to a small perturbation [1,2].Traditionally the only perturbation considered was the application of a small force.The reason is that around equilibrium steady states, changes in kinetic parameters-like the mobility of a colloidal particle or the energy-barrier along a reaction pathway-are trivial in that they have no effect: the Boltzmann distribution only depends on the energies of the system's states.Out of equilibrium, this is not the case, and kinetic perturbations not only affect the steady-state distribution, but including their effects are necessary to completely capture nonequilibrium response.
We are learning now that explicitly analyzing kinetic perturbations can lead to quantifiable insight into nonequilibrium response.Indeed, the extent to which the Fluctuation-Dissipation Theorem (FDT) is broken out of equilibrium equals a kinetic response [3,4].In addition, for Markov jump processes [3,5,6], chemical reaction networks [7], and one-dimensional diffusions [4], we can put concrete limits (or bounds) on the kinetic response.It has not been shown, however, that the response remains bounded for driven diffusions away from equilibrium in dimensions higher than one.Here, we first highlight the fact that kinetic perturbations in arbitrary diffusion processes measure the degree to which the FDT is broken.Then for homogeneous diffusions in a flat landscape we derive bounds on the kinetic response in terms of the strength of nonequilibrium driving.These results demonstrate that at least for this class of diffusions, response is indeed bounded, and allow us to speculate that similar limits may hold in general for arbitrary diffusions.
Setup.-We have in mind D-dimensional systems whose configurations evolve with time according to a periodic diffusion process.For this class of systems, their configuration x(t) = (x 1 (t), . . .x D (t)) at time t takes values in the torus Ω = [0, L] × ⋯ × [0, L], and the Fokker-Planck equation describing the time-evolution of the probability p-1 arXiv:2404.09860v1[cond-mat.stat-mech]15 Apr 2024 density p(x, t) can be parameterized as [35] (2) Borrowing language from the modeling of a colloidal particle in a viscous fluid, we have identified in the Fokker-Planck operator L a positive-definite, mobility matrix μ(x) as well as having split the force into a conservative part due to a potential U (x) and a nonconservative part F(x) (∇ × F ≠ 0).We will assume, that the system relaxes to a unique steady-state distribution π(x) given as the solution of Lπ(x) = 0.In general, this steady-state solution is not known.However, when the nonconservative force is zero (F = 0), it is straightforward to show that the steady-state distribution is π eq (x) ∝ e −U (x) (with k B T = 1), which we will identify as an equilibrium distribution.
Linear response.-A central paradigm in statistical physics is to analyze a system by how steady-state averages of observables, ⟨Q⟩ = ∫ Q(x)π(x)dx, change in response to external perturbations.
Response has traditionally been modeled by assuming that the potential depends on an external parameter λ via U (x) → U (x) − λV (x), where we call the function V (x) the conjugate coordinate.In response, averages ⟨Q⟩ change by [2] For equilibrium systems, where F = 0, the response is completely captured by how the equilibrium distribution π eq (x; λ) ∝ e −(U (x)−λV (x)) is modified.The immediate consequence of this structure is the FDT, relating the static response of the equilibrium average of the observable ⟨Q⟩ eq = ∫ Q(x)π eq (x)dx to the fluctuations [2] R eq U = ∫ V (z) where the covariance is ⟪Q, V ⟫ eq = ⟨QV ⟩ eq − ⟨Q⟩ eq ⟨V ⟩ eq .Perturbations of the kinetics via the mobility μ(x), by contrast, have no effect as it does not enter the equilibrium distribution.Now, our previous analyses [3][4][5]7] have revealed that away from thermodynamic equilibrium it is in fact useful to consider how observables change in response to perturbations of the kinetics.This is implemented by allowing the mobility to depend on the external parameter instead, μ(x) → μ(x)(1 − λV (x)), Applying this perturbation in an experimental setting, where μ is the only system parameter varied, is likely challenging but may be possible.For example, the mobility of a colloidal particle depends on the viscosity of the surrounding fluid via the Stokes-Einstein relation.One could then imagine varying the mobility by mixing fluids of differing viscosities.Nevertheless, kinetic perturbations serve as an important intermediary in our analysis of energy perturbations (4) away from equilibrium.Here, they serve as a measure of the violation of the FDT [36,37] with ⟪Q, V ⟫ the nonequilibrium covariance.Now, R U and ⟪Q, V ⟫ are commonly measured in experiments.Their difference is a purely nonequilibrium effect captured by the kinetic response.
One might then surmise that a larger kinetic response requires stronger nonequilibrium driving.Indeed, our previous work has shown that there are such quantitative trade-offs, at least for one-dimensional diffusions on the circle.In this case, the only nontrivial form for the nonconservative force is a constant f .Then we have shown that when V (x) = δ (a,b) (x) is an indicator function on the range x ∈ (a, b) and choosing Q(x) ∈ [0, 1] for ease of presentation, the response to a kinetic perturbation can be bounded by the strength of the nonequilibrium driving [4] The product ⟨Q⟩(1 − ⟨Q⟩) represents the maximum variance in the observable possible with fixed mean via the Bahtia-Davis inequality ⟪Q 2 ⟫ ≤ ⟨Q⟩(1 − ⟨Q⟩) [38].Thus, (7) can be viewed as a trade-off between the response, fluctuations, and thermodynamic driving.Viewed another way, (7) is the maximum of the kinetic response over the system parameters, µ(x) and U (x), holding fixed the thermodynamic driving f and the observable's average ⟨Q⟩.
The derivation of ( 7) relied on having a closed-form analytic solution for the steady-state distribution in one dimension for arbitrary system parameters.Such a method does not translate to diffusions in higher dimensions where no such solution is known.This naturally raises the question of whether thermodynamic force is a constraint on the kinetic response for diffusions in higher dimensions, or if unbounded response is possible?We address this question in the next section.
Bounds on kinetic response.-Without a closedform solution to the Fokker-Planck equation (1) in dimensions higher than one, we make progress by focusing on a simpler class of models: homogenous driven diffusions in a flat landscape.Specifically, for the remainder of this article we specialize to the case where the mobility matrix is constant and diagonal, μ(x) = μ with diagonal elements {µ 1 , . . ., µ D }; there is no potential U (x) = 0; and the nonconservative driving is uniform F(x) = f = (f 1 , . . ., f D ).In this case, the dynamics of the probability distribution is determined by the Fokker-Planck operator (2) with constant coefficients.The translational-symmetry implies that the steady-state solution is uniform, π(x) = 1/|Ω|, which can also be verified by direct substitution.Despite the simplicity of the steady-state distribution, spatially-dependent perturbations of the mobility can still lead to complicated changes in the steady-state distribution.It is these responses that we aim to constrain.For kinetic perturbations of driven diffusions in a flat landscape, we have derived three thermodynamic limits.In contrast to our previous work [4], we cannot optimize over all system parameters, as they are fixed.Instead, we maximize the response over the observable Q(x) and conjugate coordinate V (x).To have a well-posed problem, though, we need to constrain Q(x) and V (x) in some way.In light of our previous results (7) and taking inspiration from the FDT, we fix their fluctuations via their variances, ⟪Q 2 ⟫ and ⟪V 2 ⟫.
Our first main result is the bound where F = max j |f j |L quantifies how far the system is out of equilibrium.As a sanity check, when F = 0, the system is at equilibrium, and the bound is zero, meaning there is no violation of FDT (cf.( 6)).We find that equality is reached when Q and V are given by sine waves in one direction (shifted by a phase), and uniform in the orthogonal directions.One example of such an optimal choice is illustrated in Figs.1(b) and (c) for two dimensions.Thus, the system appears most sensitive to slowly varying perturbations that align with the thermodynamic driving.
In our second main result, we specialize to a situation where both observable and conjugate coordinate are the same, Q = V .In this case, we have shown that the response satisfies the tighter inequality The improvement comes from the denominator being 1 . This is plausible considering that the observable and conjugate coordinate are more constrained.Despite this improvement, the most sensitive response is again reached for low wave number sine waves.
The third main result is for an even further restricted situation where Q and V are indicator functions.It should be stressed that the explicit form, presented below, is based solely on numerical observations for which we have no analytic proof.To be specific, we require Q(x) = V (x) = δ S (x), where δ S (x) takes the value one for x ∈ S ⊆ Ω and zero otherwise.In this case, we observed that This upper bound corresponds to the response in the case where half of the region Ω is perturbed (in a way speci- fied below).That there is a tighter inequality is reasonable, since we are focusing on a more restricted scenario.Numerical evidence of this improvement will be provided below.
In Fig. 1, we numerically verify inequality (9) for one, two, and three dimensional driven diffusions in a flat landscape.To numerically determine the response, we exploit the translational invariance of the problem to calculate the response of the steady-state distribution in Fourier space, δπ m with m ∈ Z D , and then estimate the response from the sum over Fourier coefficients R µ = ∑ m δπ m Q −m (see the discussion leading to (18) below).As this is done numerically, the number of Fourier modes that we can keep in this analysis is finite.We implement this restriction by keeping only modes with m below a cut-off max j |m j | ≤ N F , which in this study we take to be one, two or three.Then, for each combination of N F and D, we generate 1000 random samples by varying the system parameters.We observe that in Fig. 1 the bound is saturable for every value of driving force.Moreover, we observe that potential optima are confined to the lowest space of modes with N F = 1, where both observable Q and the conjugate coordinate V are sine waves.This is why we tend to only observe saturation of the inequality when we restrict our sampling to N F = 1, and why our samples fall away from the optimal curve for higher values of N F .Our second inequality ( 10) is numerically verified in Fig. 2, again in one, two, and three dimensions.Sampling is the same as described for Fig. 1, except for fixing V = Q.All points fall below the predicted limit (10) given by the black line.
Limits to the response when the observable and conjugate coordinate are indicator functions are studied numerically in Fig. 3. To numerically specify the region S of the indicator function we divide the entire space Ω into N D B equally sized cubic regions we call blocks, with N B the number blocks along each linear dimension.We can then form a region S by combining together a collection of these blocks, setting the observable to one on the selected blocks.For each combination of D = {1, 2, 3} and N B = {2, 6, 10}, we generate 1000 random samples, choosing f and μ as before.The perturbation region S is built up by randomly adding each block to S with probability 1/2.The response-ratio |R µ |/⟪δ 2 S ⟫ is then determined numerically with N F = 20 in order to have a Fourier scale finer than the smallest block size.The data is plotted in Fig. 3(a) as a function of force F. When D = 1 and N B = 6, there are visible line structures in the plot.This is because there is a small number of different ways of choosing the blocks, which fall on a single line.As N B increases, these structures remain but are hard to dinstiguish visually.For D ≥ 2, each F corresponds to infinitely many f so the lines disappear.
The red line is the analytic bound predicted in (10).However, it is nowhere saturated.Instead all sampled responses appear to share a tighter, saturable bound, depicted by the black line.We visually inspected the samples that appear optimal, like the starred data point in Fig. 3(a), and noticed that they consistently all had a perturbation region that occupied half of the full space Ω: for example, the perturbation region for the starred data point is depicted in Fig. 3(b).Moreover, we observed that the maxima reached for all the responses is the same for every dimension tested, including 1D.This suggested to us that we could predict the maximum response possible just from the analytic solution for the 1D ring with perturbation region S = [0, 1/2], which is in fact the black line.
Discussion. -We analyzed three different kinetic perturbation schemes of the mobility for homogeneous nonequilibrium diffusions.For this class of models, we found that the response is not unlimited: it must smoothly approach zero as we near equilibrium and has a maximum arbitrarily far from equilibrium constrained by the size of the fluctuations of the observable and the perturbation's conjugate coordinate.
Natural extensions would be to have a nonuniform or nondiagonal mobility matrix, or to allow for a anisotropic perturbation.Going further and considering arbitrary potentials with a uniform nonconservative force, −∇U (x)+f , remain out of reach with the current method, as the Fokker-Planck operator L is not diagonal in the Fourier basis.Moreover, any bounds for these more general setups would likely take a different form.Indeed, we have verified that our bounds can be violated in 1D in the presence of a nonuniform potential.Thus new tools are needed, though the current results suggest progress can be made.One potential avenue is the recently developed linear-algebraic technique proposed in [6] to tackle single edge perturbations for finite-state Markov chains.An extension to diffusion processes would be intriguing, but is far from trivial.Even when allowing for multi-edge perturbations, taking the diffusive limit presents challenges [4].
Derivation of general variance bound.-In this section, we derive the bound on response in terms of variances (9).Without loss of generality, we set the length in each dimension to L = 1,to simplify the presentation.
To proceed, we explicitly allow the Fokker-Planck operator to depend on the small external parameter via with the accompanying steady-state solution π λ satisfying L λ π λ = 0. Since the unperturbed steady-state distribution π is uniform, π λ = π + λδπ is nearly uniform with δπ = ∂ λ π λ | λ=0 .With this notation, the response at λ = 0 can be expressed as This expression transfers the problem to determining how the steady-state distribution responds, δπ.A convenient way to obtain this steady-state response is to differentiate the Fokker-Planck equation at λ = 0, after using π = 1 is uniform.This is a partial differential equation for δπ [5,6,39], which we now proceed to solve.Because the operator L is linear in the Fourier basis, a compact solution to (16) can be found by Fourier transform.To this end, let us denote the Fourier basis as e m (r) = e i2πm⋅r for m ∈ Z D .Because they form an orthonormal basis, ⟨e m , e n ⟩ = ∫ e m (z)e * n (z)dz = δ mn , any periodic function G(r) can be expanded as Furthermore, for real G, the coefficients satisfy G * m = G −m .This construction is convenient, because once we have determined the Fourier coefficients of the steady-state response δπ m , the Fourier expansions in (17) can be substituted into (15) to obtain the response formula Now to solve (16), we expand both sides in the Fourier basis We first note that in this basis the Fokker-Planck operator is diagonal, Le m (r) = l m e m (r), (20) with eigenvalues l m = −4π 2 (m ⋅ μ ⋅ m) − 2πim ⋅ μ ⋅ f .Next, we evaluate the right hand side using (16), which reads Combining these observations with the orthogonality of the Fourier basis leads to a series of uncoupled linear equations for δπ m , which can be solved for m ≠ 0 with δπ 0 = 0 due to probability conservation.Now substituting this solution into (18) leads to our desired starting point for deriving bounds, with We now bound the maximum of ( 24) over all V m and Q m , with their variances constrained, To proceed, we introduce re-weighted Fourier coefficients Ṽm = √ λ m V m and Qm = √ λ −m Q m , allowing us to apply the Cauchy-Scwharz inequality to (24): ≤ max All that is left is to bound the magnitude of λ m .Let the real and imaginary parts of λ m be ρ m and σ m : λ m = ρ m + iσ m .Then the squared magnitude of λ m , which here equals its real part ρ m , can be written as To bound |λ m | 2 , we study the quantity in the denominator

Fig. 1 :
Fig. 1: Illustration of the general bound (9): (a) Random samples of the response |Rµ| plotted as a function of the force F. Samples are generated by choosing f uniformly on [0, 50] D , μ uniformly on [0, 1] D , and the Fourier coefficients of V and Q for all allowed modes uniformly on [0, 1], while maintaining proper parity so that the corresponding V and Q are real.These samples are then normalized to keep ⟪Q 2 ⟫ = ⟪V 2 ⟫ = 1.The colors label the highest Fourier mode N F possible in any direction, while the shapes label the dimension D. Each combination of color and shape contains 1000 data points.(b) & (c): Example of the optimal V and Q that saturate the bound (9) for D = 2, with f = (fx, fy), |fx| > |fy|.

Fig. 2 :
Fig. 2: Illustration of the bound (10): Random samples of the response |Rµ| plotted as a function of the force F. Samples are generated in the same manner as Fig. 1 except further restricting V = Q and normalizing ⟪Q 2 ⟫ = 1.The colors label the highest Fourier mode N F possible in any direction, while the shapes label the dimension D. Each combination of color and shape contains 1000 data points.

Fig. 3 :
Fig. 3: Illustration of the bound in (11): (a) Random sampling of μ, and perturbation region S.Each combination of color and shape contains 1000 samples, excluding trivial ones for which |Rµ| = ⟪δ 2 S ⟫ = 0. We divide Ω into N D B blocks and randomly select each block with probability 0.5 as the perturbation region.The smaller N D B , the more probable it is to form an optimal perturbation region to saturate the numerical bound.(b) The perturbation region (blue) that saturates the bound, corresponding to the data point labeled by the star in (a).