Thermodynamic bounds on spectral perturbations

In discrete Markovian systems, many important properties (such as timescales of autocorrelation, relaxation, and oscillation) depend on the spectrum of the rate matrix. We demonstrate the existence of a universal trade-off between thermodynamic and spectral properties. We first show that the entropy production rate, the fundamental thermodynamic cost of a nonequilibrium steady state, bounds the difference between the eigenvalues of a nonequilibrium rate matrix and a reference equilibrium rate matrix. Using this result, we derive thermodynamic bounds on the spectral gap and on the magnitude of the imaginary eigenvalues. We illustrate our approach using a simple model of biomolecular sensing.


I. INTRODUCTION
One of the main goals of nonequilibrium thermodynamics is to understand trade-offs between thermodynamic costs and functionality in molecular systems [1,2].The entropy production rate (EPR) is one of the most important costs, since it quantifies both thermodynamic dissipation [3] and statistical irreversibility [4].It has been found to constrain functional properties such as accuracy [5,6], sensitivity of nonequilibrium response [7], and precision of fluctuating observables [8].
Here we consider discrete-state Markovian systems, which are commonly employed in stochastic thermodynamics, biophysics, chemistry, and other fields [9][10][11].In such systems, many important properties are determined by the spectrum (set of eigenvalues) of the rate matrix.For concreteness, one may imagine a molecular system that transitions between a finite number of configurations, such as the push-pull model shown in Fig. 1(a) and studied in our examples below.In general, the system is described by a probability distribution p(t) = (p 1 (t), . . ., p n (t)) with dynamics ∂ t p(t) = W p(t), where W is the rate matrix.The system's evolution over time τ is given by p(τ ) = e τ W p(0). Assuming W is irreducible and diagonalizable [12], the time-evolution operator can be expressed as where π is the steady-state distribution while (λ α , u (α) , v (α) ) refer to the eigenvalues and right/left eigenvectors of W .In this way, relaxation toward steady state is decomposed into contributions from different eigenmodes, with mode α decaying with rate −Re λ α ≥ 0 and oscillatory frequency Im λ α /2π [14][15][16].The spectrum of the rate matrix also determines other aspects of relaxation.For instance, degenerate eigenvalues can lead to power-law decay times [17,18] as well as dynamical phase transitions in relaxation trajectories [19].* artemyk@gmail.com The spectrum also plays a role in autocorrelation functions.The steady-state autocorrelation of observable a over time lag τ is ⟨a(τ )a(0)⟩ − ⟨a⟩⟨a⟩ = i,j ([e τ W ] ji − π j )π i a j a i . (2) Combining Eqs.(1) and (2) implies that the real parts of the eigenvalues control the decay rates of autocorrelations, while the imaginary parts control their oscillations [20,21].
Given the importance of the spectrum in various functional properties, we ask whether there is any relationship between a system's spectral and thermodynamic properties.In this paper, we prove the existence of this kind of relationship and explore some of its physical consequences.As one consequence, we demonstrate that steady-state EPR bounds the size of the spec- Figure 1.(a) Our results imply thermodynamic trade-offs with relaxation speed and oscillations in discrete-state systems, such as the 4-state "push-pull" model of enzyme concentration sensing [7].(b) Our thermodynamic bounds on the spectral perturbation (7) are illustrated on randomly sampled rate matrices.The far-from-equilibrium bound κ × σ ≥ 2∥∆λ∥ tanh -1 ∥∆λ∥ is shown in green and the near-equilibrium bound κ × σ ≥ 2∥∆λ∥ 2 is shown in purple.
tral gap, implying a fundamental thermodynamic cost for fast relaxation and for rapid decay of autocorrelation.We also show that steady-state EPR bounds the size of the imaginary part of the spectrum, implying a fundamental thermodynamic cost for fast oscillations.Both consequences follow from a general relationship between thermodynamic cost and spectral perturbations, illustrated in Fig. 1(b).This general relationship shows that the steady-state EPR under the nonequilibrium rate matrix W bounds the difference between the spectrum of W and a "reference" equilibrium rate matrix W , which has the same steady state and dynamical activity as W but obeys detailed balance.

