Information-thermodynamic principles unveil the response of complex networks to perturbations

Introduction

A defining feature of complex systems is not only their structure, but also their response to perturbations [1, 2, 3, 4].

Systems with similar unperturbed trajectories may react very differently to small disturbances, revealing latent structural and dynamical properties [5, 6]. Characterizing perturbation responses therefore provides a fundamental lens on the organization and function of complex networks, since perturbation propagation reveals how information, uncertainty, and functional variability are redistributed by the interplay between structure and dynamics [4, 5, 7, 8].

Standard diffusive dynamics offers a tractable class of network processes: its evolution is linear and governed by a time-independent, often symmetric, operator, under which perturbations evolve as well. Diffusion models phenomena ranging from consensus formation [9] to epidemic spreading [10, 11] and migration [12], and often approximates nonlinear dynamics near steady states [13]. Recently, a density matrix formalism [14], analogous to equilibrium statistical mechanics, enabled a thermodynamic description of networks coupled to diffusion. This diffusive information–thermodynamic framework has been applied to network classification [15, 16], phase transitions [17], centrality [18], robustness [19], functional fragility [20], structural sparsity [8], and renormalization group approaches [21, 22].

Despite these advances, it remains unclear whether diffusion provides an exhaustive description, or whether a more general thermodynamic framework is needed to capture behaviors beyond it. This is especially relevant for real networks, which often operate in dynamical regimes far from steady states. As non-equilibrium thermodynamics and statistical mechanics [23, 24, 25] describe phenomena inaccessible to equilibrium approaches, an analogous extension may be required for network responses governed by general time-dependent and non-normal operators [26].

In this article, we explore the theoretical limits of the diffusive approach commonly adopted in the network density matrix literature, which relies on a time-independent normal evolution operator and on maximum-entropy assumptions for the prior distribution of perturbations. In particular, we show that, under these common assumptions, the entropy derivative is constrained to be non-positive. This contrasts with real systems, which may maintain or increase functional diversity. We therefore construct a general information-thermodynamic framework capable of describing the response to perturbations at any point along a complex system trajectory and with a time-dependent non-normal evolution operator. This framework is based on an information-theoretic least dissipation principle, analogous in spirit to the Onsager–Rayleigh principle, and establishes a formal correspondence between non-equilibrium thermodynamics and network response theory.

Results

Network density matrices

The network density matrix formalism [14, 20] provides an information–thermodynamic framework to describe perturbation propagation in an $N$-dimensional complex network with adjacency matrix $\mathbf{A}$, whose state $\lvert x(t)\rangle$ follows the non-linear dynamics $\frac{d}{dt}\lvert x(t)\rangle= \mathbf{v}(\lvert x(t)\rangle;\mathbf{A})$. Perturbations are small shifts $\lvert\delta x(t)\rangle$ of the system state $\lvert x(t)\rangle$ (see Appendix A). Denoting by $\lvert i\rangle$ the canonical basis vector associated with node $i$, the response of individual nodes to a perturbation mode $|\delta x(t)\rangle$ is quantified through its projection onto this basis. The absolute response intensity is $|\delta x_i(t)|=\sqrt{\langle\delta x(t)|i\rangle^2}$, while the dynamical correlation between nodes $j$ and $i$ is $\delta x_j(t)\delta x_i(t)=\langle j|\delta x(t)\rangle\langle\delta x(t)|i\rangle$. Both are encoded in the local propagator $\mathbf{u}(t)=|\delta x(t)\rangle\langle\delta x(t)|$, which describes the spreading of a single perturbation mode.

Symbol	Dynamics	$\dot{x}_i=$
E	Epidemics (SIS)	$-Bx_i+R\sum_j A_{ij}(1-x_i)x_j$
NV	Noisy voter	$A-Bx_i+\frac{C}{k_i}\sum_j A_{ij}x_j$
$\mathbb{P}$	Population	$-Bx_i^b+R\sum_j A_{ij}x_j^a$
N	Neuronal	$-x_i+R\sum_{j=1}^{N}A_{ij}\tanh(x_j)$
$\mathbb{B}$	Biochemical	$F-Bx_i-R\sum_j A_{ij}x_ix_j$

A single perturbation generally cannot characterize the full network response, since different initial directions may produce distinct dynamics and the actual perturbation is often unknown. We therefore use an ensemble of $M$ possible initial perturbations, $\Psi=\{\delta x^\alpha(t_0)\}_{\alpha=0,\ldots,M-1}$, where each mode $\delta x^\alpha$ occurs at time $t_0$ with probability $p_\alpha$ and evolves under the Jacobian $\mathbf{J}=\frac{\partial\mathbf{v}}{\partial\mathbf{x}}$ (see Tab. I). To span all perturbation directions, we choose $M=N$ orthogonal modes.