II. PRELIMINARIES
Given a rate matrix W with steady-state π, the steady-state EPR is defined as [3,9] For simplicity, we assume that W is irreducible and "weakly reversible" (W ji > 0 whenever W ij > 0) ensuring that π is unique and σ is finite.The steady-state EPR vanishes when forward and backward probability fluxes balance for all i, j, π i W ji = π j W ij .This condition is termed detailed balance (DB), and we identify it with equilibrium.
The quantity σ can be interpreted as the rate of increase of thermodynamic entropy when the master equation obeys "local detailed balance" (LDB).Formally, LDB holds if for each transition i → j, the log-ratio ln(W ji /W ij ) can be identified with the corresponding increase of thermodynamic entropy in the environment.(See Ref. [32] for more details and a microscopic derivation of LDB.) LDB can also be considered for coarse-grained systems [33,34] and systems in contact with multiple thermodynamic reservoirs [35], in which case σ provides a lower bound on the rate of increase of thermodynamic entropy.If LDB does not hold, σ can still serve as a useful information-theoretic measure of time irreversibility [4].
In our first result below, we relate σ to the difference between the spectrum of the rate matrix W and that of a "reference" rate matrix Wij := (W ij +W ji π i /π j )/2.The rate matrix W obeys DB while having many of the same dynamical properties as W : the same steady state π, same escape rates Wii = W ii , and same dynamical activity across each transition.The dynamical activity across transition i → j is the rate of back-and-forth jumps [36], The rate matrix W is the equilibrium analogue of W , and it is equal to W if and only if W obeys DB.It is sometimes called the "additive reversibilization" of W in the literature [37,38].
For any rate matrix such as W , we use λ W to indicate the vector of eigenvalues of W sorted in descending order by real part, Re λ W 1 ≥ • • • ≥ Re λ W n .We will study the spectral perturbation between W and W , i.e., the difference between their eigenvalues.There are some subtleties involved in quantifying this difference, since in general there is no canonical one-to-one mapping between the eigenmodes of these two matrices.Following previous work [39,40], we quantify it in terms of the distance between the sorted eigenvalue vectors: ∥∆λ∥ := ∥λ W − λ W ∥. Importantly, ∥∆λ∥ is the minimal distance across all possible one-to-one mappings between the eigenvalues, ∥∆λ∥ = min Π ∥λ W − Πλ W ∥, where Π runs over all permutation matrices [41].In physical terms, ∥∆λ∥ compares the slowest modes of W with the slowest modes of W , and the fastest modes with the fastest modes.
Finally, we introduce two rate constants that will play a central role in our analysis.The first rate constant, which we term the mixing rate, is defined as The mixing rate quantifies the maximum dynamical activity across any edge, after an appropriate normalization by the steady-state probabilities.As shown in Appendix A, it can be understood as the maximum speed with which probability can flow between different regions of state space.The second rate constant is It reflects the overall magnitude of the normalized dynamical activity, where the contribution of each transition is weighted by the steady-state distribution.These two constants depend only on the dynamical activity and steady-state distribution and can be equivalently defined either in terms of W or W .

III. THERMODYNAMIC BOUND ON SPECTRAL PERTURBATIONS
In our first main result, we demonstrate that the steadystate EPR bounds the magnitude of the spectral perturbation The first bound is derived using a classic spectral perturbation theorem from linear algebra [39], see derivation in Appendix C 1.The second bound follows from x tanh -1 x ≥ x 2 , which is tight for small x.The two bounds become equivalent near equilibrium, when σ ≈ 0, W ≈ W , and ∥∆λ∥ ≈ 0.
Eq. ( 7) implies that there is a universal thermodynamic cost for having the eigenvalues of W be different from those of its equilibrium analogue W . Importantly, all terms that appear in our bounds (σ, κ, η W , and ∥∆λ∥) are functions of the system's rate matrix and they all have units of "per time".Moreover, all of these terms increase in the same linear manner when the rate matrix is scaled as W → αW .Therefore, the bounds ( 7) are invariant to the overall scaling of the rate matrix, depending instead only on its structure.
Observe that steady-state EPR (3) depends both on the system's kinetic properties, via the fluxes π i W ji − π j W ij , and its thermodynamic properties, via the dimensionless ratios ln(π i W ji /π j W ij ) (which are called "thermodynamic forces" in the literature [35]).The rate constants η W and κ serve as normalization terms for the kinetics, allowing EPR to be compared to the spectral perturbation in a dimensionless manner.Near equilibrium, only the mixing rate κ plays a role, while the rate constant η W becomes relevant in the far-from-equilibrium regime.The dimensionless ratio η W /κ in (7) measures the disorder of the normalized dynamical activity, and it reaches its maximum when the values of A ji /2π i π j are equal for all allowed transitions.
We illustrate our result in Fig. 1(b).We generate 10 5 random 5-by-5 rate matrices, scaled to have η W = 1.We plot κ × σ versus two lower bounds that follow from ( 7): 2∥∆λ∥ tanh -1 ∥∆λ∥ and 2∥∆λ∥ 2 .Observe that the first bound can be arbitrarily tight far from equilibrium.In fact, one can verify that it is saturated by unicyclic rate matrices with uniform rates.
We can also invert Eq. ( 7) to derive an upper bound on the spectral perturbation as a function of the EPR.Define the function Φ y : R ≥0 → R ≥0 as so that (7) implies κσ ≥ Φ η W (∥∆λ∥).We can then bound the spectral perturbation using EPR as where Φ −1 η W is the inverse function (which can be computed numerically from Φ η W ).