The propagation of information through the network is then described by the statistical propagator $\mathbf{U}(t)=\sum_{\alpha=0}^{N-1}p_\alpha\mathbf{u}^\alpha(t)$, where $\mathbf{u}^\alpha(t)$ is the local propagator of mode $\alpha$. Its trace, $Z(t)\equiv\operatorname{Tr}(\mathbf{U}(t))$, equals the total weighted ensemble strength. The normalization of $\mathbf{U}(t)$ gives the network density matrix

This positive, unit-trace operator encodes the distribution of the perturbation ensemble across network nodes and modes.

Without prior information on $p_\alpha$, the maximum entropy principle [27, 28] assigns probability $1/N$ to each mode. A standard choice for $\Psi$ is the canonical basis, with perturbations initially localized on single nodes [8, 14, 17, 20, 21]. To isolate the effect of the evolution operator, throughout the main text we adopt this maximum-entropy prior. for the initial perturbation ensemble, $\rho(t_0) = I/N$, unless otherwise specified. This keeps the initial ensemble fixed when comparing different evolution operators and avoids introducing an additional directional bias through the perturbation prior. The role of non-uniform priors, including their possible effect on the sign of the entropy rate, is discussed in Appendix § 1.5.

For diffusive dynamics, this framework is formally connected to equilibrium statistical mechanics and the Gibbs ensemble of densities [29] (see Appendix B). Moreover, by selecting this uniform prior for the initial perturbation distribution, the density matrix formally coincides with the normalized left Cauchy–Green tensor of classical dynamical systems theory [30]. Therefore, the density-matrix representation does not introduce additional microscopic information with respect to tangent-space dynamics, but changes the descriptive level by recasting the collective response to perturbations as a statistical state. This representation makes it possible to summarize high-dimensional response dynamics through a small set of macroscopic thermodynamic observables and to apply statistical and variational tools to characterize the joint effects of network structure and dynamics across scales.

The density matrix flow equation and entropy evolution

The first step toward a response theory for complex interconnected systems along general trajectories is to understand how the evolution deforms the density operator $\rho(t)$ describing response modes. For perturbations evolving under a time-dependent, generally non-symmetric $\mathbf{J}(\mathbf{x}(t))$, differentiating the density matrix (Appendix C) yields the equation governing its flow on the manifold of all admissible density matrices:

The first two terms describe the local deformation of $\rho$ induced by $\mathbf{J}(\mathbf{x}(t))$, including stretching and directional changes of perturbation modes. The third term preserves $\operatorname{Tr}(\rho) = 1$ and thus contains $\operatorname{Tr}(\rho\mathbf{J})$, linked to the local stretching rate of an infinitesimal volume around the trajectory.

We then exploit the flow equation to explore the role of symmetric and antisymmetric dynamical features [31].

Decomposing $\mathbf{J}$ into its symmetric part $\mathbf{D}$ and antisymmetric part $\mathbf{C}$ gives

where $\{\cdot,\cdot\}$ and $[\cdot,\cdot]$ denote anticommutator and commutator. Eigenvalue and eigenvector perturbation theory [32] (Appendix C) shows that $\mathbf{D}$ affects both eigenvalues and eigenvectors of $\rho$, while $\mathbf{C}$ only rotates eigenvectors unitarily. Thus, $\mathbf{D}$ represents the dissipative part of the evolution, changing the relative relevance of perturbations, whereas $\mathbf{C}$ is conservative. Consistently, only $\mathbf{D}$ enters the normalization term, since $\operatorname{Tr}(\rho\mathbf{J}) = \operatorname{Tr}(\rho\mathbf{D})$. The diffusive case is recovered by setting $\mathbf{C} = 0$ and assuming time independence of $\mathbf{D}$.

The conservative part can always be removed by rotating $\rho(t)$ and $\mathbf{D}(t)$ through $\mathbf{O}(t,t_0) = \mathcal{T}\exp\left[\int_{t_0}^{t} \mathrm{d}s\,\mathbf{C}(s)\right]$, where $\mathcal{T}$ is the time-ordered exponential. In this dissipative frame, $\sigma(t) \equiv\mathbf{O}^{-1}\rho(t)\mathbf{O}$ evolves only under $\mathbf{Q}(t) \equiv\mathbf{O}^{-1}\mathbf{D}(t)\mathbf{O}$:

This rotation preserves relevant thermodynamic quantities such as the Von Neumann entropy $S$ and $\operatorname{Tr}(\rho\mathbf{J})$, since they are invariant under unitary transformations. Despite this formal mapping, the time dependence and the non-normality of the evolution operator in the general case mark a fundamental difference from the time-independent and symmetric diffusive case.