A. Intermediate bounds
The constant η W depends on the overall pattern of dynamical activity across all transitions, which can be difficult to access in practice.Below, it will be useful to derive more general bounds that do not depend on this constant.We define another rate matrix G, (10) This rate matrix obeys DB and has the same steady-state distribution and the same graph topology as W (the same pattern of allowed transitions with non-zero rates).At the same time, it does not depend on the dynamical activity under W .The constant η W can then be bound as where η G is defined as in ( 6) but using G instead of W (see Appendix B).Plugging these inequalities into (9) and using that Φ −1 y (x) is increasing in y gives a hierarchy of inequalities which fall in-between the two expressions in (9): Within this hierarchy, tighter bounds reflect finer-grained information about the rate matrix W .The weakest bounds depend only on the mixing rate κ, the intermediate bound also depends on the steady-state distribution and the graph topology (via η G ), while the tightest bound depends on the overall pattern of dynamical activity (via η W ). In principle, there are systems in which all the bounds expressed in terms of Φ −1 y are tight far-from-equilibrium.

B. Real and imaginary eigenvalues
The spectral perturbation can be decomposed into contributions from the real and imaginary parts of each eigenvalue, where ∆Re In this expression, we used λ W 1 = λ W 1 = 0 and that the eigenvalues of W are real valued (Appendix C).
In the following, we use our general results to study separate trade-offs between EPR versus decay timescales and EPR versus oscillations.Interestingly, (12) also implies a three-way trade-off between EPR, decay timescales, and oscillations, exploration of which we leave for future work.

IV. SPECTRAL GAP
We now derive thermodynamic bounds on the spectral gap |Re λ W 2 |, which controls the slowest decay timescale of relaxation and autocorrelation.
It is often useful to bound the spectral gap of W itself, rather than the increase of the spectral gap of W relative to W .To do so, we combine Φ −1 η W (κσ) ≤ Φ −1 κη G (κσ), as follows from (11), with the following inequalities on the spectral gap of W (Appendix C 3): where |λ G 2 | is the spectral gap of the rate matrix G from (10).Combining with (13) and rearranging gives These bounds depend only on the mixing rate κ, steady state π, and graph topology of the rate matrix.They include one term that is a thermodynamic cost, which depends on the EPR and vanishes in equilibrium, and a second "baseline" term, which does not vanish in equilibrium.Weaker but simpler bounds can also be derived using ( 11) and (15).For instance, we may derive bounds which depend only on the mixing rate as Eqs. ( 16) and ( 17) together form our second main result.