We exploit the flow equation to study the evolution of the Von Neumann entropy $S = -\operatorname{Tr}(\rho\log\rho)$, which measures the functional diversity of the response through the relative importance of its modes. Our goal is to identify differences from the time-independent diffusive case. Substituting the flow equation into the entropy derivative yields

where $\operatorname{Cov}_{\rho}(\cdot,\cdot)$ denotes the covariance with respect to $\rho$ [33]. By the Cauchy–Schwarz inequality, we then obtain

This is an upper bound on the entropy rate [34] and sets a scale-dependent limit to changes in functional diversity.

Given a scale $t$, we define the instantaneous Gibbsian reference as $\gamma(t)=\frac{e^{2t\mathbf{D}(t)}}{\operatorname{Tr}(e^{2t\mathbf{D}(t)})}$, obtained by freezing the symmetric operator $\mathbf{D}(t)$ at that time. We call $\gamma(t)$ the instantaneous diffusive reference, interpreting it as the effective diffusive process locally approximating the real dynamics. The bound in Eq. (6) is saturated when $\log\rho=a\mathbf{D}+b\mathbf{I}$, with $a,b\in\mathbb{R}$. As a special case, this occurs when the actual density operator coincides with $\gamma(t)$ at the corresponding scale.

Aside from the comparison with its instantaneous diffusive counterpart, we also contrast the entropy derivative of a system evolving under non-linear dynamics with that of a network with the same adjacency matrix evolving under Laplacian diffusion (Figure 1–a).

they do not seem to involve features that are, in principle, theoretically incompatible with a time-independent and symmetric diffusive approximation.

However, one phenomenon cannot be captured in such a framework: entropy increase, i.e. a positive entropy derivative $\dot{S} > 0$. In real systems, we observe phases where the diversity of responses to perturbations increases, due to both structural and dynamical effects. In the time-independent and symmetric diffusive case, instead, the entropy derivative is always non-positive. Indeed, for a constant and symmetric evolution operator $\mathbf{D}$ and a maximum–entropy prior, Eq. (5) reduces to $\dot{S} = -4t\operatorname{Var}_{\rho}(\mathbf{D})$; since the variance is non-negative, $\dot{S} \le0$ (Appendix § 1.4). This is analogous to the Boltzmann H-theorem [37], but with an inverted sign due to the different interpretation of entropy in our context, where we start from a maximum–entropy state. This bound has a clear physical interpretation: if the evolution operator is constant and symmetric and we assume a maximum–entropy prior, the diversity of the response to perturbations of a complex network cannot increase in time. For this reason, a non-diffusive information-thermodynamic framework is required to fully explain the response to perturbations along a general dynamical trajectory, where a non-uniform perturbation prior, time dependence, or non-normality of the evolution operator may each yield $\dot{S} > 0$. Indeed, Eq. (5) shows that the sign of the entropy rate is not determined by the instantaneous symmetric operator $\mathbf{D}$ alone, but by its relation to the current density operator, which contains the accumulated effect of both the initial perturbation ensemble and the preceding dynamics. Let $\rho= \sum_i p_i |p_i\rangle\langle p_i|$ and define $d_i^{(\rho)} \equiv\langle p_i|\mathbf{D}|p_i\rangle$ and $\langle\mathbf{D}\rangle= \operatorname{Tr}(\mathbf{D}\rho)$, where $d_i^{(\rho)}$ is the instantaneous stretching rate along the ith response direction. Then, equation (5) can be rewritten as

Hence, the sign of $\dot{S}$ is controlled by the correlation between the accumulated relevance of each response mode, encoded by $\log p_i$ and shaped by both the initial prior and the previous evolution along the trajectory, and its instantaneous stretching, represented by $d_i^{(\rho)}$. In particular, $\dot{S} > 0$ when these two quantities are negatively correlated. Physically, this occurs when modes that were previously weakly represented undergo comparatively stronger instantaneous stretching, thereby redistributing the response toward previously suppressed directions. This mech...anism is absent in the standard diffusive case, for which $\log\rho= 2tD - (\log Z)\,I$, so that mode relevance remains locked to the fixed stretching rates and cannot become anticorrelated with them, whereas such an anticorrelation may occur whenever this relation no longer holds (Appendices § 1.5– § 1.6– § 1.7).

To test whether this thermodynamic signature can be associated with concrete dynamical features, we consider the C. elegans neuronal network with nonlinear neuronal dynamics (N). In the low-activity regime, the system has a linear-stability threshold around $R_c \simeq0.1$, above which the Jacobian can acquire eigenvalues with positive real part, signaling the onset of locally unstable modes in the linearized dynamics. We therefore examine whether the occurrence of $\dot{S} > 0$ changes across this threshold. In this system, the total time spent with $\dot{S} > 0$ becomes appreciable above $R_c$ and increases at stronger coupling, whereas it remains identically zero for Laplacian diffusion (Figure 1b). This association does not imply a general one-to-one correspondence between local instability and positive entropy production, whose sign depends on the instantaneous correlation between the response density and the symmetric part of the Jacobian. Rather, it shows that, in this nonlinear system, $\dot{S} > 0$ captures a change in perturbation-response organization occurring across a dynamical stability threshold that the diffusive description cannot reproduce.

Informational least-dissipation principle

Because of the intrinsic limits of the diffusive approach, we develop an information-thermodynamic framework to describe responses to perturbations along general trajectories. We provide a theoretical foundation for the relevant thermodynamic quantities, working for simplicity in the fully dissipative frame; the conservative part can be restored by an appropriate rotation, and the fundamental thermodynamic quantities are invariant under this rotation.

The evolution of all thermodynamic quantities is determined by the flow of the density operator $\sigma$. Thus, we seek a foundational principle from which it can be derived (Appendix § 1.8). For this purpose, we consider the manifold of density matrices endowed with the Bures metric [38, 39], defining the infinitesimal distance between $\sigma$ and $\sigma+ \mathrm{d}\sigma$:

where $\mathbf{K}_{\sigma}(\mathrm{d}\sigma)$ satisfies $\sigma\mathbf{K}+\mathbf{K}\sigma=\mathrm{d}\sigma$. We interpret this metric as the informational cost of distinguishing neighboring density operators. Unit trace preservation restricts admissible variations to traceless $\mathrm{d}\sigma$, equivalently $g_{\sigma}(\sigma,\mathrm{d}\sigma)=0$. We assume that the system’s intrinsic evolution drives its response to perturbations. Accordingly, we introduce the potential $U(\sigma)=-2\operatorname{Tr}(\mathbf{Q}\sigma)$, measuring alignment between the intrinsic evolution operator $\mathbf{Q}$ and the density operator. This generalizes the energy-like role that the Laplacian $\mathbf{L}$ has in the diffusive limit. We now state the following variational principle: at any scale $t$ associated with a driving $\mathbf{Q}(t)$, the density operator follows the trace-preserving path with the least informational cost that maximizes the instantaneous alignment with $\mathbf{Q}(t)$. This is obtained by minimizing

under $g_{\sigma}(\sigma,\dot{\sigma})=0$, enforced by the Lagrange multiplier $\alpha$. Imposing $\frac{\delta}{\delta\dot{\sigma}}\mathcal{R}=0$ recovers the density matrix flow equation. The principle determines the dissipative structure, while the conservative part is a gauge choice recoverable by rotating back to the original frame.

This is an informational least dissipation principle, based on the cost of distinguishing neighboring density operators. It is analogous to the Onsager–Rayleigh least dissipation principle [24, 25], and $\mathcal{R}$ has the same formal structure as the Rayleigh dissipation functional [40]. Moreover, the variational formulation rewrites $\dot{\sigma}$ as a gradient flow driven by the gradient of the time-dependent potential $U$ with respect to the Bures metric (Appendix I). This supports the extension of $U(t)=-2\operatorname{Tr}(\mathbf{J}\rho)$ as a thermodynamic potential governing the density flow, interpreted as the average contraction ($U>0$) or expansion rate ($U<0$) of perturbations (Figure 2).

*Generalized free energy—*After the introduction of the generalized framework, we examine its relation to the diffusive approximation at scale $t$. Following non-equilibrium thermodynamics [24, 25], we combine the potential $U$ with the Von Neumann entropy into the generalized free energy $\mathcal{F}\equiv U-\frac{1}{t}S$. This yields a measure of the distance from the instantaneous diffusive reference $\gamma(t)$, quantifying how well the local perturbation response is approximated by diffusion. Indeed, with some algebra,

where $F_{\gamma}=-\frac{1}{t}\log[\operatorname{Tr}(e^{2tD})]$ is the Gibbsian free energy of $\gamma$ and $D_{KL}(\rho\|\gamma)$ is the Kullback–Leibler divergence between $\rho$ and $\gamma$. Thus, $D_{KL}$ emerges as a thermodynamically motivated measure of the distance between a system density operator and its Gibbsian approximation.

We evaluate this divergence on real-world networks with appropriate dynamics, comparing them with configuration-model and Erdős–Rényi counterparts with the same number of nodes and comparable average degree (Figure 2). Two distinct comparisons are involved. For each structure–dynamics pair, $D_{KL}(\rho\|\gamma)$ measures the departure of the actual perturbation response from its own instantaneous diffusive reference. These distances can then be compared across network ensembles. Although no strict ordering holds at every time, we observe a recurrent tendency for empirical networks to depart more strongly from their diffusive references than configuration-model networks, and for the latter to depart more strongly than Erdős–Rényi networks.