A. Example
We illustrate our bounds on the spectral gap by studying the thermodynamic cost of response speed in biomolecular sensing [7,[51][52][53][54][55].We consider a 4-state system with a cyclic topology, Fig. 1(a), inspired by the "push-pull" enzymatic sensor in the regime of low substrate concentrations [7].Here, a substrate molecule can be either unmodified S or modified S * .Modification is catalyzed by enzyme E 1 via bound state E 1 S, and demodification is catalyzed by enzyme E 2 via bound state E 2 S.
The steady-state distribution depends on the concentrations of E 1 and E 2 and serves as the sensor readout.Following a change of enzyme concentrations, the sensor's response speed is determined by the rate of relaxation toward the new steady state, which can be quantified using the spectral gap |Re λ W 2 |.The same idea applies not only to the push-pull system but also other molecular sensors where the steady state reflects environmental parameters.
We use numerical optimization to find 4-by-4 cyclic rate matrices that have the largest spectral gap for a given steady-state EPR, given a fixed steady-state distribution π and mixing rate κ.For concreteness, we fix π = (0.4,0.1, 0.4, 0.1) and κ = 1.Fig. 2(left) compares our bounds on the spectral gap (16) with the largest gap found by numerical optimization at different levels of EPR.Fig. 2(right) also shows that increased steady-state EPR is actually associated with faster relaxation.Specifically, using the rate matrices identified using numerical optimization, we plot D(p(t)∥π) = i p i (t) ln[p i (t)/π i ].This is the Kullback-Leibler (KL) divergence between the time-dependent probability distribution p(t) (starting from p(0) = (1, 0, 0, 0)) and the steady state π.As seen in Fig. 2, our bounds are always satisfied but only saturate at σ = 0. Interestingly, the spectral gap plateaus above σ ≈ 0.5, when it reaches the largest value achievable for a given topology and steady state.Beyond this point, additional EPR leads to changes in other eigenvalues and increases the imaginary part of λ 2 , as evidenced by oscillatory relaxation dynamics that emerge at large EPR in Fig. 2(right).
In future work, it may be interesting to relate our thermodynamic bound on response speed to existing thermodynamic bounds on sensitivity (ability of small parameter changes to elicit large changes in the steady state) [7,[51][52][53][54].

V. IMAGINARY EIGENVALUES
In our final set of main results, we show that the steady-state EPR bounds the magnitude of the entire imaginary part of the spectrum of W , where . This follows immediately from ( 9) and (12).Intermediate bounds which do not depend on the constant η W can be derived using (11).
We can also bound the imaginary part of any particular eigenmode (e.g., the slowest or the fastest).Consider the size of the imaginary part of any eigenvalue, |Im λ W α | for α ∈ {2, . . ., n}.Since nonreal eigenvalues come in conjugate pairs, there must be another eigenmode , which is tight for n ≤ 4 and becomes increasingly weak for systems with many states.Plugging into (18) gives

A. Example
As an illustration, we study the trade-off between steadystate EPR and oscillatory behavior, as quantified by the imag- inary part of the slowest eigenmode, |Im λ W 2 |.We again consider the 4-state cyclic system studied above, Fig. 1(a).Although oscillatory behavior is not typically desired in biomolecular sensors, it is necessary for other types of systems modeled using similar cyclic rate matrices, such as biochemical clocks [14].
We use numerical optimization to find 4-by-4 cyclic rate matrices that maximize |Im λ W 2 | for a given EPR, given a fixed steady-state distribution and mixing rate κ (chosen as above).In Fig. 3(left), we compare our bounds (19) with the largest imaginary part found by numerical optimization at different levels of steady-state EPR.In Fig. 3(right), we select some rate matrices identified using numerical optimization and plot Observe that relaxation dynamics exhibit oscillatory "pulsing", whose frequency increases with steady-state EPR.

VI. DISCUSSION
In this paper, we proved a relationship between the thermodynamic and spectral properties of a rate matrix.Our results on imaginary eigenvalues contribute to the extensive literature on the relationship between dissipation and oscillations [14,15,21,[27][28][29][30][31].
More surprising, perhaps, is our thermodynamic bound on the real eigenvalues, including the increase of the spectral gap.This is because EPR quantifies the violation of time-reversal symmetry, which at first glance appears to be unrelated to the real part of the eigenvalues.These results contribute to the growing interest in time-symmetric phenomena [36,56,57] and relaxation timescales [18,19,50,58,59] as signatures of nonequilibrium.
In future work, it may be interesting to consider nonlinear systems, e.g., by studying the relationship between EPR and the Jacobian spectrum in deterministic chemical reaction networks.It is also interesting to extend our approach to continuous-state systems (Fokker-Planck generators) and open quantum systems (Lindblad generators).Interestingly, a recent preprint derived a different type of thermodynamic bound on relaxation in Fokker-Planck dynamics [59]; comparison with our approach is left for future work.
To show the equivalence of ( 5) and (A1), we first bound the numerator in (A1) as In the last line, we used the definition of dynamical activity (4).We also restricted the maximization to i ̸ = j by using the fact that A ii ≤ 0 (since W ii ≤ 0).The denominator of the objective in (A1) is given by Combining gives the following upper bound on the right side of (A1), where the last expression is the definition of κ from (5).It can be verified that this upper bound is achieved by choosing X = {i * } and Y = {j * } in (A1), where Therefore, ( 5) and (A1) are equivalent.We first derive the following bound: In the second line, we used the assumption W is weakly reversible (W ij > 0 whenever W ji > 0), and in the last line we used the definition of κ from (5).Next, observe that i̸ =j:Wji>0 where we used the definition of the rate matrix G (10).Thus, we recover the first inequality in (11).
The second and third inequalities follow from i̸ =j:Wji>0

APPENDIX C: SPECTRAL INEQUALITIES
We begin by reviewing some of the notation and definitions introduced in the main text, along with a few useful results from linear algebra.
For any n-by-n matrix M , we use the notation λ M to indicate the vector of eigenvalues of matrix M sorted in descending order by real part, Re λ M 1 ≥ • • • ≥ Re λ M n .Given the rate matrix W , we define the reference equilibrium matrix Assuming W has a unique steady-state distribution π, the first eigenvalue is λ W 1 = 0, and Re λ W α < 0 for all α ∈ {2, . . ., n}.The right eigenvector corresponding to λ W 1 is π = (π 1 , . . ., π n ) T so that W π = 0.The corresponding left eigenvector is 1 = (1, . . ., 1) T so that 1 T W = 0. Any other right eigenvectors are orthogonal to 1, and any other left eigenvectors are orthogonal to π.The same statements apply to W because it has the same unique steady-state distribution.
For convenience, we define the matrix We write the Hermitian and anti-Hermitian parts of K as Note that L is Hermitian, thus λ L and λ W are real-valued.In addition, the second eigenvalue of L obeys the variational principle [60,Sec. 4.2] where 1. Derivation of first main result (7) We first plug the definitions of L and M from (C2) and (C3) into the definition of steady-state EPR, which appears as (3) in the main text.After a bit of rearranging, this gives The second line uses √ π i π j ≥ L ji /κ, where κ is defined in (5).The third line applies Jensen's inequality to the convex function tanh -1 x (for x ≥ 0) with weights Applying the Cauchy-Schwarz inequality then gives where we used η W from (6) and the Frobenious norm ∥ • ∥ F .We now combine (C6) and (C7), while using that x tanh -1 (y/x) is decreasing in x for x, y ≥ 0. After a bit of rearranging, this gives Next, we will use the following classic theorem by Kahan [39], which we state here without proof and using our own notation.

Theorem. For Hermitian
Observe that L + M = K, which follows from our definitions (C3).Moreover, M is anti-Hermitian, so ∥(M + M † )/2∥ F = 0 and ∥(M − M † )/2∥ F = ∥M ∥ F .Combining with the theorem above and rearranging gives In the last line, we used that Im λ L = 0 since L is Hermitian, and then used Eq.(C4).
The result (7) follows by combining Eqs.(C8) and (C9), while using that x tanh -1 x is decreasing in x ≥ 0. (14) Let u be the right eigenvector of K corresponding to λ K 2 , normalized so that ∥u∥ = 1.Using Ku = λ K 2 u and u † K T = (λ K 2 ) * u † , we can express the real part of the eigenvalue as

Derivation of the spectral gap inequality
Since K = D −1 W D and 1 T W = 0, K has a left eigenvector (1 T D) T = √ π with the corresponding eigenvalue λ K 1 = 0.The vector u is orthogonal to this eigenvector, u ⊥ √ π.Therefore, u satisfies the constraints of the maximization (C5).Combining (C4), (C5), and (C10) then gives Finally, to derive (14), note that W and W are rate matrices, thus λ W 2 ≤ 0 and Re λ W 2 ≤ 0.

Derivation of spectral gap inequalities
where we used that L ij = Wij π j /π i , as follows from our definitions (C1), (C2), and (C3).Next, we introduce the variables and the norm This allows us to rewrite where the last equality follows from 1 2 ij |ϕ i −ϕ j | 2 π i π j = 1 for any ϕ that satisfies i ϕ i π i = 0 and i |ϕ i | 2 π i = 1.This recovers the inequalities (15).

Figure 2 .
Figure 2. Thermodynamic bound on the spectral gap illustrated on the biomolecular sensing model from Fig. 1(a).Left: Two bounds from (16) (solid lines) versus the maximal spectral gap (dotted line) identified by numerical optimization at varying levels of steady-state EPR.Right: larger EPR is associated with faster relaxation to steady state, shown using the decay of KL divergence D(p(t)∥π) over time.

Figure 3 .
Figure 3. Thermodynamic bound on imaginary eigenvalues illustrated on the model from Fig. 1(a).Left: Two bounds from (19) (solid lines) versus the maximal imaginary part of the second eigenvalue (dotted line) identified by numerical optimization at varying levels of steady-state EPR.Right: larger EPR is associated with faster oscillatory "pulsing" during relaxation.

)
Introducing the diagonal matrix D ji = δ ji √ π i , we can express K and L as similarity transformations of W and W respectively, K = D −1 W D and L = D −1 W D. It then follows that they share the same spectra [60, Cor.1.3.4],