This tendency admits a natural information-theoretic interpretation. In maximum-entropy approaches, complexity can be quantified by the information not captured by reconstructions based on a restricted set of constraints [41]. The three network classes considered here define an analogous hierarchy of structural information: Erdős–Rényi networks retain essentially only network size and mean density, the configuration model additionally preserves the degree sequence, while empirical networks retain correlations and higher-order organization removed by randomization. Our framework provides a thermodynamically motivated dynamical counterpart of such complexity measures. Since the response state $\rho(t)$ is generated by $\mathbf{J}[\mathbf{x}(t);\mathbf{A}]$, structural information retained in $\mathbf{A}$ can be reflected in the organization of the perturbation response. Accordingly, $D_{KL}(\rho\|\gamma)$ quantifies the response information that is not reducible to the corresponding instantaneous Gibbsian approximation. The recurrent hierarchy in Fig. 2 therefore suggests that progressively richer structural organization is accompanied, in the systems considered here, by progressively less reducible dynamical responses. We stress that this is an information-theoretic interpretation of the observed statistical tendency rather than as a universal pointwise inequality holding for every realization and time.

Conclusions and outlook

Real-world complex networks exhibit heterogeneous responses to perturbations. Using information-thermodynamic arguments, we showed that standard diffusive approaches with time-independent evolution operators cannot capture this richness, as they fail to describe transient increases in response diversity. This motivates a trajectory-based formalism covering broad dynamical regimes, including far from steady state ones.

We developed an information-thermodynamic framework based on a least dissipation variational principle, establishing a correspondence with Onsager non-equilibrium thermodynamics analogous to the link between the diffusive-like case and equilibrium thermodynamics. This paves the way for the extension of many tools previously applied successfully in the diffusive-like case to a broader dynamical landscape. Moreover, current renormalization group techniques for complex networks [21, 22] are based on the preservation of the dissipative structure of Laplacian diffusion across time scales. Our formalism manages to encode the dissipative structure of many dynamics at any scale $t$ in $\mathcal{R}(t)$, opening the possibility for a network renormalization group that incorporates both structural and dynamical features governing information spreading: a key step toward a general renormalization theory for interacting complex systems.

The normalization-preserving factor in the flow equation is consistent with analogous terms in probability-density equations defined on systems whose intrinsic dynamics are not probability preserving, such as the classical Liouville equation for non-Hamiltonian systems [23] or the Lindblad equation for open quantum systems [44]. A structure combining intrinsic evolution with a normalization term linear in the distribution and proportional to the generator average is also typical of replicator dynamics [45]. Thus, equation (2) can be viewed as a generalized matrix replicator equation [46], where $\mathbf{J}(\mathbf{x}(t))$ acts as a fitness matrix for $\rho$, aligning perturbation ensembles with instantaneous evolution directions. These links may connect the framework to Bayesian inference [47] and support new network inference techniques.

Author contributions. A.F.B. and M.D.D. designed the study, performed the theoretical analysis and wrote the manuscript. A.F.B. performed the numerical experiments.

Competing interests. The authors declare no competing interests.

Appendices

Response to perturbations of a dynamical system

We consider a complex system composed of a set of $N$ interacting units, or nodes. At any time $t$, each node $i$ is associated with a state variable $x_i(t)$, and the collective state of the system is represented by the $N$-dimensional vector $\lvert\mathbf{x}(t)\rangle$. The time evolution of the system is determined by both its structure and its dynamics. Structural interactions are encoded in an adjacency matrix $\mathbf{A}$, whose entry $A_{ij}$ represents the strength of the interaction from node $j$ to node $i$, while the dynamics are specified by a generally nonlinear vector field $\mathbf{v}(\lvert x(t)\rangle;\mathbf{A})$.

Throughout this letter, we focus on systems governed by first-order differential equations, $\frac{d}{dt}\lvert x(t)\rangle=\mathbf{v}(\lvert x(t)\rangle;\mathbf{A})$.

To probe the system’s sensitivity to perturbations, we consider a small displacement $\lvert\delta x(t)\rangle$ of the state vector.

Linearizing the dynamics around the trajectory $\lvert x(t)\rangle$, the evolution of the perturbation is governed by the Jacobian matrix $\mathbf{J}$, with elements $J_{ij}=\frac{\partial v_i}{\partial x_j}$, according to

The Jacobian is evaluated along the trajectory generated by $\mathbf{v}(\lvert x(t)\rangle;\mathbf{A})$ together with the chosen initial condition.

Our aim is to characterize the response to perturbations of a complex system along any trajectory generated by $\mathbf{v}(\lvert x(t)\rangle;\mathbf{A})$, given an initial condition $\lvert x(t_0)\rangle=\lvert x_0\rangle$.

Correspondence between the diffusive case and equilibrium statistical mechanics

For diffusive dynamics, the evolution of the system is determined by a linear vector field, $\mathbf{v}(\lvert x(t)\rangle;\mathbf{A})=-\mathbf{L}\lvert x(t)\rangle$, where, in the common case of Laplacian diffusion, $\mathbf{L}$ is a time-independent Laplacian matrix [48], with entries $L_{ij}=k_i\delta_{ij}-A_{ij}$ and $k_i$ denoting the degree of node $i$. In the case of Laplacian diffusion on an undirected network, equation (A.1) can then be solved with constant $\mathbf{J}=-\mathbf{L}$, yielding the density matrix at any propagation scale $t$ as

Formally, the density matrix defined in equation (B.1) is equivalent to a canonical Gibbs state associated with the network Laplacian [49], with $\mathbf{L}$ playing the role of a Hamiltonian-like operator and $t$ acting as a scale parameter, as the classical inverse temperature $\beta$ of equilibrium statistical mechanics. This correspondence establishes a direct link between diffusion on networks and equilibrium statistical mechanics, enabling the construction of an information–thermodynamic framework for diffusive dynamics [29]. In this framework, $\rho$ acts as a density operator and $Z$ plays the role of a partition function, from which a free energy $F=-\frac{1}{t}\log(Z)$ can be defined.

This quantity has been interpreted as a measure of the speed of information propagation across the network [50].

The associated Von Neumann entropy, $S=-\operatorname{Tr}(\rho\log\rho)$, quantifies how evenly perturbation modes contribute to the system’s response. High entropy corresponds to a broad distribution of contributing modes, while low entropy indicates the dominance of a small number of modes. The entropy is maximized when all modes contribute equally, yielding $S=\log N$, which provides a natural normalization such that the normalized entropy satisfies $S\in[0,1]$. Within this interpretation, the entropy captures the functional diversity of the system’s response to perturbations.

The free energy $F$, entropy $S$, and Laplacian $\mathbf{L}$ satisfy the information-theoretic relation [51]

where $U=2\operatorname{Tr}(\mathbf{L}\rho)$ defines an average internal energy consistent with the role of $-\mathbf{L}$ as an evolution operator.

At the scale $t=0$, this energy has a clear physical interpretation as it is proportional to the number of edges in the network.

Density flow equation

Given an ensemble of orthogonal perturbations $\Psi_t=\{\delta x^\alpha(t)\}_{\alpha=0,\ldots,N-1}$, we build the associated density matrix as

where the $\{p_\alpha\}_{\alpha=0,\ldots,N-1}$ are the initial probabilities assigned to each perturbation mode and

By using (A.1), we derive that

and

Then, exploiting (C.3) and (C.4), we directly recover the flow equation by applying the time derivative to the definition in (C.1):

Now, consider an ensemble of perturbations $\Psi_{t_0}=\{\delta x^{\alpha}(t_0)\}_{\alpha=0,\ldots,N-1}$ that is an orthonormal basis at time $t_0$. We associate to it a density operator $\rho(t_0)$, whose evolution is given by equation (2). Given the initial perturbation density operator $\rho(t_0)$, the solution of equation (2) is

in which $\Phi(t,t_0)=\mathcal{T}\left\{\exp\left[\int_{t_0}^{t} \mathrm{d}s\,J(s)\right]\right\}$ is the time-ordered exponential built from the Jacobian, and $\Phi^{T}$ is the anti-time-ordered exponential built from its transpose. This general solution yields a density matrix that is a matrix functional $\rho[x(t);t]$ of the trajectory $x(t)$ and allows for the direct extension of many of the applications previously cited for the diffusive case to the general one.

We want to investigate how the symmetric part $D$ and the antisymmetric part $C$ of $J$ influence the evolution of the eigenvalues and eigenvectors of (C.6). We begin by substituting this separation into the flow equation, obtaining

According to Rayleigh’s matrix perturbation formulas, given a matrix $M$ with starting eigenvalues $m_i$ and associated eigenvectors $|m_i\rangle$, the shifted matrix $B=M+\delta M$ has eigenvalues

and eigenvectors

for a small enough shift $\delta M$. We can apply these formulas to the density matrix evolution considering the shift $\delta\rho= \delta t \frac{d\rho}{dt}$, and substitute the explicit form of the flow equation. By doing so and indicating with $p_i(t)$ and $\lvert p_i(t)\rangle$ respectively the eigenvalues and the eigenvectors of $\rho$ at time $t$, we get for the eigenvalues

and, for the eigenvectors,

This shows that the symmetric part $D$ affects both the eigenvalues and eigenvectors, while the antisymmetric part $C$ only acts on the eigenvectors, generating a unitary rotational evolution.

Network H-theorem

In this section, we characterize the conditions under which the response entropy is guaranteed to be non-increasing. For a density operator $\rho(t)$ evolving according to the previously discussed flow equation, the entropy rate is

Therefore, at any given time, the general instantaneous condition to have a non-increasing entropy is

In particular, the sign of $\dot{S}$ is not determined by the instantaneous symmetric operator $D(t)$ alone, but by its relation to the current density operator $\rho(t)$, which contains the accumulated effect of both the initial perturbation ensemble and the preceding dynamics.

The physical content of the H-theorem lies in identifying classes of evolution for which the features of the underlying dynamics and of the chosen prior always forbid an increase of entropy. A relevant class of this kind is represented by dynamics with a time-independent normal evolution operator $\mathbf{J}$ starting from a maximum entropy prior $\rho(0) = \mathbf{I}/N$. Indeed, in this case,

hence,

and Eq. (D.1) gives

The time-independent symmetric diffusive case discussed in the main text corresponds to the special case $\mathbf{J} = \mathbf{D}$.

By imposing a few additional constraints on this class, the theorem can be extended to dynamics governed by a symmetric time-dependent evolution operator. In particular, we require a symmetric time-dependent operator $\mathbf{D}(t)$ satisfying

and assume that the ordering of its eigenvalues,

is preserved throughout the interval of interest. In this case, starting from $\rho(0) = \mathbf{I}/N$ and defining the accumulated stretching rates as

the eigenvalues of the density operator are

Since the ordering of the instantaneous stretching rates is preserved, it follows that

Writing the covariance in equation (D.1) in the common eigenbasis of $\mathbf{D}(t)$ and $\boldsymbol{\rho}(t)$, and using $\log p_i(t)=2A_i(t)-\log Z(t)$, we obtain

Using $\sum_i p_i=1$ and symmetrizing, this becomes

Indeed, each term in the sum is non-negative because the instantaneous and accumulated stretching rates have the same ordering. Therefore, the entropy remains non-increasing throughout any interval over which this ordering is preserved.

Aside from the theorem’s range of applicability, it is important to investigate its limitations and the mechanisms that may induce $\dot{S}(t)>0$. In the following sections, we therefore discuss explicitly how entropy can increase when the assumptions on the initial perturbation ensemble or the evolution operator are relaxed. Specifically, we start from the time-independent, symmetric diffusive case with a maximum-entropy prior and examine separately the effects of introducing a biased prior, non-normality, and time dependence, showing through examples how each of these elements can induce time intervals with $\dot{S}>0$. Indeed, the examples will identify three distinct mechanisms that influence the sign of the entropy derivative. A non-uniform prior changes the initial relation between the density operator and the stretching directions even for ordinary diffusion. Time dependence can reverse the accumulated concentration of response modes even when the instantaneous evolution operator is symmetric. Non-normality can dynamically rotate response directions relative to the corresponding symmetric stretching operator even for a constant evolution generator.

Dependence on the perturbation prior

The maximum-entropy prior is essential to the monotonicity result that we obtain in the main text for the standard diffusive case. To show this, we consider the time-independent symmetric diffusion operator $\mathbf{J}=\mathbf{L}$, but choose the non-uniform initial density

Since $\rho(0)$ commutes with $\mathbf{L}$, its evolution is

Therefore,

Since $\mathbf{D}=-\mathbf{L}$, Eq. (D.1) becomes

For every nontrivial $\mathbf{L}$, it follows that

At $t=c/2$, $\rho=\mathbf{I}/N$ and the entropy reaches its maximum, after which the usual diffusive entropy decrease is recovered.

Thus, diffusion alone does not determine the sign of the entropy rate. A biased perturbation prior may initially evolve toward a more uniform response distribution and therefore produce $\dot{S}>0$ even under a constant symmetric diffusion operator. The maximum-entropy prior adopted in the main text removes this independent source of directional bias and allows us to isolate the effect of the evolution operator.

Time dependence

Time dependence provides a distinct mechanism for breaking entropy monotonicity, even when the evolution operator is symmetric and commutes with itself at different times. To see this, consider

The generator is symmetric for every $t$, and

Starting again from $\rho(0)=\mathbf{I}/2$, define

The density operator becomes

where $\sigma_z=\operatorname{diag}(1,-1)$. Its entropy is

and therefore

Since $A(t)>0$ for $0<t<1$,

whereas

At $t = 1$, $A(1) = 0$ and the maximum-entropy state is recovered.

This example also clarifies the conditions required for a time-dependent extension of the $H$-theorem. Pairwise commutation and a common eigenbasis are not sufficient if the ordering of the instantaneous stretching rates changes during the evolution.

Non-normality of the evolution operator

Time dependence of the evolution operator and a biased prior are not required to obtain a positive entropy derivative; indeed, non-normality can also generate intervals with $\dot{S} > 0$. To isolate its effect, consider the constant two-dimensional linear evolution operator

For $\alpha\ne1$,

so that $\mathbf{J}$ is non-normal. The optional isotropic contribution $-\gamma\mathbf{I}$, with $\gamma> 0$, makes the linear dynamics asymptotically stable, but cancels out of the normalized density operator and therefore does not affect $S(t)$.

Since

we have

Starting from the maximum-entropy state $\rho(0) = \mathbf{I}/2$, let

The density operator is then

The determinant of $\boldsymbol{\Phi}(t)\boldsymbol{\Phi}^{\mathsf{T}}(t)$ is unity, while its trace is

from which it follows that $\tau(t) \ge2$. The two eigenvalues of $\rho(t)$ are therefore

and the entropy is

Using the chain rule, we write the time derivative of the entropy as

and observe that its sign can change only through the factor $\dot{\tau}$. Indeed, whenever $\dot{S} \ne0$, and hence $\tau> 2$, we have

and

Thus, since

the entropy decreases for $0 < t < \pi/2$, but increases for

At $t=\pi$, $\boldsymbol{\Phi}=-\mathbf{I}$, and the maximum-entropy state $\rho=\mathbf{I}/2$ is recovered.

This example shows that a system with a constant evolution operator and starting from the same maximum-entropy prior used in the main text can exhibit finite intervals of increasing response entropy solely as a consequence of non-normality. In terms of the decomposition $\mathbf{J}=\mathbf{D}+\mathbf{C}$, the antisymmetric contribution $\mathbf{C}$ does not directly change the eigenvalues of $\rho$ at a given instant. For a non-normal generator, however, $[\mathbf{D},\mathbf{C}] \ne0$, and the rotation induced by $\mathbf{C}$ continuously changes the orientation of the response eigenvectors relative to the stretching directions of $\mathbf{D}$. This modifies their subsequent evolution and can reverse the sign of the covariance in Eq. (D.1).

Derivation of the flow equation from the least dissipation principle

We consider the manifold $\mathcal{D}$ of full-rank density operators $\sigma$ with $\operatorname{Tr}\sigma= 1$, endowed with the Bures metric. Given a tangent vector $X$ at $\sigma$, we define the operator $K_{\sigma}(X)$ as the unique solution of the equation

Given two tangent vectors $X$ and $Y$, the Bures metric is then defined as:

which is symmetric and positive definite on the space of trace-preserving variations. Since $\operatorname{Tr}\sigma= 1$ at any scale $t$, admissible velocities satisfy $\operatorname{Tr}\dot{\sigma} = 0$. Using (H.1) with $X = \sigma$, we find $K_{\sigma}(\sigma) = \frac{1}{2}I$, hence we write this constraint as

Moreover, varying the tangent vector at fixed $\sigma$, we have

Then, exploiting (H.1), we write

Applying the cyclic property of the trace to the second term of the right hand side of this equation, we get to

Subsequently, noticing that $\sigma\delta K_\sigma(X) + \delta K_\sigma(X)\sigma= \delta X$ and considering together (H.5) and (H.6), we obtain

which yields $\delta g_\sigma(X,X) = 2\operatorname{Tr}(K_\sigma(X)\delta X)$ and thus

We introduce the informational potential

and the dissipation functional

where the term with $\alpha$ enforces the constraint (H.3) and $\dot{U} = -2\operatorname{Tr}(Q\dot{\sigma}) - 2\operatorname{Tr}(\dot{Q}\sigma)$. Varying (H.10) with respect to $\dot{\sigma}$ at fixed $(\sigma,t)$ and then imposing $\frac{\delta}{\delta\dot{\sigma}}\mathcal{R} = 0$, yields

or equivalently

Using (H.1) with $X = \dot{\sigma}$ and (H.12),

Imposing $\operatorname{Tr}\dot{\sigma} = 0$ fixes $\alpha= 4\operatorname{Tr}(Q\sigma)$, leading to

Density flow equation as a gradient flow

The differential of the time-dependent informational potential $U$ along a traceless direction $\mathbf{X}$ is

On the subspace of traceless directions, the Riemannian gradient $\nabla_{g}U$ is defined by

This is satisfied by

which yields

Therefore, the flow (H.14) can be written as

showing that the density operator follows a gradient flow with respect to the Bures metric with a time-dependent potential $U$.

References and Notes