Stochastic Counterdiabatic Driving via Biorthogonal Liouvillian Eigenmodes

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Stochastic Counterdiabatic Driving via Biorthogonal Liouvillian Eigenmodes

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Stochastic Counterdiabatic Driving via Biorthogonal Liouvillian Eigenmodes Sandeep Suresh Cranganore,1 Sebastian Lehner,1 Johannes Brandstetter,1, 2 and Max Welling3, 4

arXiv:2607.24393v1 [physics.comp-ph] 27 Jul 2026

Institute for Machine Learning (Ellis Unit), Johannes Kepler University, Linz 2 Mistral AI, Paris, France 3 AMLab, Informatics Institute, University of Amsterdam, The Netherlands 4 CuspAI, Cambridge, England (Dated: July 28, 2026)

Finite-time driving of stochastic systems generates excess dissipation, causing the evolving probability distribution to lag behind the instantaneous equilibrium, and consequently degrading the convergence of nonequilibrium free energy estimators based on the Jarzynski equality. Escorted free energy simulations address the non-adiabatic lag by engineering control fields u that eliminate the lag, enforcing the trajectory-wise equality Wu = ∆F, and yielding zero-variance estimators. However, constructing the escorting field in closed form remains a challenge, approached variously through flow-field methods, targeted free energy perturbation, or learned diffeomorphisms. In this work, we construct a complementary numerical framework based on gauge-type transforms instead of generalized coordinate transforms for perfect escorting based on the exact spectral decomposition of the time-dependent Fokker–Planck generator. The biorthogonal decomposition of the Liouville operator directly yields a counterdiabatic correction whose action on the instantaneous equilibrium distribution exactly cancels the non-adiabatic lag at arbitrary driving speed in formal analogy with shortcuts-to-adiabaticity techniques such as Berry’s transitionless driving for quantum systems. Numerical verification for simulations of an overdamped particle in a time-varying double-well potential and harmonic traps confirms that the counterdiabatic condition is satisfied to machine precision, with the non-adiabatic lag suppressed by roughly twelve orders of magnitude in total variation distance and sixteen orders in KL divergence relative to the unescorted dynamics. As a diagnostic, we demonstrate vanishing dissipated work Wdiss (t) ≈ 0 for the deterministically propagated Fokker–Planck density across all protocol speeds.

INTRODUCTION

The adiabatic theorem underpins a wide range of phenomena in physics, from quantum state preparation [1] to classical non-equilibrium thermodynamics [2]. In both settings, adiabatic following fails whenever there is a finite-time switching, i.e., when protocol speed is finite: a quantum system leaks into excited states, while a classical diffusing particle lags behind the instantaneous equilibrium distribution. To systematically quantify this lag and relate microscopic fluctuations to macroscopic laws, stochastic thermodynamics has emerged as a foundational framework.This microscopic approach has since found widespread applications across multiple disciplines, ranging from molecular dynamics and soft-matter design to nonperturbative quantum chromodynamics [3]. Crucially, it provided the arena for the groundbreaking work of Jarzynski [4], which marked a major breakthrough in non-equilibrium physics. Later identified as an instance of annealed importance sampling (AIS) [5], this framework allows for the estimation of equilibrium free energy differences from non-equilibrium work via the celebrated Jarzynski Equality (JE): ⟨e−βW ⟩ = e−β∆F ,

(1)

where β = 1/(kB T ) is the inverse temperature. By Jensen’s inequality applied to the convex exponential, ⟨e−βW ⟩ ≥ e−β⟨W⟩ , Eq. (1) immediately yields ⟨W⟩ ≥ ∆F, recovering the second law for isothermal processes: the mean dissipated work ⟨Wdiss ⟩ = ⟨W⟩ − ∆F ≥ 0 is

non-negative, with equality achieved only in the quasistatic limit τ → ∞. The JE is an instance of a broader class of exact nonequilibrium relations known as fluctuation theorems [6]. In particular, the Crooks fluctuation theorem [6] relates the probability distributions of work in the forward and reverse protocols, PF (W)/PR (−W) = eβ(W−∆F ) , from which the Jarzynski equality follows by integration over the forward work distribution. These relations hold arbitrarily far from equilibrium and for any protocol duration τ . However, the practical utility of Eq. (1) as a free energy estimator is severely bottlenecked by its convergence properties under finite-time protocols. Driven by a time-varying external potential V (x, ζt ) with a control parameter ζt ≡ ζ(t), the system develops a structural ”lag” relative to the instantaneous Boltzmann distribution πeq (x, ζt ) = exp[−β(V (x, ζt ) − F (ζt ))] due to the non-adiabatic nature of the driving. These accumulated diabatic errors scale with the protocol speed, driving the system far from equilibrium and manifesting macroscopically as an excess dissipated work Wdiss (x, ζt ). Consequently, the exponential average in the JE becomes heavily dominated by rare, microscopic trajectories that manage to counteract this dissipation, leading to severe sample variance and poor statistical convergence. To mitigate this limitation, a variety of seminal strategies have been developed to enhance convergence and optimize free energy estimation. The idea of eliminating dissipation via an escorting control field was first derived by Vaikuntanathan and Jarzynski [7], who showed that a flow field u(x, t) can

2 be engineered that can reduce or even purge this lag entirely, yielding the escorted version of Eq. (1), which is often called the escorted Jarzynski equality (EJE): −βWu

⟨e

−β∆F

⟩u = e

(2)

where, the conventional work is modified by the contributions from the escorting term: Wu (x, ζt ) := W(x, ζt )+u(x, t)·∇V (x, ζt )−β −1 ∇·u(x, t) , (3) and the non-equilibrium work is defined as W(x, ζ ) = t D E Rt ′ ∂V (x,ζt′ ) ′ . Under perfect escorting u(x, t), dt ζ̇ t ∂ζ 0 ρ(t′ )

the additional terms in Eq. (3) ensure that the system always remains in equilibrium, with respect to V (x, ζt ). Thus, for every trajectory, it holds that Wu = ∆F, resulting in a zero variance estimate of the free energy difference [7, 8]. Thus, one can arbitrarily accelerate the convergence relative to Eq. (1), by designing the escorting fields such that the system remains near equilibrium throughout the finite-time process. The strategy of engineering auxiliary escorting fields to enforce tracking of a target distribution – as if the system were evolving adiabatically – falls under the overarching framework of shortcuts to adiabaticity (STA) [9, 10]. Historically, the concept of a shortcut was pioneered in the context of steering quantum states, introduced through the seminal frameworks of transitionless quantum driving by Berry [11] and assisted adiabatic passage by Demirplak and Rice [12]. When these engineered fields are explicitly designed to cancel out internal non-adiabatic contributions, the protocol is designated as counterdiabatic driving (CD) [13]. Initially formalized for closed quantum systems, the Berry-type formulation and quantum adiabatic theorem has already been adapted in the context of stochastic dynamics [14], while counterdiabatic principles in general have since been rigorously extended to classical stochastic thermodynamics architectures, providing a powerful toolkit for minimizing dissipation in fluctuating, finite-time regimes [13–18]. A common thread in prior work on classical stochastic STA is that the counterdiabatic field is constructed from either a series of generalized coordinate transforms (diffeomorphisms) [15, 19, 20] (see Appendix I), a probability-current continuity equation [17], or MBARbased estimators [21]. Owing to the computational cost of these approaches, machine-learning methods have become increasingly prominent, including flowmatching [22–24], (stochastic) normalizing flows (NFs) [25], learned switching protocols [26] and virtual escorted trajectories [8] all targeting efficient estimation of free energy differences. Controlled fluctuation theorems have similarly motivated neural samplers that derive variational objectives directly from the Jarzynski and Crooks identities and effectively optimise escorting fields via Radon–Nikodym derivative-based objectives. These include Controlled Monte Carlo Diffusion [27], the Non-

Equilibrium Transport Sampler [28], and discrete extensions [29]. In this work, we implement a spectral framework for the counterdiabatic driving (CD) [13], which we call Liouvillian counterdiabatic driving (LCD) for classical stochastic systems governed by the Fokker-Planck equation. The underlying Fokker-Planck generator (Liouville operator) L̂(t) when discretized over the spatial coordinates is intrinsically a non-symmetric rate matrix that admits an exact biorthogonal decomposition [30, 31] whose structure directly encodes all the information of the non-adiabatic lag accumulated during finite-time driving. Exploiting this intrinsic feature, we derive a closed-form spectral formula to construct the counterdiabatic correction L̂CD (t) using the Liouville operator itself. We further show that this spectral correction collapses, under detailed balance, to a rank-one operator determined in closed form by the equilibrium distribution; the biorthogonal machinery becomes equivalent to the known closed-form escorting for equilibrium problems, and its value is primarily structural—establishing the transitionless-driving analogy, exposing the spectralgap origin of the lag. Our LCD systematically eliminates this lag, escorting the probability distribution along the instantaneous Boltzmann equilibrium πeq (t) at arbitrary protocol speeds. Consequently, our deterministic verification shows that LCD enforces instantaneous equilibrium tracking ensuring zero dissipated work (Wdiss = 0) for arbitrary protocol speeds (not to be confused with free-energy estimation over sampled trajectories, done via the (E)JE [Eqs. (1, 2)] which is not performed in this work) . Previous methods predominantly operate without any explicit reference to the spectral features of L̂(t), our approach fills this gap by working directly within its time-dependent biorthogonal eigenbasis, based on a similar approach developed in the context of quantum mechanics [11]. This yields an exact, optimization-free construction on the discrete generator that requires no coordinate transformations or learned representations. To this end, we introduce our framework that has the following features: • Alternative to coordinate-based methods. Prior constructions of the escorting field rely on finding a series of bijective coordinate reparametrisations (diffeomorphisms) of phase space [15, 19]. Our counterdiabatic correction is obtained directly from the spectral decomposition of the Fokker–Planck generator itself, via a time-dependent gauge-like transformation of the probability density, without any reference to Jacobian determinants. • Formal correspondence to Berry/Demirplak–Rice methodology. The spectral correction L̂CD (t) mirrors the structural form of the quantum counterdiabatic Hamiltonian, adapted to the non-symmetric Fokker–Planck generator via biorthogonal decomposition. This enables the reduction (elimination)

3 of non-adiabatic lag and as a result excess dissipated work.

• Systematic spectral truncation and numerical scalability. The counterdiabatic correction L̂CD (t) naturally admits a controlled low-rank approximation by restricting the biorthogonal expansion to the M < N relaxation modes of the generator. Although it does not correspond to the perfect escorting, this offers a highly predictable numerical trade-off, where the number of spectral modes required to reach a target accuracy systematically scales with the driving protocol speed and more importantly the inverse spectral gap. This feature enables computationally efficient implementations in large-scale systems where tracking the full spectrum is unfeasible.

The remainder of this paper is organized as follows. Section II introduces the matrix representation of the Fokker-Planck equation obtained via a probabilityconserving finite-difference spatial discretization, establishing the resulting non-symmetric structure of the discrete Liouville operator. Section III reformulates the finite-time Fokker-Planck dynamics within an adiabatic (co-moving) frame via a time-dependent gauge-type transformation, explicitly decoupling the evolution into an instantaneous-equilibrium component and its nonadiabatic component. Section IV develops the biorthogonal decomposition of L̂(t) and derives the fundamental structural properties of the instantaneous eigenbasis {rn (t), ℓTn (t)}, showcasing how the right zero mode tracks the Boltzmann distribution r0 (t) = πeq (t) while the left zero mode satisfies ℓT0 = 1T as a direct consequence of probability conservation. Our main result—the classical analogue to Berry’s transitionless quantum driving [32] for stochastic systems governed by the Fokker– Planck equation is derived in Section V, yielding an exact formulation for the counterdiabatic correction L̂CD (t) constructed directly from the biorthogonal eigenpairs of L̂(t). Section VII presents numerical demonstrations across two prototypical model systems: a time-varying double-well (DW) potential and a harmonically confined Brownian particle in the overdamped setting. We assess tracking fidelity across varying protocol durations using the total variation distance, the Kullback-Leibler divergence (DKL ), and the trajectory-wise vanishing of the dissipated work Wdiss . In Section VII D, we investigate the numerical scalability of a spectrally truncated approximation to L̂CD (t) using a restricted subset of M < N relaxation modes, establishing a systematic trade-off for improving convergence. Lastly, in Section VII E, we identify potentials, where our spectral CD methodology has clear disadvantages, mainly tied to vanishing spectral gaps.

II.

MATRIX REPRESENTATION OF THE FOKKER-PLANCK EQUATION

Consider an overdamped Brownian particle in a timedependent potential V (x, ζ(t)), where ζ(t) is an externally controlled parameter. The probability density ρ(x, t) evolves according to the Fokker–Planck equation written in the Liouville form [33] (in 1D space): ∂ρ(x, t) ∂ ∂V (x, t) ∂ 2 ρ(x, t) , (4) = ( ) ρ(x, t) + β −1 ∂t ∂x ∂x2 } {z |∂x L(x, t) ρ(x, t)

∂ ∂V (x, t) ∂ where, L(x, t) := ( ·) + β −1 is the Li∂x ∂x ∂x ouville operator. The gradient of the time-varying po(x,ζt ) tential, corresponds to the drift velocity: ∂V ∂x = −v(x, t). Driven by a time-varying external potential V (x, ζt ) with a control parameter ζt , finite-time protocols inevitably induce a lag in the form of non-adiabatic (diabatic) excitations that prevent the time-evolved FokkerPlanck probability density from tracking a target distribution with high fidelity. This lag manifests thermodynamically as an excess dissipated work Wdiss = W − ∆F, whose magnitude scales with the protocol speed |ζ̇t |. In the quasi-static limit (ζ̇t → 0), the system evolves reversibly along a sequence of instantaneous equilibria without diabatic transitions, yielding a vanishing dissipation Wdiss → 0. As considered in other works, we have the same goal as in [7], to track the instantaneous equilibrium distribution e−βV (x,ζ(t)) := e−β V (x,ζt )−F (ζt ) , (5) πeq (x, t) = Z(t) at all times during the Fokker–Planck evolution for any arbitrary finite-time protocol speeds. The instantaneous equilibrium distribution defined in Eq. (5), is a stationary solution of the Liouville operator satisfies L(x, t) πeq (x, t) = 0 at every time-instant [7, 34].

Discretisation: Turning the PDE into a Matrix-Valued ODE

We discretise x on N equally spaced grid points −1 {xi }N i=0 with spacing ∆x. The continuous density ρ(x, t) can be written in terms of a column vector, with the spatial points discretized: T ρ(t) = ρ0 (t) · · · ρN −1 (t) ∈ RN , ρi (t) ≈ ρ(xi , t) ∆x, (6) P with normalisation i ρi = 1. The differential Liouville operator L(x, t) described in Eq. (4) becomes an N × N rate matrix L̂(t), and the scalar PDE (4) becomes ∂ρ(t) = L̂(t) ρ(t). ∂t

(7)

4 Structure of the rate matrix. The entries of L̂(t) are transition rates between neighbouring grid points, obtained from the Sasa–Tasaki discretisation [35] that describes local transition rates between adjacent lattice sites i and i + 1 in the form: Li+1,i (t) =

Li,i+1 (t) =

1 e−β∆Vi (t)/2 , β∆x2

(8a)

1 e+β∆Vi (t)/2 , β∆x2

(8b)

where ∆Vi = V (xi+1 , ζt ) − V (xi , ζt ), with reflecting boundary conditions. It is important to note that L̂(t) is not a symmetric matrix by construction, but can be made symmetric (see Section IV for details). Its diagonal enP tries enforce probability conservation: Lii = − j̸=i Lji , so that every column of L̂ sums to zero: N −1 X

Lij (t) = 0

for all j, t.

(9)

i=0

By the Perron–Frobenius theorem, L̂(t) has a unique zero eigenvalue, all other eigenvalues are strictly negative, and the zero-eigenvalue right eigenvector is exactly πeq (t). The rate matrix L̂(t) is a highly sparse matrix ∼ O(d N ), in our case amounting to only (2d + 1)N non-zero elements. Thus, for despite much finer discretizations, the rate matrix exhibits a sparse representation. We cover this in detail in the Appendix VII F.

III.

NON-ADIABATIC LAG VIA FRAME TRANSFORMATIONS

In order to quantify the non-adiabatic excitations induced by finite-time protocols, we transform the probability vector via a time-dependent gauge-like transformation, ρ̃(t) = D −1 (t) ρ(t).

(10)

Operationally, D −1 (t) projects ρ(t) ∈ RN onto a spectral domain, changing the basis from a natural spatial grid basis, where components denote local site probabilities to the instantaneous (time-dependent) eigenbasis of the system’s dynamics. Crucially, unlike external coordinateparametrizations used in normalizing flows, the underlying coordinate space remains untransformed; what varies is strictly the representation of the state vector itself, allowing us to cleanly isolate non-adiabatic transitions during the time-varying external potential V (x, ζt ). Thus, this transformation maps the system’s evolution directly into the adiabatic frame of reference [9, 11]. The explicit role of D(t) as the diagonalizing operator emerges naturally when examining the dynamics in the adiabatic frame. Differentiating this transformation [Eq. (10)] with

respect to time and substituting the Liouvillian evolution equation from Eq. (7) yields: ∂ ρ̃ ∂ρ ∂D −1 = D −1 +( )ρ ∂t ∂t ∂t ∂D −1 = D −1 LD ρ̃ + ( )D ρ̃ ∂t ∂D = Λ ρ̃ − D −1 ρ̃. ∂t

(11)

where the similarity transformation D−1 LD = Λ maps the Liouville operator to its diagonal spectrum of relaxation rates (see Section IV B for numerical corrobora−1 −1 tion). The last equality uses ∂(D∂t D) = 0 ⇒ ( ∂D ∂t )D = −D −1 ( ∂D ∂t ). The equation of motion in the adiabatic frame then reads ∂ ρ̃ = Λ(t) ρ̃ − Γ(t) ρ̃ , (12) | {z } | {z } ∂t adiabatic

non-adiabatic

where Λ(t) is diagonal matrix in this representation, while the latter term in Eq. (12) is identified as the adiabatic gauge connection (AGC) pertaining to the FokkerPlanck equation, a term often used in the STA community [9, 10]. ∂D(t) . (13) ∂t Physical interpretation. In the quasi-static limit where the parameter velocity |ζ̇t | vanishes (Γ → 0), the spectral modes completely decouple, reducing the adiabatic frame dynamics to ∂ ρ̃n /∂t = λn ρ̃n . Here, the zero mode (λ0 = 0) acts as an invariant manifold while all higher-order excited modes relax to zero. Under finitetime driving (Γ ̸= 0), however, the off-diagonal elements encapsulated within the AGC dynamically couple the zero mode to these transient relaxation modes. This spectral mixing constitutes the exact mathematical origin of the non-adiabatic lag: because the basis is timedependent, the probability density cannot instantly relax to the moving equilibrium profile. Instead, the driving continuously scatters probability density into higher relaxation modes. This is formalized in Eq. (11) where a clean orthogonal split of the stochastic evolution into an instantaneous equilibrium (ground-state) contribution and non-adiabatic transitions out of the instantaneous equilibrium manifold occurs. This geometric deviation directly generates the irreversible dissipated work Wdiss along a trajectory and must be exactly canceled or reduced at each time instant by the engineered counterdiabatic field. Γ(t) ≡ D −1 (t)

IV.

BIORTHOGONAL DECOMPOSITION OF THE LIOUVILLE OPERATOR

The Liouville operator L̂(t), written in rate-matrix form as in Eq. (8), is a real-valued matrix that is generally non-symmetric. The non-symmetry of L̂ means

5 that left and right eigenvectors are distinct objects that must be treated separately. We define the right and left eigenproblems as L̂(t) rn (t) = λn (t) rn (t),

(14)

ℓTn (t) L̂(t) = λn (t) ℓTn (t),

(15)

where rn (t) ∈ RN is a column vector and ℓTn (t) ∈ R1×N is a row vector. The eigenvalues {λn } are real and nonpositive, 0 = λ0 > λ1 ≥ λ2 ≥ · · · ≥ λN −1 ,

(16)

a consequence of the Perron–Frobenius theorem applied to the rate matrix L̂(t). As the columns of D(t) are the right eigenvectors rn (t), so the transformation ρ̃(t) = D −1 (t)ρ(t) re-expands ρ from the natural grid basis {ei } (components = probability at site i) into the instantaneous eigenbasis {rn (t)} of L̂(t) (components = amplitude of relaxation mode n). Biorthogonality. Because L̂(t) is non-symmetric due to the spatial discretization, its left and right eigenvectors form a biorthogonal system satisfying ℓTn rm =

N X (ℓn )i (rm )i = δnm .

(17)

i=1

Collecting the right eigenvectors as columns of the matrix D = [r0 r1 · · · rN −1 ], and the left eigenvectors as rows of D −1 , the biorthogonality relation (17) is equivalent to D −1 D = I. The associated completeness relation (resolution of the identity) reads I =

N −1 X

rn ℓTn

⇐⇒

−1

= I,

(18)

n=0

so that any probability vector ρ(t) ∈ RN admits the expansion ρ(t) =

N −1 X

ρ̃n (t) rn (t),

ρ̃n (t) = ℓTn (t) ρ(t),

Physically, ℓT0 = 1T expresses probability conservation: since 1T ρ = 1 for all normalised ρ, and 1T ∂ρ/∂t = 1T L̂ρ = 0T ρ = 0, the total probability is conserved exactly. Moreover, it is important to emphasize that the detailed-balance condition holds throughout this paper; the non-symmetric nature of the transition rate matrix L̂(t) stems entirely from the spatial discretization of the continuous system. Consequently, the dynamics can be cast into a symmetric gauge via the similarity transformation

(19)

−1/2 1/2 L̂symm (t) = πeq (t) L̂(t) πeq (t).

(22)

This transformation reduces the complexity of the spectral problem from a biorthogonal framework to a standard orthogonal eigendecomposition, albeit evaluated with respect to a weighted inner product under the met−1 ric πeq (t). This symmetrized version is fully equivalent to our non-symmetric version. Under such a similarity transformation, one can view it analogous to the Hamiltonian [11, 12] in the quantum setting, which is an Hermitian operator, unless dealing with non-Hermitian Lindbladians in open quantum systems [37, 38]. Because we choose to work directly with the raw, spatially discretized rate matrix rather than its symmetrized counterpart, we need to perform a biorthogonal decomposition where the left and right zero modes do not coincide, i.e., T (t). This structural asymmetry highlights ℓT0 = 1T ̸= πeq the probability-preserving nature of the unskewed classical representation. We later cover the computational advantages of working in a symmetric representation of the rate matrix [Eq. (22)] in Section VII E. Probability conservation condition in the adiabatic frame. Expanding the laboratory frame Fokker– Planck distribution ρ(t) in terms of the biorthogonal eigenbasis via the resolution of identity described in Eq. (18): ρ(t) =

N −1 X

ρ̃n (t) rn (t),

ρ̃n (t) ≡ ℓTn (t) ρ(t).

(23)

n=0

where, this corresponds to a change of basis in the space of probability vectors, not a change of coordinates in physical space. The zero mode and its left partner. The zero eigenvalue λ0 = 0 carries a distinguished physical meaning. Its right eigenvector is the instantaneous equilibrium distribution [36]: r0 (t) = πeq (t),

L̂(t) r0 (t) = 0.

(20)

The corresponding left eigenvector satisfies ℓT0 L̂ = 0T . P

Since L̂ is a rate matrix with zero column sums, i Lij = 0 for all j, the all-ones row vector 1T = (1, 1, . . . , 1) is annihilated from the left: ℓT0 ≡ 1T ,

1T L̂(t) = 0T .

(21)

Applying the all-ones vector 1T = (1, · · · , 1) to both sides and using 1T rn = δn0 ([Eq. (21)] and biorthogonality): T

1 ρ(t) =

N −1 X n=0

ρ̃n (t) 1 rn =

N −1 X

ρ̃n (t) δn0 = ρ̃0 (t). (24)

n=0

Since ρ(t) is a probability vector, 1T ρ(t) = at all times, and therefore: ρ̃0 (t) = ℓT0 ρ(t) = 1T ρ(t) = 1,

i ρi (t) = 1

∀ t ∈ [0, τ ].

(25)

This is not a dynamical equation but an identity expressing probability conservation: it holds at every instant under any dynamics, regardless of whether ρ(t) tracks πeq (t). The zeroth component ρ̃0 is pinned always to unity by probability conservation, while the components

6 ρ̃n (t) for n ≥ 1 are are mode amplitudes (not probabilities) that can take any real value and unconstrained in this frame. In particular, N −1 X n=0

ρ̃n (t) = 1 +

N −1 X

ρ̃n (t) ̸= 1.

(26)

n=1

Thus, as derived in Eq. (26),in the adiabatic frame, the P condition ρ̃ = 1 need not hold in general. Infact, i i the continuous version of Eq.(26) and deviations from equilibrium distribution under time-dependent protocols using biorthonormal scheme was already shown in the context of the adiabatic theorem adapted to stochastic systems [14].

This is the analogue of the transitionless quantum driving approaches of [32, 39] adapted to the Fokker–Planck evolved instantaneous equilibrium distribution. The Eq. (29) has also appeared in the context of flowfields [17] and counterdiabatic driving currents on discrete graphs [36]. Moreover, one can cast Eq. (29) in the time-ordered exponential: Z t πeq (t) = T exp

N −1 X ∂ ρ̃n (t) = λn (t) ρ̃n (t) − Γnm (t) ρ̃m (t), ∂t m=0

(27)

Spectral solution via biorthogonal decomposition. To solve (29), we differentiate the stationarity condition L̂(t) πeq (t) = 0 [7, 34] with respect to t:

CD driving of the instantaneous equilibrium distribution

We seek an operator L̂CD (t) such that the modified dynamics i ∂ρ(t) h = L̂(t) + L̂CD (t) ρ(t), ∂t

(28)

admits πeq (t) as an exact solution for all t. Substituting ρ = πeq (t) into Eq. (28) and using the fact that the equilibrium distribution is always a stationary solution of the dynamics L̂ πeq = 0, gives the counterdiabatic condition to follow the instantaneous equilibrium distribution at all times: ∂πeq (t) = L̂CD (t) πeq (t). ∂t

L̂(t)

∂πeq (t) ∂ L̂(t) =− πeq (t). ∂t ∂t

(31)

i P ∂πeq ∂πeq = 0 (normalisation is preserved), i ∂t ∂t lies in the range of L̂ restricted to the non-zero modes. Projecting (31) onto the biorthogonal basis via the resolution of identity (18) and using ℓTn L̂ = λn ℓTn , we obtain

Since

X ℓT (t) (∂ L̂/∂t) r0 (t) ∂πeq (t) n =− rn (t). ∂t λn (t)

(32)

n̸=0

where Γnm = [D −1 (∂D/∂t)]nm = ℓTn (t)(∂rm (t)/∂t). For n = 0, i.e. corresponding to the instantaneous equilibrium distribution: since the eigenvalue is λ0 = 0 the corresponding matrix elements of the non-adiabatic contributions are Γ0m = ℓT0 (t)∂rm (t)/∂t = 1T ∂rm /∂t = ∂(1T rm )/∂t = 0 for all m (because 1T rn = δn0 is timeρ̃0 independent), Eq. (27) gives ∂∂t = 0, consistent with ρ̃0 = 1 for all t. For n ≥ 1: the first term drives ρ̃n → 0 at rate |λn | (relaxation), while Γnm sources ρ̃n from mode m (non-adiabatic coupling).

(30)

V. ENGINEERING THE PERFECT CD TERM FOR THE INSTANTANEOUS EQUILIBRIUM DISTRIBUTION

The dynamics of ρ̃n (t) as derived in Eq. (12) can be obtained by substituting the expansion [Eq. (23)] into the adiabatic-frame equation [Eq. (12)] and projecting onto mode n by applying ℓn (t) from the left:

L̂CD (s) ds πeq (0).

(29)

The n = 0 term vanishes because 1T (∂πeq /∂t) = 0. The Eq. (32) can also be obtained from the component form derived in Eq. (27). The matrix element (n, 0) of Γ gives the rate at which the zero mode (instantaneous equilibrium distribution) leaks into other modes n: ∂D ∂r0 ℓT (∂ L̂/∂t) r0 ]n0 = ℓTn =− n , ∂t ∂t λn (33) where the last step follows from differentiating L̂r0 = 0. The non-adiabatic leakage is therefore proportional to the matrix element ℓTn (∂ L̂/∂t) r0 and inversely proportional to the spectral gap |λn |. The spectrum {λn (t)} of the Fokker-Planck generator L̂(t) characterizes an ensemble of relaxation rates. These characterizing the fluctuations about the instantaneous equilibrium, 0 = λ0 > λ1 ≥ λ2 ≥ · · · ≥ λN −1 , carrying units of inverse time (s−1 ). Accordingly, each instantaneous eigenmode rn (t) decays exponentially as eλn t , establishing |λn |−1 as an intrinsic relaxation timescale on which the system clears fluctuations away from the instantaneous equilibrium distribution. While this relaxation spectrum becomes continuous in the thermodynamic continuum limit (N → ∞)—structurally akin to the continuous Laplacian −∂x2 on R—it is discretized here by our finite spatial grid (Nx ). Crucially, the slowest non-zero mode λ1 (t) defines the spectral gap of the generator, which governs the dominant, macroscopically observable relaxation timescale: τrelax (t) = |λ1 (t)|−1 . Γn0 (t) = [D −1

7 One can now construct the Liouvillian counterdiabatic operator L̂CD required for adiabatic tracking of Eq. (29) as a sum of rank-1 operators that, when acting on πeq , reproduces (32). Using ℓTn πeq = δn0 (biorthogonality) and ℓT0 πeq = 1T πeq = 1, the minimal-rank solution, in the resolvent form: X ℓT (∂ L̂/∂t) r0 n L̂CD (t) = − rn ℓT0 . (34) λn (t) n̸=0

Each term in (34) is a rank-1 matrix rn ℓT0 weighted by the matrix element ℓTn (∂L/∂t)r0 /λn (see derivation in Section III). This corrective field is a biorthogonal adiabatic gauge connection matrix of our nonsymmetric Fokker–Planck generator, in exact formal analogy with the transitionless quantum driving methodologies of Berry [32], Demirplak and Rice [39]. This general form admits a substantial simplification in equilibrium, which we derive next. As shown in Eq. (35), 1T L̂CD = [1T (∂πeq /∂t)]1T = [∂(1T πeq )]1T = 0. (or equivalently, it follows from 1T rn = 0, n ̸= 0, by biorthogonality (ℓT0 rn = δ0n = 0)), the corrected generator L̂eff = L̂ + L̂CD conserves probability exactly. We emphasize, however, that L̂eff is not a stochastic rate matrix: L̂CD = (∂πeq /∂t)1T is dense with sign-indefinite off-diagonal entries, so L̂eff generally has negative offdiagonals. It is a probability-conserving linear operator that renders πeq (t) an instantaneous solution, not the generator of a Markov jump process. While Iram et al. [14] has derived the mathematical deviation from the instantaneous equilibrium distribution via the Fokker-Planck adiabatic driving, their lag removal concerns which physically accessible control realizes the counterdiabatic generator (involving no diagonalization) in an evolutionary/biological systems setting, whereas we treat its numerical construction and stability on a discretized Fokker–Planck operator. What Eq. (34) adds is an explicit construction of that generator from the spectral data of L̂(t), rather than an implicit definition or coordinate transforms. This is the step that makes the correction computable without an analytic target, exposes the inverse spectral gaps λ−1 n that control its conditioning, and admits the systematic low-rank truncation for escorting the Fokker–Planck equation and removing dissipated work. B.

The counterdiabatic operator is rank one and has a closed form

The spectral correction Eq. (34) admits a striking simplification for the instantaneous equilibrium problem. Even though, it is written as a sum over all N − 1 relaxation modes, every term carries the same left factor ℓT0 = 1T . Factoring it out and recognising the bracket as the spectral representation of ∂πeq (t)/∂t [Eq. (32)], the sum collapses to a single rank-one operator, L̂CD (t) = ∂πeq /∂t 1T , (35)

which satisfies the counterdiabatic condition [Eq. (29)] exactly, L̂CD πeq = (∂πeq /∂t)(1T πeq ) = ∂πeq /∂t, since 1T πeq = 1. Using πeq ∝ e−βV (x,ζt ) , the source term is available analytically, ∂πeq ∂V ∂V = −β ζ̇t πeq ⊙ −⟨ ⟩π , ∂t ∂ζ ∂ζ eq

∂V dF ⟩π = , ∂ζ eq dζ (36) where ⊙ is the Hadamard product over grid points and V(ζt ) ∈ RN is the column vector representation of the external potential evaluated at each spatial grid point. For the equilibrium (detailed-balance) problem the exact counterdiabatic generator therefore can be computed without involving eigenvalues nor the eigenvectors of L̂(t): it is fixed at O(N ) cost by the Boltzmann weight and its parametric derivative alone. Eq. (36) is the matrix form of the classical counterdiabatic term obtained in the continuum by Guéry-Odelin et al. [40], ⟨

∂πeq (x, ζt ) ∂W ∂W = −β −⟨ ⟩ πeq (x, ζt ). ∂t ∂t ∂t Eqs. (34) and (35) are the same operator, not two competing constructions. Differentiating the stationarity condition L̂πeq = 0 yields the singular linear system L̂ (∂πeq /∂t) = −(∂ L̂/∂t) πeq [Eq. (31)], solvable because 1T (∂ L̂/∂t)πeq = 0 and with a unique solution once the zero-mass condition 1T (∂πeq /∂t) = 0 is imposed. The spectral expansion (32) solves it through the pseudoP −1 T −1 inverse L̂† = n̸=0 λn rn ℓn , so the inverse rates λn are the eigen-weights of that pseudo-inverse rather than a property of L̂CD . The spectral expansion and the closed form (35) coincide mathematically; numerically, however, they behave differently. The closed form follows from direct parametric differentiation of the Boltzmann distribution πeq ∝ e−βV , using only elementary, well-scaled operations on πeq —no eigendecomposition, no operator inversion, and no division by vanishing scalars—and is therefore stable at any gap. The spectral representation instead exposes the inverse relaxation rates λ−1 n , the eigenvalues of the Liouvillian pseudo-inverse, which make its assembly severely ill-conditioned as the gap closes. Crucially, the closed form removes the inverse gaps only from the construction of L̂CD ; the physical spectral gap persists in the dynamics, where it continues to set the stiffness of the propagation and the relaxation rate of any numerical deviation from πeq . Thus, the spectral approach enables the following (i) makes the transitionless-driving analogy manifest (see also [14] for quantum adiabatic theorem adapted to Fokker–Planck equation), (ii) exposes the spectral-gap origin of the non-adiabatic lag induced during finite-time protocols – Eq. (34) (detailed in Section VII E) (iii) allows systematic truncation for engineering imperfect escorting (see Section VII D).

10 1

0.00

0.05

0.10

0 2 1

4 40

(a) Time varying double-well potential.

1 0.00

0.05

0 0.10

˙0

ζ̇

Potential profile, V(x, 0 (t))

Potential profile, V(x, ζ(t))

5 0 3

(b) Time varying harmonic potential.

FIG. 1: Time-varying experimental potentials and control protocol profiles. (a) Double-well potential Vdw (x, ζt ) = x4 − 2x2 + ζt x with tilt parameter ζt driven from ζi = −1 to ζf = +1 over t ∈ [0, 0.1]. The potentials start at t = 0.0 (violet coded) and modulates to t = 0.1 (orange coded) at the end of the protocol. Inset: protocol ζ(t) and driving rate ζ̇(t). (b) Harmonic trap Vh (x, κ0 (t)) = 12 κ0 (t)x2 with stiffness κ0 (t) ≡ ζt driven from κi = 1 to κf = 4 over t ∈ [0, 0.1]. Inset: protocol ζ(t) [κ0 (t)] and driving rate ζ̇(t) [κ̇0 (t)]. Both protocols follow the smoothstep form [Eq. (S29)], ensuring zero first and second derivatives at t = 0 and t = τ .

VI.

NUMERICAL SETUP

The system [Eq. (28)] is stiff, with fast modes at rates proportional to |λN −1 | for a protocol acting on a timescale τ . Explicit integrators (Euler, RK4) would require ∆t < |λN −1 |−1 , making long simulations prohibitively expensive. We instead evolve ρ(t) via the matrix exponential [41], ρk+1 = exp Ω(tk , tk + ∆t) ρk , (37) where Ω is the fourth-order Magnus exponent (Appendix V C), evaluated at two Gauss–Legendre nodes √ t1,2 = tk + (1/2 ∓ 3/6) ∆t. This retains the unconditional stability of the matrix exponential while achieving O(∆t4 ) convergence. Both potentials are driven by the smooth protocol ζ(t) = ζi + (ζf − ζi ) s3 (6s2 − 15s + 10),

s = t/τ, (38)

whose first and second derivatives vanish at t = 0 and t = τ . We consider two model systems: a double-well (DW) potential as shown in Fig. 1a, with Vdw (x, ζt ) = x4 − 2x2 + ζt x with tilt parameter driven from ζi = −1 to ζf = +1, discretised on N = 80 grid points over x ∈ [−2.5, 2.5] (∆x ≈ 0.063); and a harmonic trap Vh (x, κ0 ) = 12 κ0 (t) x2 with stiffness κ0 (t) ≡ ζt driven from κi = 1 to κf = 4, discretized on N = 80 grid points over x ∈ [−4.0, 4.0] (∆x ≈ 0.101). Both systems are evolved at inverse temperature β = 1. As already described in Section II A, we employ the

Sasa–Tasaki discretisation [35, 42] to construct the rate matrix (for both the bare and CD-evolved rate matrix). To minimize temporal discretization errors across the wide range of protocol durations τ ∈ [10−3 , 10], the integration time-step is dynamically scaled between ∆t = 10−7 and ∆t = 10−4 . By choosing smaller timesteps for shorter protocol durations, we ensure high temporal resolution in the strongly driven regime, where the parameter velocity |ζ̇t | is large, thereby maintaining a uniformly high numerical tracking fidelity across all driving speeds. Additionally, in the numerics, the P right zero eigenvector r0 (t) is normalized to unit sum, i (r0 )i = 1, uniquely identifying it with the instantaneous Boltzmann distribution πeq (xi , t) ∝ e−βV (xi ,ζt ) .

Computation of the dissipated work during evolution

The thermodynamic work performed on the system during a finite-time protocol ζt is defined in our spatially discretized setup as the time-integral of the instantaneous power: Z t W(t) = 0

∂V(ζs ) ∂s

Z tX N

T ρ(s)ds,

∂V (xi , ζs ) ρi (s)ds, ∂s 0 i=1

(39)

Probability density

t = 0.00 τ

t = 0.25 τ

t = 0.50 τ

t = 1.00 τ

t = 0.75 τ

max diff: 7.1e-16 max diff: 7.6e-16 max diff: 7.7e-15 max diff: 2.2e-15 max diff: 6.3e-15

0.10

r0 (t) πeq (t)

0.05 0.00

2.5

0.0

2.5

0.0

2.5

0.0

2.5

0.0

2.5

0.0

2.5

(a) Double-well potential zero-mode verification.

Probability density

t = 0.00 τ

t = 0.25 τ

t = 0.50 τ

t = 0.75 τ

t = 1.00 τ

2.5 0.0 2.5

max diff: 4.5e-15 max diff: 1.6e-15 max diff: 2.5e-15 max diff: 2.8e-15 max diff: 2.4e-15 0.05 0.00

r0 (t) πeq (t)

2.5 0.0 2.5

(b) Harmonic trap potential zero-mode verification.

FIG. 2: Spectral verification of the discrete stationary zero-mode. The numerically diagonalized right eigenvector r0 (t) (dashed grey profile) corresponding to the zero eigenvalue is coincides upto machine precision with the analytical instantaneous Boltzmann distribution πeq (t) (dotted green profile).

where, Vi (ζt ) = V (xi , ζt ). The average is taken over the time-evolved distribution ρ(t′ ) at each instant. Depending on the context, ρ(t′ ) corresponds to either the bare distribution [Eq. (7)] or the counterdiabatic quantity evolving under the escorted generator [Eq. (28)]. In the latter case, W(t) represents the total escorted work. The excess dissipated work is subsequently quantified as Wdiss (t) = W(t) − ∆F (ζt ),

(40)

where ∆F(ζt ) = F (ζt ) − F (ζ0 ). For a perfect counterdiabatic driving, ρ(t′ ) = πeq (t′ ), i.e., the excess dissipation vanishes identically throughout the protocol, since in ∂V this case ⟨ ∂V ∂t ⟩ρ = ⟨ ∂t ⟩πeq = dF/dt at each instant. For numerically integrating Eq. (39), we use the Simpson’s method, which is accurate to O(∆t4 ) as detailed in Section V E and consistent with the fourth-order Magnusbased evolution (see Section V C). VII.

MAIN EXPERIMENTAL RESULTS

To demonstrate the efficacy of the proposed Liouvillian counterdiabatic driving [Eq.(28)] framework, we study an overdamped Brownian particle in one dimension subjected to the two aforementioned driving fields: a time-varying double-well potential and

a modulated harmonic potential (see Section VI and Fig. 1). For both landscapes, the control parameter is swept across a wide range of protocol durations τ ∈ {0.001, 0.01, 0.1, 0.25, 0.5, 0.75, 1.0, 2.0, 4.0, 6.0, 10.0}, allowing us to evaluate the methodology from the highly non-equilibrium, fast-switching regime (τ = 10−3 ) to the near-adiabatic limit (τ = 10.0). We systematically quantify the performance of our spectral framework against three rigorous operational criteria: (i) Zero-mode tracking fidelity—verifying that the right zero-mode eigenvector r0 (t) coincides with the instantaneous equilibrium distribution πeq (t) up to machine precision at all times; (ii) Elimination of non-adiabatic lag—evaluating the ability of the driven dynamics to precisely track πeq (t) at each time instant, irrespective of the protocol speed |ζ̇t |; (iii) and Suppression of dissipation—confirming the vanishing of the excess dissipated work at every instant during the entire Fokker–Planck evolution, Wdiss (ζt ) → 0 as achieved by the flow field in EJE (see [Eqs. (2, 3)] bounded only by numerical discretization errors of the matrix-valued evolution. Additionally, we systematically investigate a spectrally truncated approximation of the counterdiabatic correction in Eq. (34) and also identify specific class of potentials where our spectral methodology breaks down.

t = 0.00 τ

t = 0.10 τ

t = 0.25 τ

0.075

0.04

0.050

0.02

0.025

0.00

t = 0.50 τ

0.06

ρ(x, t)

0.06

t = 1.00 τ

t = 0.75 τ

0.04 0.02 0.00

t = 0.00 τ

t = 0.10 τ

t = 0.25 τ

t = 0.50 τ

t = 0.75 τ

t = 1.00 τ

0.000 0.075 0.050 0.025

Bare

LCD

πeq (t)

(a) FP distribution evolution under DW potential.

0.000

2.5

0.0

2.5

Bare

2.5

0.0

LCD

2.5

0.0

2.5

πeq (t)

(b) FP distribution evolution under harmonic potential.

FIG. 3: Evolutionary snapshots of the probability density distributions. Spatiotemporal density profiles ρ(x, t) comparing the bare driving dynamics (solid orange line) and our spectral-based Liouvillian CD driving (LCD) (grey) against the instantaneous target Boltzmann distribution πeq (x, t) (dashed green) across selected temporal fractions of the protocol for (a) the double-well landscape and (b) the harmonic confinement trap.

Zero-mode alignment with instantaneous equilibrium

As established in Section IV, a fundamental requirement of the matrix-valued Liouvillian framework is that the right zero-mode eigenvector identically mirrors the instantaneous equilibrium distribution. This condition is robustly corroborated by our numerical experiments: r0 (t) coincides with the instantaneous equilibrium distribution πeq (t) at all times, irrespective of the driving rate or potential landscape. To verify the algebraic definition in Eq. (20), we evaluate this spectral baseline across uniform temporal snapshots, t/τ ∈ {0, 0.25, 0.5, 0.75, 1.0}, for both the double-well – Figure 2a and harmonic – Figure 2b landscapes. In all cases, the states exhibit absolute structural convergence down to machine precision. The global coordinate-wise deviation remains strictly bounded by maxt ∥r0 (t) − πeq (t)∥∞ = O(10−15 ), confirming that the discrete zero-mode manifold (λ0 = 0) precisely maps to the instantaneous equilibrium target under arbitrary driving profiles. B.

Metrics and Tracking Fidelity

Propagating the Fokker-Planck density under the effective counterdiabatic (CD) Liouvillian L̂eff (t) = L̂(t) + L̂CD (t) [Eq. (28)] is engineered to eliminate the nonadiabatic lag induced by finite-time driving. The spectral P ℓT (∂ L̂(t)/∂ζ)r correction term L̂CD (t) = −ζ˙t n̸=0 n λn (t) 0 rn ℓT0 , perfectly escorts the state vector along the instantaneous equilibrium manifold. To quantify the fidelity of this tracking shortcut, we evaluate the total variation distance (TVD), TVD ρ(t), πeq (t) = 12 ∥ρ(t) − πeq (t)∥1 , (41)

and the Kullback-Leibler (KL) divergence, DKL ρ(t)∥πeq (t) , between the time-evolved density and the instantaneous target. Double-well potential.—The tracking metrics across different protocol durations are compiled in Table I. Compared to the uncontrolled bare evolution [Eq. (7)], the CD framework yields an overwhelming improvement: the maximum TVD remains strictly bounded within O(10−12 )–O(10−11 ), while the corresponding KL divergence sits at the machine-precision floor of O(10−16 ). The time-resolved profiles of these metrics are shown in Figs. 4a and 4b for a representative protocol (τ = 0.1), demonstrating continuous, uniform lag suppression. Harmonic potential.—As shown in Table I, the numerical CD framework exhibits equally robust performance under a driven harmonic potential. To establish a rigorous benchmark, we compare our numerical results against the exact analytical CD driving scheme derived in Section V D. The numerical framework recovers the analytical target up to 10–12 orders of magnitude over the bare dynamics. The minor discrepancy between the numerical and exact analytical tracking metrics represents minor time-stepping and discretization artifacts arising from propagating the full matrix-valued generator rather than scalar-valued evolution used for the analytic CD simulations. Crucially, the KL divergence remains pinned at O(10−16 ), matching the exact analytical solution to within numerical noise. This near-perfect overlap and the systematic elimination of the non-adiabatic lag are visually confirmed by the temporal evolution of the metrics for τ = 0.1 in Figs. 4c and 4d, as well as the spatial snapshots provided in Fig. 3b. We note that, KL is dominated by the far tails where πeq is exponentially small, so it saturates near machine precision, whereas TVD reflects the bulk error set by the

11 TABLE I: Maximum KL divergence and maximum TVD for bare FP and LCD driven FP at several protocol durations τ for both the double well and harmonic trap potentials. Double Well Protocol duration (τ )

bare DKL

LCD DKL

10−3 (∆t = 10−7 ) 0.01 (∆t = 10−6 ) 0.10 (∆t = 10−5 ) 0.25 (∆t = 10−5 ) 0.50 (∆t = 10−5 ) 0.75 (∆t = 10−5 ) 1.00 (∆t = 10−5 ) 2.00 (∆t = 10−5 ) 4.00 (∆t = 10−4 ) 6.00 (∆t = 10−4 ) 10.0 (∆t = 10−4 )

1.40 1.39 1.27 1.13 0.97 0.85 0.75 0.48 0.25 0.15 0.07

4.68 × 10−16 3.63 × 10−16 3.88 × 10−16 3.98 × 10−16 4.13 × 10−16 4.26 × 10−16 4.25 × 10−16 4.76 × 10−16 4.19 × 10−16 4.36 × 10−16 4.27 × 10−16

Harmonic Trap

bare

TVD

0.68 0.67 0.66 0.64 0.60 0.57 0.53 0.44 0.32 0.25 0.18

TVD

3.69 × 10−12 1.71 × 10−12 2.14 × 10−12 1.20 × 10−12 1.08 × 10−12 9.74 × 10−13 8.76 × 10−13 1.01 × 10−12 1.28 × 10−12 7.14 × 10−13 5.97 × 10−13

quadrature/construction. C.

LCD

Trajectory-Resolved Dissipated Work

Beyond the terminal values, a complete assessment of the escort requires the ensemble-averaged dissipated work resolved continuously along the driving path, Wdiss (t) = W(t) − ∆F(t). The thermodynamic and distributional pictures are tied by a single identity: because ⟨∂V /∂t⟩πeq = dF/dt at every instant [Section VI A], perfect tracking of the instantaneous equilibrium, ρ(t) = πeq (t), is equivalent to vanishing dissipation, Wdiss (t) = 0. The trajectory-resolved Wdiss (t) is therefore the thermodynamic corollary of the distributional tracking established in Section VI A; reporting it additionally confirms that the escorting equality holds for the work functional. While the LCD isolates the system from non-adiabatic excitations, ensuring that the dissipated work vanishes, this tracking is bought at the expense of a transient energetic cost. The counterdibatic L̂CD (t) must perform work along the protocol path to suppress mode excitations, a thermodynamic cost of this shortcut that scales aggressively as the protocol duration τ approaches the limits of the system’s spectral gap (see for instance [43]). Computing this excess dissipated work during the Fokker–Planck evolution is detailed in Section VI A. Double-well potential.—The dissipation profile for a rapid switch (τ = 0.1) is shown in Fig. 5a. In the bare evolution the accumulated dissipated work climbs steadily as the landscape is modulated; under the LCD it is suppressed continuously. As reported in Table II, the maximal transient deviation maxt |Wdiss (t)| remains at the O(10−12 ) floor set by the fourth-order work quadrature and Magnus evolution [Section VI A]—roughly four orders of magnitude above the double-precision floor— confirming that the lag is canceled continuously rather than only at the protocol endpoints. Harmonic potential.—The cancellation is further

bare DKL

LCD DKL

TVDbare

TVDLCD

0.80 0.76 0.52 0.32 0.19 0.12 0.09 0.03 0.01 0.006 0.002

2.25 × 10−16 4.93 × 10−16 4.26 × 10−16 4.33 × 10−16 4.56 × 10−16 4.28 × 10−16 4.60 × 10−16 4.76 × 10−16 4.21 × 10−16 4.70 × 10−16 4.85 × 10−16

0.32 0.32 0.27 0.23 0.18 0.15 0.13 0.09 0.05 0.04 0.02

2.3 × 10−12 2.2 × 10−12 1.95 × 10−12 1.61 × 10−12 1.30 × 10−12 1.01 × 10−12 9.34 × 10−13 5.88 × 10−13 4.05 × 10−13 2.61 × 10−13 1.78 × 10−13

benchmarked against the exact analytic CD solution for the driven harmonic trap (Section V D). Figure 5b shows |Wdiss (t)| at τ = 0.1 for the bare, analytic-CD, and spectral-LCD protocols; broader metrics are collected in Table II. The bare dynamics accumulates substantial dissipation, whereas both the analytic and spectral CD drives maintain |Wdiss (t)| ≲ 10−12 at every instant, so the spectral construction reproduces the continuum shortcut to within the time-discretisation accuracy of the Simpson integral. The equilibrium free-energy difference for this system is ∆F = 12 ln(κf /κi ) = ln 2 ≈ 0.6931, independent of the protocol duration τ (derivation in Section VI A). Figure 7 shows the total work W(τ ) across protocol durations. Under LCD, W LCD (τ ) = ∆F to within ∼ 10−12 for all τ ∈ [10−3 , 10] (Table S3), whereas the uncontrolled work ranges from W ≈ 1.50 at τ = 10−3 to W ≈ 0.73 at τ = 10, approaching ∆F only in the quasi-static limit. This equality is the expected consequence of exact tracking: with ρ(t) = πeq (t), the work Rτ reduces to W(τ ) = 0 ⟨∂V /∂t⟩πeq dt = ∆F. We stress a methodological distinction from trajectory-sampling approaches: rather than estimating ∆F from a Jarzynski average ⟨e−βW ⟩ over stochastic realizations, we propagate the full density ρ(t) deterministically, so W(τ ) is the exact ensemble-averaged work for each τ —a single number with no statistical error, making the comparison with ∆F exact rather than asymptotic.

Spectral Truncatation Performance

Our Liouvillian CD correction term Eq. (34) admits a truncation, where one can restrict to a first set of spectral modes M < N (or even M ≪ N ). We detail this in Section VII D. Although, the efficacy of such a truncation is fundamentally governed by two quantities in Eq. (34) (i) the speed of the driving protocol |ζ̇t |, which affects the relaxation to other higher modes and (ii) the instantaneous

10 2

10 5

Bare LCD

10 8 10 11 10 14

Bare LCD

10 8

10 11 10 14

0.00

0.25

0.50

t/τ

0.75

0.00

1.00

(a) TVD against normalized time evolved under DW potential.

10 1

TVD(ρ, πeq )

D KL (ρkπeq )

10 2

0.25

0.50

t/τ

0.75

1.00

(b) KL divergence against normalized time evolved under DW potential.

10 2

10 4

Bare Analytic CD LCD (spectral)

10 7 10 10 10 13

10 5

D KL (ρkπeq )

TVD(ρ, πeq )

Bare Analytic CD LCD (spectral)

10 8

10 11 10 14

10 16 0.00

0.25

0.50

t/τ

0.75

1.00

0.00

0.25

0.50

t/τ

0.75

1.00

(d) KL divergence against normalized time evolved under an harmonic trap potential.

FIG. 4: Informational metrics comparison across potentials. Total variational distance (TVD) and KL divergence between the time-evolved FP density and the instantaneous equilibrium distribution. Top row (a, b) shows results of the TVD and KL when evolved under a DW potential respectively, whereas the bottom row (c, d) details the same for the harmonic potential. Bare (uncontrolled) evolution (solid orange line), LCD (dashed grey line) and if available, Analytic CD (solid green line). −1 inverse relaxation gaps ∆−1 = n (t) ≡ |λ0 (t) − λn (t)| −1 |λn (t)| , which dictate a closing spectral gap and the presence of nearly-degenerate modes. Although for simpler potentials such as DW or harmonic traps, omitting highly excited modes yields an incomplete cancelation of the non-adiabatic lag (imperfect escorting), as shown in Figures 6a [DW] and 6b [harmonic] and Table S6 for τ = 0.1. As already noted in Section V B, the full counterdiabatic operator is rank one, the spectral decomposition of that rank-one correction into individual mode

contributions identifies which spectral modes dominate the non-adiabatic lag and how many must be resolved for a given level of imperfect escorting. This information is inaccessible in the closed-form expression but essential for understanding the physics of finite-time driving. We show this for varied protocol speeds in Tables (S7, S8), it circumvents the computationally prohibitive requirement of full matrix diagonalization while achieving modest improvements in tracking fidelity. This enables a favorable numerical trade-off in terms of com-

10 1 10 4 10 7 10 10 10 13 10 16

10 2 10 5 10 8 10 11 10 14 10 17

0.00

0.25

0.50

t/τ

0.75

Bare Analytic CD LCD (spectral)

|W diss (t)|

Bare LCD

1.00

(a) Dissipated work for overdamped FPE evolved under DW potential for τ = 0.1

0.00

0.25

0.50

t/τ

0.75

1.00

(b) Dissipated work for overdamped FPE evolved under harmonic potential for τ = 0.1

FIG. 5: Dissipation during Fokker–Planck time-evolution. Dissipated work magnitude plotted against normalized time t/τ . Bare (uncontrolled) evolution (solid orange line), LCD (dashed grey line) and if available, Analytic CD (solid green line). TABLE II: Maximum dissipated work magnitude max |Wdiss (ζt )| calculated for the bare and LCD-driven protocols across the double-well and harmonic potentials at varying driving periods. Double-well potential Protocol duration (τ ) −3

−7

10 (∆t = 10 ) 0.01 (∆t = 10−6 ) 0.10 (∆t = 10−5 ) 0.25 (∆t = 10−5 ) 0.50 (∆t = 10−5 ) 0.75 (∆t = 10−5 ) 1.00 (∆t = 10−5 ) 2.00 (∆t = 10−5 ) 4.00 (∆t = 10−4 ) 6.00 (∆t = 10−4 ) 10.0 (∆t = 10−4 )

Bare 1.40 1.40 1.37 1.32 1.26 1.20 1.15 0.99 0.76 0.62 0.43

LCD −12

8.08 × 10 1.94 × 10−12 2.85 × 10−12 1.40 × 10−12 1.02 × 10−12 1.88 × 10−12 2.08 × 10−12 1.55 × 10−12 1.23 × 10−13 3.95 × 10−13 9.82 × 10−13

puting free-energy simulations, where the spectral resolution needed to achieve a target tracking fidelity can be reduced with improved convergence compared to fully uncontrolled schemes.

Harmonic trap Bare

Spectral framework and vanishing spectral gap

As a stress test of the spectral framework, we apply the LCD to the quartic coalescence potential [8], p V (x, ζt ) = x4 − 16(1 − ζt )x2 , whose two minima at ± 8(1 − ζt ) merge at the origin as ζt → 1, annihilating a barrier of height [8(1 − ζt )]2 . The exact free energy difference is ∆F = 62.9407... at β = 1. Table S5 reports the terminal

0.80 0.79 0.72 0.61 0.48 0.39 0.33 0.19 0.10 0.07 0.04

LCD

Analytic −12

6.80 × 10 5.70 × 10−12 5.15 × 10−12 3.97 × 10−12 3.31 × 10−12 2.48 × 10−12 2.27 × 10−12 1.28 × 10−12 7.07 × 10−13 4.33 × 10−13 2.68 × 10−13

9.43 × 10−13 6.22 × 10−14 8.02 × 10−14 4.14 × 10−14 5.21 × 10−14 3.59 × 10−14 2.03 × 10−14 2.92 × 10−14 1.20 × 10−14 8.77 × 10−15 6.55 × 10−15

work and dissipation for protocol durations τ ∈ [0.1, 0.5]: the LCD consistently drives W(τ ) closer to ∆F than the bare dynamics, with improvements ranging from ∼ 33× at τ = 0.1 to ∼ 1400× at τ = 0.5. Nevertheless, the residual Wdiss for LCD remains nonnegligible even at moderate driving speeds, in sharp contrast to the double-well and harmonic potentials where |Wdiss | ≲ 10−12 . The origin of this degradation is visible in Fig. 8: the spectral gap |λ1 (t)| drops below 10−12 for t/τ ≈ 0.33–0.4, reflecting the exponentially slow inter-well relaxation in the presence of a barrier exceeding 16 kB T . In this regime the spectral coefficients ∼ 1/|λn | diverge, the eigenvector matrix D(t) becomes ill-conditioned (cond(D) > 108 ), and the biorthogonal

10 3

10 5

10 7

10 5 10 7

10 7

|W diss (t)|

D KL (ρkπeq )

TVD(ρ, πeq )

10 9

10 11

10 13

10 15

10 15 0.0

0.2

M=5

0.4

t/τ

0.6

0.8

M = 10

1.0

10 15 0.0

M = 15

0.2

0.4

M = 20

t/τ

0.6

0.8

1.0

M = 25

0.0

0.2

M = 30

0.4

t/τ

0.6

0.8

1.0

Full

M = 35

10 5

10 7

|W diss (t)|

D KL (ρkπeq )

TVD(ρ, πeq )

(a) Truncated LCD performance for DW potential.

10 9

10 11

10 13

10 15

10 15 0.0

M=5

0.2

0.4

t/τ

0.6

M = 10

0.8

1.0

M = 15

0.0

0.2

M = 20

0.4

t/τ

0.6

0.8

M = 25

1.0

0.0

M = 30

0.2

0.4

t/τ

0.6

M = 35

0.8

1.0

Full

(b) Truncated LCD performance for harmonic potential.

FIG. 6: Convergence metrics of the spectrally truncated counterdiabatic expansion as a function of the active mode cutoff M . (a) Evolution under the double-well potential and (b) the driven harmonic trap. For both physical landscapes, the tracking fidelity—quantified via the total variation distance (TVD), Kullback-Leibler divergence (DKL ), and transient dissipated work Wdiss (t)—exhibits rapid exponential convergence toward the machine-precision floor as M approaches the full basis limit (N = 80). This uniform decay demonstrates that high-fidelity non-adiabatic lag suppression is preserved even under substantial modal reduction, rendering large-scale system tracking computationally viable.

decomposition can no longer resolve the nearly degenerate slow modes at double precision. Here, we essentially highlight and distinguish the breakdown of our spectral LCD construction as opposed to the counterdiabatic driving itself. As the barrier grows the ground-state gap |λ1 (ζ)| ∼ e−β∆Vb (ζ) [44] falls below double precision for t/τ ≲ 0.4 [Figure 8]; the numerically returned λ1 ∼ 10−12 –10−13 there is round-off dominated rather than the true (exponentially smaller) gap, and the near-degenerate eigenvectors become ill-conditioned (cond(D) > 108 ). The modal coefficient ℓ1 (∂ L̂/∂ζ)r0 /λ1 is then a quotient of two round-off quantities and corrupts the assembled L̂CD [Eq. (34)]. Crucially, the exact closed form counterdiabatic operator (35), by contrast remains bounded throughout this protocol and tracks πeq (t) directly; the divergence is therefore a breakdown of the spectral representation of L̂CD , not of the counterdiabatic control itself. For the symmetric coalescence potential the closing mode is antisymmetric, so ℓ1 (∂ L̂/∂ζ)r0 = 0 identically by symmetry and the inverse spectral gap 1/λ1 (t) term multiplies a vanishing

numerator (see Section VII C for details). The symmetric coalescence thus exposes a failure of the spectral representation, not a physical limit on engineering a counterdiabatic control.

VIII.

DISCUSSION

We have developed an explicit spectral construction of the counterdiabatic generator for classical stochastic systems governed by the Fokker–Planck equation. The underlying biorthogonal formalism—the eigenbasis {rn (t), ℓTn (t)} with ℓT0 = 1T (probability conservation) and r0 = πeq —follows [14]; our contribution is in the assembly L̂CD (t) from these eigenpairs (different from the approach of [14], its exact reduction to a rank-one closed form [Eq. (35)] under detailed balance, and the numerical characterization of when that spectral assembly is reliable, especially for reducing/eliminating dissipated work. The resulting spectral formula for the counterdiabatic correction L̂CD (t) is formally identical to the Berry–

Inverse relaxation rate, 1/|λn (t)|

1015 1012 109

Bare FP LCD (spectral) ∆F (exact) Analytic CD

106 103

1/|λ1 (t)| 1/|λ2 (t)| 1/|λ3 (t)|

100

10 2

10 1

100

101

FIG. 7: Work versus protocol duration for the harmonic trap. Total work W(τ ) as a function of protocol duration τ for the bare (red, solid) and LCD (blue, square) evolutions, compared with the analytic CD (green, star). The horizontal line marks the exact free energy difference ∆F = 21 ln 4 ≈ 0.6931. Under bare dynamics, W(τ ) > ∆F for all finite τ , approaching ∆F only in the quasi-static limit. Both the spectral LCD and the analytic CD yield W LCD (τ ) = ∆F to within ∼ O(10−12 ) across all protocol durations tested, confirming that the escorting condition Wdiss (τ ) ≈ 0 holds irrespective of driving speed.

Demirplak–Rice transitionless driving prescription (also analytically shown in [14]) and earns its place as an alternative approach escorting the Fokker–Planck density and a diagnostic that ties the lag to the spectral gap. The numerical demonstrations on the double-well and harmonic trap confirm two key claims. First, under LCD the figure of merits such as TVD (∼ 12-orders of magnitude improvement) and KL (∼ 16-orders of magnitude improvement – machine precision in FP64) clearly demonstrate that we track the instantaneous equilibrium distribution at all times compared to the uncontrolled evolution. Moreover, we show that W LCD (t) = ∆F(ζt ) holds at every instant t ∈ [0, τ ], making the dissipated work Wdiss (t) ≈ 0 [O(10−12 )] throughout the protocol at arbitrary driving speed. Secondly, we present a spectral truncation scheme, with the number of required spectral modes increasing modestly for faster protocols and vanishing spectral gap potentials, establishing spectral truncation as a systematic way to reduce (if not completely eliminate) the non-adiabatic lag. Unlike standard flow-field protocols which obscure the underlying spectral features, our biorthogonal Liouvillian formulation explicitly exposes the physical source of dissipation: near critical points (such as the quartic coalescence), the spectral coefficients/construction become ill-conditioned

0.00

0.25

1/|λ4 (t)| 1/|λ5 (t)|

ζ = 0.5

W(τ)

1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 10 3

0.50

t/τ

0.75

1.00

FIG. 8: Inverse relaxation rate for the quartic potential. Inverse gap 1/|λn (t)| for the five slowest modes along the quartic coalescence protocol. The ground (n = 1) mode spans twelve orders of magnitude, from a round-off-dominated plateau (t/τ ≈ 0.33–0.4, where the true gap is below machine precision) down to O(1) as ζ → 1. Modes n = 2–5 remain well-conditioned (O(0.3–0.5)) throughout, confirming that only the ground-state gap ∆1 (t) ≡ |λ1 (t) − λ0 (t)| (λ0 = 0) is the bottleneck.

as the gap closes, while the exact rank-one correction stays bounded. This establishes that the CD driving control stays bounded and only the spectral representation fails due to numerics (conditioning-related).

IX.

CONCLUSION

We introduced Liouvillian counterdiabatic driving (LCD), a spectral method, adapting Berrys’ framework for eliminating non-adiabatic lag and performing escorted dynamics in classical stochastic systems driven by finite-time protocols. The method constructs the counterdiabatic correction by promoting to a matrixvalued approach and extracts the non-adiabatic contributions during the Fokker–Planck evolution directly from the biorthogonal eigenpairs (spectral properties) of the Fokker–Planck generator, by moving to the instantaneous eigenbasis (adiabatic frame) using a timedependent gauge-type transform. Thus, one requires no generalized coordinate transforms, learned representations, or iterative optimisation, and yields a closed-form spectral sum that admits systematic truncation to the slowest relaxation modes, depending on the potential and protocol speeds. Numerical verification on a double-well potential, a harmonic trap, confirms tracking of the instantaneous equilibrium to high precision, vanishing dissipated work Wdiss (t) ≈ 0 at all times, and robustness

16 across protocol speeds. The spectral framework retains three roles that Eq. (36) does not make manifest, and which motivate it as the organizing object of this numerical approach: (i) it exposes the AGC Γ = D −1 (∂D/∂t) as the mechanism that pumps probability out of the instantaneous equilibrium manifold (zero mode), and hence the spectral-gap origin of the lag; (ii) it yields a systematically truncable low-rank approximation; and (iii) offers an escorting method with natural extension to systems where closed form solutions are non-analytic [45] and must be obtained numerically. Thus, our framework provides a complementary route to escorting the Fokker–Planck equation alongside existing flow-field and diffeomorphism-based approaches.

LIMITATIONS AND FUTURE DIRECTIONS

The spectral LCD formula contains inverse relaxation rates 1/λn (t) in its coefficients and the spectral coefficients/construction become ill-conditioned as the gap closes, while the exact rank-one correction stays bounded. This is demonstrated explicitly on the quartic coalescence potential (Section VII E), where the LCD still reduces Wdiss by one to two orders of magnitude but cannot reach the machine-precision escorting achieved on the double-well and harmonic potentials, whose gaps remain well separated throughout the protocol. We note that closed-form solutions [40] or diffeomorphism-based approaches such as normalizing flows [25] which operate entirely in configuration space are therefore insensitive to the spectral gap of L̂(t); the quartic coalescence regime that limits the present spectral framework may thus remain tractable for flow-based methods, suggesting that the two approaches have complementary domains of applicability. Moreover, the O(N 3 ) cost of full diagonalization is mitigated by three complementary strategies: exploiting the symmetric form of the rate matrix [Eq. (22)], which admits tridiagonal solvers at 12–15× speedup (Section VII E) for our 1D use cases; spectral truncation to M < N (or even M ≪ N ) modes (Section VII D) at the expense of varied degrees of imperfect escorting accuracies depending on the mode truncation; and sparse Lanczos eigensolvers [46, 47] that exploit the (2d + 1)N nearest-neighbor sparsity of the Sasa–Tasaki discretization at a cost of O(M dN ) (Section VII F). Together, these establish a viable path toward higher-dimensional systems where full diagonalisation is infeasible. Looking ahead, several directions remain open: the most important direction is to extend our spectral methodology in the absence of the detail-balance condition. For such non-equilibrium steady state systems, the steady-state density πss (x, t) is not known analytically, and requires to be computed numerically. In this work, we took the first steps by implementing the spectral framework to known baselines, subsequently to tackle complex nonequilibrium steady state (NESS) sys-

tems as a future direction, which require careful consideration of the non-conservative forces and housekeeping (Hatano-Sasa) heat terms [48, 49]. Secondly, potentials which leads to vanishing spectral gaps, wherein timedependent near-degeneracies cause the biorthogonal decomposition to become ill-conditioned and demand regularisation strategies beyond the standard eigensolver, requiring diagonalization-free methodologies operating in Krylov spaces [50, 51], together with a careful analysis of the thermodynamic cost and interpretation of counterdiabatic driving in this regime. Thirdly, establishing the precise functional relationship between the eigenbasis gauge connection D −1 (∂D/∂t) and the escorting flow field u(x, t) [7] would clarify the similarities and differences between the two approaches. Finally, extending the framework to two- and three-dimensional potentials with non-trivial barrier structure, where sparse Lanczos-based truncated LCD can be benchmarked in regimes where full eigendecomposition is no longer available, is the natural next step toward practical applications.

Reproducibility statement

All simulations were carried out in Python (3.9.25) using NumPy (1.26.4), SciPy (1.13.1), and Matplotlib (3.9.3). NumPy’s dense eigendecomposition (numpy.linalg.eig), used throughout for the biorthogonal decomposition of the Liouville generator L̂(t), was linked against OpenBLAS 0.3.23 (dynamic architecture dispatch, DYNAMIC ARCH), which parallelizes the underlying LAPACK calls across available cores; wall-clock timings reported in this work (e.g. the dense-vssymmetric-tridiagonal speedup benchmarks) reflect this multi-threaded BLAS backend and will vary with core count and BLAS implementation on other hardware, though the relative scaling/speedup trends are expected to hold generally. Computations were run on a dual-socket workstation with two Intel Xeon Platinum 8468 processors (48 cores/socket, 2 threads/core; 192 logical CPUs total) and 2 TiB RAM, running Rocky Linux 9.8 (kernel 5.14.0687.17.1.el9 8.x86 64). No GPU acceleration was used; all reported computations are CPU-only. Code, a curated set of example results, and a demonLCD-FPE stration notebook are available at

ACKNOWLEDGEMENTS

The ELLIS Unit Linz, the LIT AI Lab, the Institute for Machine Learning, are supported by the Federal State Upper Austria. We thank the projects FWF AIRI FG 9-N (10.55776/FG9), AI4GreenHeatingGrids (FFG- 899943), Stars4Waters (HORIZON-CL6-2021CLIMATE-01-01). We thank NXAI GmbH, Audi AG, Silicon Austria Labs (SAL), Merck Healthcare KGaA,

17 GLS (Univ. Waterloo), TÜV Holding GmbH, Software Competence Center Hagenberg GmbH, dSPACE GmbH, TRUMPF SE + Co. KG. Sandeep Suresh

Cranganore was supported by the FWF Bilateral Artificial Intelligence initiative under Grant Agreement number 10.55776/COE12.

[1] T. Albash and D. A. Lidar, Adiabatic quantum computation, Rev. Mod. Phys. 90, 015002 (2018). [2] C. Jarzynski, Equalities and inequalities: irreversibility and the second law of thermodynamics at the nanoscale, Annu. Rev. Condens. Matter Phys. 2, 329 (2011). [3] T. Schäfer and E. V. Shuryak, Instantons in qcd, Rev. Mod. Phys. 70, 323 (1998). [4] C. Jarzynski, Nonequilibrium equality for free energy differences, Phys. Rev. Lett. 78, 2690 (1997). [5] R. M. Neal, Annealed importance sampling (1998), arXiv:physics/9803008 [physics.comp-ph]. [6] G. E. Crooks, Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences, Phys. Rev. E 60, 2721 (1999). [7] S. Vaikuntanathan and C. Jarzynski, Escorted free energy simulations: Improving convergence by reducing dissipation, Phys. Rev. Lett. 100, 190601 (2008). [8] S. Lee and C. Jarzynski, Estimating free energy differences with virtually escorted trajectories (2026), arXiv:2606.30451 [cond-mat.stat-mech]. [9] D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martı́nez-Garaot, and J. G. Muga, Shortcuts to adiabaticity: concepts, methods, and applications, Rev. Mod. Phys. 91, 045001 (2019). [10] M. Kolodrubetz, D. Sels, P. Mehta, and A. Polkovnikov, Geometry and non-adiabatic response in quantum and classical systems, Phys. Rep. 697, 1 (2017). [11] M. V. Berry, Transitionless quantum driving, Journal of Physics A: Mathematical and Theoretical 42, 365303 (2009). [12] M. Demirplak and S. A. Rice, Adiabatic population transfer with control fields, J. Phys. Chem. A 107, 9937 (2003). [13] A. del Campo, Shortcuts to adiabaticity by counterdiabatic driving, Phys. Rev. Lett. 111, 100502 (2013). [14] S. Iram, E. Dolson, J. Chiel, J. Pelesko, N. Krishnan, Ö. Güngör, B. Kuznets-Speck, S. Deffner, E. Ilker, J. G. Scott, and M. Hinczewski, Controlling the speed and trajectory of evolution with counterdiabatic driving, Nature Physics 17, 135 (2021). [15] S. Vaikuntanathan and C. Jarzynski, Escorted free energy simulations, The Journal of Chemical Physics 134, 054107 (2011). [16] C. Jarzynski, Generating shortcuts to adiabaticity in quantum and classical dynamics, Phys. Rev. A 88, 040101(R) (2013). [17] A. Patra and C. Jarzynski, Shortcuts to adiabaticity using flow fields, New J. Phys. 19, 125009 (2017). [18] G. Li, H. T. Quan, and Z. C. Tu, Shortcuts to isothermality and nonequilibrium work relations, Phys. Rev. E 96, 012144 (2017). [19] C. Jarzynski, Targeted free energy perturbation, Phys. Rev. E 65, 046122 (2002). [20] P. Wirnsberger, A. J. Ballard, G. Papamakarios, S. Abercrombie, S. Racanière, A. Pritzel, D. Jimenez Rezende, and C. Blundell, Targeted free energy estimation via

learned mappings, The Journal of Chemical Physics 153, 144112 (2020). [21] M. R. Shirts and J. D. Chodera, Statistically optimal analysis of samples from multiple equilibrium states, The Journal of Chemical Physics 129, 124105 (2008). [22] Y. Lipman, R. T. Q. Chen, H. Ben-Hamu, M. Nickel, and M. Le, Flow matching for generative modeling (2023), arXiv:2210.02747 [cs.LG]. [23] D. J. Rezende and S. Mohamed, Variational inference with normalizing flows (2016), arXiv:1505.05770 [stat.ML]. [24] D. Nielsen, P. Jaini, E. Hoogeboom, O. Winther, and M. Welling, Survae flows: Surjections to bridge the gap between vaes and flows (2020), arXiv:2007.02731 [cs.LG]. [25] H. Wu, J. Köhler, and F. Noé, Stochastic normalizing flows, in Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 612, 2020, virtual , edited by H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin (2020). [26] L. Holdijk, N. M. Anand, M. M. Bronstein, and M. Welling, Learning escorted protocols for multistate free-energy estimation, in The Fourteenth International Conference on Learning Representations (2026). [27] F. Vargas, S. Padhy, D. Blessing, and N. Nüsken, Transport meets variational inference: Controlled monte carlo diffusions, in The Twelfth International Conference on Learning Representations (2024). [28] M. S. Albergo and E. Vanden-Eijnden, Nets: A nonequilibrium transport sampler (2025), arXiv:2410.02711 [cs.LG]. [29] P. Holderrieth, M. S. Albergo, and T. Jaakkola, Leaps: A discrete neural sampler via locally equivariant networks, in Proceedings of the 42nd International Conference on Machine Learning (2025). [30] J. Dieudonné, On biorthogonal systems., Michigan Mathematical Journal 2, 7 (1953). [31] D. C. Brody, Biorthogonal quantum mechanics, Journal of Physics A: Mathematical and Theoretical 47, 035305 (2013). [32] M. V. Berry, Transitionless quantum driving, J. Phys. A 42, 365303 (2009). [33] C. Jarzynski, Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach, Phys. Rev. E 56, 5018 (1997). [34] G. Hummer and A. Szabo, Free energy reconstruction from nonequilibrium single-molecule pulling experiments, Proceedings of the National Academy of Sciences 98, 3658 (2001), https://www.pnas.org/doi/pdf/10.1073/pnas.071034098. [35] S.-I. Sasa and H. Tasaki, Steady state thermodynamics, Journal of Statistical Physics 125, 125 (2006). [36] E. Ilker, O. m. c. Güngör, B. Kuznets-Speck, J. Chiel, S. Deffner, and M. Hinczewski, Shortcuts in stochastic systems and control of biophysical processes, Phys. Rev. X 12, 021048 (2022).

1 [37] A. C. Santos and M. S. Sarandy, Generalized transition[47] G. Golub and C. Van Loan, Matrix Computations, Johns less quantum driving for open quantum systems, Phys. Hopkins Studies in the Mathematical Sciences (Johns Rev. A 104, 062421 (2021). Hopkins University Press, 1996). [38] G. Vacanti, R. Fazio, S. Montangero, G. M. Palma, [48] U. Seifert, Entropy production along a stochastic traM. Paternostro, and V. Vedral, Transitionless quantum jectory and an integral fluctuation theorem, Phys. Rev. driving in open quantum systems, New Journal of Physics Lett. 95, 040602 (2005). 16, 053017 (2014). [49] T. Hatano and S.-i. Sasa, Steady-state thermodynamics [39] M. Demirplak and S. A. Rice, Assisted adiabatic passage of langevin systems, Phys. Rev. Lett. 86, 3463 (2001). revisited, The Journal of Physical Chemistry B 109, 6838 [50] S. Morawetz and A. Polkovnikov, Universal counterdia(2005). batic driving in krylov space, PRX Quantum 6, 040320 [40] D. Guéry-Odelin, C. Jarzynski, C. A. Plata, A. Prados, (2025). and E. Trizac, Driving rapidly while remaining in control: [51] A. W. Shrestha, B. Bhattacharjee, and A. del Campo, classical shortcuts from hamiltonian to stochastic dynamShortcuts to adiabaticity for non-hermitian systems in ics, Reports on Progress in Physics 86, 035902 (2023). krylov space (2026), arXiv:2607.07802 [quant-ph]. [41] A. Zhong and M. R. DeWeese, Beyond linear response: [52] W. Magnus, On the exponential solution of differential Equivalence between thermodynamic geometry and opequations for a linear operator, Commun. Pure Appl. timal transport, Phys. Rev. Lett. 133, 057102 (2024). Math. 7, 649 (1954). [42] S. ichi Sasa and H. Tasaki, Steady state thermodynamics [53] S. Blanes, F. Casas, J. A. Oteo, and J. Ros, The Magnus (2005), arXiv:cond-mat/0411052 [cond-mat.stat-mech]. expansion and some of its applications, Phys. Rep. 470, [43] S. Campbell and S. Deffner, Trade-off between speed and 151 (2009). cost in shortcuts to adiabaticity, Phys. Rev. Lett. 118, [54] N. J. Higham, The scaling and squaring method for 100601 (2017). the matrix exponential revisited, SIAM Journal on [44] M. G. Evans and M. Polanyi, Some applicaMatrix Analysis and Applications 26, 1179 (2005), tions of the transition state method to the calhttps://doi.org/10.1137/04061101X. culation of reaction velocities, especially in so[55] G. Baker and P. Graves-Morris, Padé Approximants, Enlution, Transactions of the Faraday Society cyclopedia of Mathematics and its Applications (Cam31, 875 (1935), https://pubs.rsc.org/tf/articlebridge University Press, 1996). pdf/doi/10.1039/TF9353100875/1608808/tf9353100875.pdf. [56] I. A. Martı́nez, A. Petrosyan, D. Guéry-Odelin, E. Trizac, [45] B. Derrida, Non-equilibrium steady states: fluctuations and S. Ciliberto, Engineered swift equilibration of a and large deviations of the density and of the current, brownian particle, Nature Physics 12, 843 (2016). Journal of Statistical Mechanics: Theory and Experi[57] K. J. Laidler and M. C. King, Development of transitionment 2007, P07023 (2007). state theory, The Journal of Physical Chemistry 87, 2657 [46] C. Lanczos, An iteration method for the solution of the (1983). eigenvalue problem of linear differential and integral op[58] H. Oberhofer, C. Dellago, and P. L. Geissler, Biased samerators, Journal of Research of the National Bureau of pling of nonequilibrium trajectories: can fast switching Standards 45, 255 (1950). simulations outperform conventional free energy calculation methods?, The journal of physical chemistry. B 109, 6902 (2005).

SUPPLEMENTAL MATERIAL I.

ZERO-VARIANCE FREE ENERGY ESTIMATION VIA PERFECT BIJECTIVE MAPPINGS

For an arbitrary set of invertible, bijective coordinate mapping functions {Mi }, the statistical efficiency of utilizing escorted simulations to estimate the free energy difference ∆F depends crucially on minimizing the relative entropy between the driven and target ensembles. Because the convergence of exponential averages—such as the free energy estimator in Eq. (1)—deteriorates rapidly with increasing non-equilibrium dissipation [7, 15], the statistical variance is fundamentally governed by the phase-space lag between the instantaneous non-equilibrium distribution and the corresponding equilibrium target. Consequently, an optimal choice of coordinate transformations must actively suppress this lag to maximize phase-space overlap and optimize estimator convergence. To formalize this optimization framework, one analyzes the limiting case of a “perfect” set of mapping functions, denoted by {M∗i }, which entirely eliminates the non-equilibrium lag. Formally, for an ensemble of trajectories initiated ζ0 from the canonical equilibrium state π0 ∼ πeq (x0 ) [since the evolution is described for discrete time-stamps ζti = ζi , ζ0 one succinctly denotes πeq (x0 , ζ0 ) = πeq (x0 )] and evolved under the escorted dynamics, a perfect map ensures that ζi the subsequent microstates xi strictly trace the instantaneous equilibrium manifold πeq (xi ) for all discrete steps 1 ≤ i < N . Geometrically, this requires the target and escorted probability densities to coincide identically. Under a bijective coordinate transformation M∗i : x → x′ , conservation of probability requires the initial distribution ζi πeq (x) to transform under the push-forward map as [19, 25]: ζi+1 (x′ ) = πeq

ζi (x) πeq , ∗ Jζi (x)

(S1)

where Jζ∗i (x) = det |∂x′ /∂x| denotes the Jacobian determinant of the transformation. ζt (x) ≡ πeq (x, ζt ) = e−β V (x,ζt )−F (ζt ) into Eq. (S1) and taking the Substituting the canonical Boltzmann weight πeq natural logarithm of both sides yields the exact microscopic energy balance for a perfect mapping: δWiesc ≡ Vζi+1 (x′ ) − Vζi (x) − β −1 log Jζ∗i (x) = Fζi+1 − Fζi .

(S2)

where W esc corresponds to the escorted work. Summing this relation over the full discrete protocol reveals that the total accumulated work with the Jacobian determinant scaled density (escorting) along any arbitrary phase-space trajectory γ collapses identically to the total free energy difference, W esc [γ] = ∆F. Thermodynamically, this condition suppresses all sample-to-sample fluctuations in the work distribution, forcing the non-equilibrium work probability density to contract to a singular Dirac delta function: PF (W ) = δ(W − ∆F).

(S3)

Crucially, for a perfect set of mappings {M∗i }, the exponential estimator [Eq. (2)] yields a deterministic, zero-variance evaluation of ∆F from a single trajectory invocation. Under these conditions, i.e., a perfect sequence of mappings for the forward protocol, an identical zero-variance condition is structurally preserved in the time-reversed process.

II.

DISCRETIZATION FRAMEWORK AND THE CONTINUUM LIMIT

We consider the overdamped stochastic evolution of a classical system in a time-varying potential V (x, t). In the continuum limit, the probability density ρ(x, t) evolves according to the Fokker-Planck equation (FPE), which can be written compactly in terms of the continuous Liouville operator as: ∂ρ(x, t) ∂ = ∂t ∂x

∂V (x, t) ∂x

∂ 2 ρ(x, t) ρ(x, t) + β −1 , ∂x2

(S4)

where β = 1/(kB T ) represents the inverse temperature of the thermal reservoir. To treat this evolution numerically −1 and analyze its spectral properties, we project the continuous spatial coordinate onto a uniform lattice {xi }N i=0 with an associated grid spacing ∆x. The continuous probability density is mapped onto a discrete state P vector ρ(t) = (ρ0 , . . . , ρN −1 )T ∈ RN , where ρi (t) ≈ ρ(xi , t)∆x, subject to the strict normalization constraint i ρi = 1.

3 Under this projection, the spatial partial differential equation maps onto a discrete Master Equation governing the probability exchange via local fluxes: ∂ρi = kfwd [i − 1]ρi−1 + kbwd [i]ρi+1 − kfwd [i] + kbwd [i − 1] ρi , ∂t

(S5)

where kfwd [i] and kbwd [i] represent the transition rates out of site i to neighboring sites i + 1 and i − 1, respectively. To preserve local detailed balance with respect to the instantaneous Boltzmann distribution πeq (x, t) ∝ exp[−βV (x, t)], we select Sasa-Tasaki rates defined as: 1 β exp − ∆Vi , kfwd [i] = (S6) β∆x2 2 β 1 exp + ∆Vi , kbwd [i] = (S7) β∆x2 2 where ∆Vi ≡ Vi+1 − Vi denotes the forward potential difference across the lattice bond. To rigorously demonstrate that the discrete version structurally converges to the continuous FPE, we perform a multivariable Taylor expansion of the fields ρ(x, t) and V (x, t) around the reference grid site xi up to third order in ∆x. Let ∂f /∂x ≡ f ′ and ∂ 2 f /∂x2 ≡ f ′′ . The non-local spatial shifts and potential gradients map to local expansions according to: ρi±1 = ρ ± ∆xρ′ +

∆x2 ′′ ρ ± O(∆x3 ), 2

∆x2 ′′ V + O(∆x3 ), 2 ∆x2 ′′ ∆Vi−1 = ∆xV ′ − V + O(∆x3 ). 2 ∆Vi = ∆xV ′ +

(S8) (S9) (S10)

Expanding the exponential terms in the Sasa-Tasaki rate definitions to second order in ∆x yields the incoming rate operators feeding probability into site i: β ′′ β 2 ′ 2 1 β∆x ′ 2 kfwd [i − 1] = V + ∆x V + (V ) , (S11) 1− β∆x2 2 4 8 1 β ′′ β 2 ′ 2 β∆x ′ 2 kbwd [i] = V + ∆x V + (V ) . (S12) 1 + β∆x2 2 4 8 Taking the respective products of these incoming rates with the shifted state densities ρi±1 , we evaluate the total incoming probability flux up to order O(∆x0 ): kfwd [i − 1]ρi−1 + kbwd [i]ρi+1 β ′′ β2 ′ 2 1 2 ′′ ′ ′ = 2ρ + ∆x (ρ + βV ρ + V ρ + (V ) ρ β∆x2 2 4 +O(∆x). Correspondingly, we evaluate the outgoing escape rates moving away from site i: 1 β∆x ′ β∆x2 ′′ β 2 ∆x2 ′ 2 kfwd [i] = 1 − V − V + (V ) , β∆x2 2 4 8 1 β∆x ′ β∆x2 ′′ β 2 ∆x2 ′ 2 kbwd [i − 1] = V − V + (V ) . 1 + β∆x2 2 4 8 Summing these rate coefficients establishes the total escape rate acting on the local state density ρi : kfwd [i] + kbwd [i − 1] ρi 1 β ′′ β2 ′ 2 2 = 2ρ + ∆x − V ρ + (V ) ρ . β∆x2 2 4

(S13)

(S14) (S15)

(S16)

4 Finally, substituting the explicit incoming flux expansion Eq. (S13) and outgoing escape expansion Eq. (S16) back into the master equation structure Eq. (S5), the singular leading-order term 2ρ/(β∆x2 ) cancels out identically. Grouping the surviving structural components yields: ∂ρ = V ′ ρ′ + V ′′ ρ + β −1 ρ′′ . ∂t Recognizing that V ′ ρ′ + V ′′ ρ ≡ ∂x (V ′ ρ), we reconstruct the derivative components to obtain: ∂ ∂V (x, t) ∂ 2 ρ(x, t) ∂ρ(x, t) . = ρ(x, t) + β −1 ∂t ∂x ∂x ∂x2

(S17)

(S18)

This confirms that the discrete grid network updates perfectly recover the continuous Fokker-Planck dynamics in the continuum limit ∆x → 0. III.

SPECTRAL REPRESENTATION OF THE GAUGE CONNECTION

To establish the exact equivalence between the the AGC term D −1 (∂D/∂t) and the spectral-gap-dependent Berrytype formulation, we begin with the instantaneous eigenvalue equation for the right eigenvectors of the Liouvillian operator, L̂(t)rm (t) = λm (t)rm (t). Differentiating Eq. (S19) with respect to time t yields ∂ L̂/∂t rm + L̂ ∂rm /∂t = ∂λm /∂t rm + λm ∂rm /∂t .

(S19)

(S20)

We now project Eq. (S20) onto the dual space by left-multiplying with the left eigenvector ℓTn (t) corresponding to the index n ̸= m. This yields ℓTn ∂ L̂/∂t rm + ℓTn L̂ ∂rm /∂t = ∂λm /∂t ℓTn rm + λm ℓTn ∂rm /∂t . (S21) Invoking the left-eigenvalue relation ℓTn L̂ = λn ℓTn and the biorthogonality condition ℓTn rm = δnm , the first term on the right-hand side vanishes for n ̸= m. Simplifying the left-hand side leads to ℓTn ∂ L̂/∂t rm + λn ℓTn ∂rm /∂t = λm ℓTn ∂rm /∂t . (S22) Rearranging terms directly isolates the projection of the eigenvector velocity onto the reciprocal basis, ℓTn ∂ L̂/∂t rm T ℓn ∂rm /∂t = (n ̸= m). λm − λn

(S23)

The left-hand side of Eq. (S23) defines the off-diagonal elements of the geometric gauge connection matrix in the adiabatic moving frame, Γnm ≡ ℓTn ∂t rm . This result highlights a key physical insight: the gauge connection Γ(t) = D−1 ∂t D is fundamentally a spectral-gapdependent object. When computed numerically, the inverse relaxation gaps (λm − λn )−1 are implicitly encoded in the gradients of the state-space trajectory ∂rm /∂t, which diverge near avoided crossings or spectral degeneracies. To express the full non-adiabatic correction in the laboratory frame, we map the off-diagonal components of the connection back via the similarity transformation DΓoff D−1 . Utilizing the completeness relation of the biorthogonal basis, we obtain X X ℓTn ∂ L̂/∂t rm rn ℓTm . (S24) DΓoff D−1 = Γnm rn ℓTm = λm − λ n n̸=m

n̸=m

The spectral counterdiabatic driving term L̂CD (t) corresponds to the instantaneous equilibrium distribution m = 0 column of the full gauge connection. Given that λ0 = 0, setting m = 0 in Eq. (S24) recovers the spectral correction X ℓTn (t) ∂ L̂/∂t r0 (t) L̂CD (t) = − rn (t)ℓT0 , (S25) λn (t) n̸=0

which is mathematically identical to the localized Berry formulation, as adapted here to the Fokker–Planck evolution.

5 IV. A.

NUMERICS SPECIFICS

Eigendecomposition of the Liouville Operator

At each instant t, one can perform exact diagonlization of the rate matrix: L̂(t) = D(t) Λ(t) D −1 (t),

(S26)

where D(t) ∈ RN ×N contains the right eigenvectors rn as its columns, D −1 (t) has the left eigenvectors as its rows, and Λ(t) = Diag(λ0 , λ1 , . . . , λN −1 ),

0 = λ0 > λ1 ≥ λ2 ≥ · · ·

(S27)

We verify numerically that, D −1 L̂D − Λ

∞

< 10−6

∀ protocol values visited during the anneal,

(S28)

confirming exact diagonalisation. We report the numerics for the exact diagonalized Liouville operator, for both systems studied in this work, below.

Exact diagonalization verification Double well

Problem setup. For our experimental model: the external time-varying potential V (x, ζt ) = x4 − 2x2 + ζt x on a grid of N = 80 points over x ∈ [−2.5, 2.5], β = 1. The protocol is a smooth sigmoid from ζi = −1 to ζf = +1: ζ(t) = ζi + (ζf − ζi ) s3 (6s2 − 15s + 10),

s = t/τ,

(S29)

evaluated at t/τ = 0, 0.1, 0.25, 0.5, 0.75, 1. Since ζ(t) depends only on s = t/τ , this check is independent of the total anneal duration τ . Verification of D −1 L̂D = Λ. Table S1 reports the off-diagonal norm of D −1 L̂D and the biorthogonality residual ∥D −1 D−I∥∞ at the six snapshots. Both are at floating-point precision (10−7 –10−11 ), confirming Eq. (S28). Figure S1 shows L̂(t) itself (top row – the tridiagonal birth–death structure is visible throughout, magnitude growing toward the ends of the anneal where the wells are most separated) alongside log10 |D −1 L̂D − Λ| (bottom row), confirming the residual shows no structure beyond numerical noise at every snapshot. TABLE S1: Double well: numerical verification of the biorthogonal diagonalisation at six points along the anneal. εΛ = ∥D −1 L̂D − Λ∥∞ (reconstruction residual); εI = ∥D −1 D − I∥∞ (biorthogonality residual). t/τ 0.00 0.10 0.25 0.50 0.75 1.00

ζ −1.0000 −0.9829 −0.7930 0.0000 0.7930 1.0000

εΛ 1.00 × 10−7 1.65 × 10−7 7.85 × 10−8 2.80 × 10−8 3.93 × 10−8 2.49 × 10−8

εI 1.60 × 10−10 1.00 × 10−10 6.35 × 10−11 2.55 × 10−11 1.51 × 10−11 1.17 × 10−11

t = 0.00 τ ζ = − 1.00

t = 0.10 τ ζ = − 0.98

t = 0.25 τ ζ = − 0.79

t = 0.50 τ ζ = 0.00

t = 0.75 τ ζ = 0.79

t = 1.00 τ ζ = 1.00

1000

L(t)

500 0 500 1000

log 10 |D −1 LD − Λ|

6 7 8 9 10 11 12 13 14

FIG. S1: Double well: L̂(t) (top) and the diagonalization residual log10 |D −1 L̂D − Λ| (bottom) at six points along the anneal. The residual is uniformly at the 10−14 –10−6 noise floor at every snapshot, with no visible structure.

Harmonic trap

Problem setup. The harmonic trap, V (x, t) = 12 κ0 (t)x2 , on a grid of N = 80 points over x ∈ [−4, 4], β = 1, annealed via the same smoothstep schedule Eq. (S29) with κ0 : 1 → 4, evaluated at the same six t/τ points. Verification of D −1 L̂D = Λ. Table S2 reports the same two residuals at the same six snapshots. Both remain at floating-point precision throughout (10−8 –10−15 ), though εΛ and εI grow by roughly four orders of magnitude from κ0 = 1 to κ0 = 4 – visible directly in Fig. S2’s residual row, which brightens monotonically left to right. This tracks the generator’s overall matrix norm growing with κ0 (a stiffer trap gives larger transition rates kf , kb ∝ Dc e∓β∆V /2 ), which sets the scale of floating-point rounding in D −1 L̂D; the residual remains many orders of magnitude below any physically relevant scale throughout. TABLE S2: Harmonic trap: numerical verification of the biorthogonal diagonalisation at six points along the protocol. t/τ 0.00 0.10 0.25 0.50 0.75 1.00

κ0 1.0000 1.0257 1.3105 2.5000 3.6895 4.0000

εΛ 4.36 × 10−12 5.71 × 10−12 1.17 × 10−11 6.55 × 10−10 3.88 × 10−8 9.82 × 10−8

εI 5.24 × 10−15 4.30 × 10−15 9.60 × 10−15 5.65 × 10−13 1.74 × 10−10 2.04 × 10−10

t = 0.10 τ 0 = 1.03

t = 0.25 τ 0 = 1.31

t = 0.50 τ 0 = 2.50

t = 0.75 τ 0 = 3.69

t = 1.00 τ 0 = 4.00

200 100

L(t)

t = 0.00 τ 0 = 1.00

0 100 200

log 10 |D −1 LD − Λ|

6 7 8 9 10 11 12 13 14

FIG. S2: Harmonic trap: L̂(t) (top) and the diagonalization residual log10 |D −1 L̂D − Λ| (bottom) at six points along the anneal. The residual grows visibly but remains at floating-point precision throughout, tracking the growth of κ0 . V.

DYNAMICS IN THE LABORATORY FRAME

The transformation from the laboratory frame to the adiabatic frame of reference via a change of basis, ρ̃(t) = D −1 (t) ρ(t), does not produce a probability vector. Its components ρ̃n (t) = ℓTn ρ(t) are projections of ρ onto the left eigenvectors of L̂(t) and are not constrained to be non-negative or to sum to unity. Thus, the normalization condition, P i.e. i ρ̃i = 1, need not hold in this representation. Nevertheless, the normalisation of ρ in the original frame imposes a precise constraint on ρ̃, which we now derive. A.

Dynamics under change of basis

Separating the n = 0 term in Eq. (23) and using ρ̃0 = 1 and r0 (t) = πeq (t): N −1 X ρ̃n (t) rn (t). ρ(t) = πeq (t) + | {z } n=1

(S30)

ρ̃0 =1

Thus, the vector-valued density ρ(t) can be viewed as a linear combination of relaxation modes weighted by their excitation amplitudes ρ̃n (t). The first term is the instantaneous equilibrium distribution; the sum is over higher modes that will have larger excitation amplitudes induced by finite-time protocols ζt during the Fokker–Planck evolution. The leakage to relaxation modes other than πeq (t) are responsible for non-adiabatic lag. Since 1T rn = 0 for n ≥ 1 (biorthogonality relation), the deviation term carries zero net probability: 1T

N −1 X n=1

ρ̃n (t)rn =

N −1 X

ρ̃n (t) 1T rn = 0,

(S31)

n=1

confirming that ρ(t) in the laboratory frame remains normalised for any values of ρ̃n . The component ρ̃n (t) = ℓTn ρ(t) measures the projection of the current distribution onto the n-th relaxation mode, weighted by the corresponding left eigenvector. By biorthogonality Eq. (17): ρ̃n (t) = 0 ⇐⇒ ρ(t) ⊥bio rn ,

(S32)

where ⊥bio denotes orthogonality in the biorthogonal sense. In particular, if ρ(t) = πeq (t) = r0 (t), then: ρ̃n (t) = ℓTn r0 = δn0 . Thus, for n ≥ 1 ρ̃n = ℓTn ρ is the projection of ρ onto relaxation mode n. T

(S33)

It can be positive, negative, or zero. It is not a probability — it is a mode amplitude. So ρ̃(t) = (1, 0, 0, . . . , 0) when and only when ρ(t) = πeq (t).

8 B.

Counterdiabatically evolved Liouville operator

Bare dynamics (no CD correction): The adiabatic gauge potential term in Eq. (27) Γnm is non-zero for n ̸= m, continuously pumping excitation from the zero mode into the relaxation modes, thus yielding an excess non-zero work in the form of dissipation, i.e. Wdiss = W − ∆F > 0. In the long-time driven steady state, ρ̃n ̸= 0 for n ≥ 1, meaning ρ(t) ̸= πeq (t) and the system lags behind equilibrium: X ρ(t) = πeq (t) + ρ̃n (t)rn (t), ρ̃n ̸= 0. (S34) n≥1

The scalar-valued version of Eq. (S34) was introduced in [14]. On adding the matrix-valued CD control terms LCD , one exactly cancels Γ in Eq. (27) which hinders the adiabatic tracking of the instantaneous equilibrium distribution. This yields the following evolution in the adiabatic frame: ∂ ρ̃n = λn ρ̃n , n ≥ 1. (S35) ∂t For a system initialized in equilibrium, ρ(0) = πeq (0), Eq. (S35) then indicates that ρ̃n (t) = 0 for all t ∈ [0, τ ], and therefore: ρ(t) = πeq (t),

∀ t ∈ [0, τ ].

(S36)

Under finite-speed bare driving, ρ̃n≥1 (t) ̸= 0 develops as the protocol pumps amplitude into the relaxation modes; the counterdiabatic correction cancels this pumping exactly, preserving ρ̃n≥1 (t) = 0 for all t ∈ [0, τ ]. Thus, the CD correction preserves the initial condition ρ̃ = (1, 0, . . . , 0)T exactly, ensuring perfect tracking of the instantaneous equilibrium at arbitrary driving speed. Probability conservation of L̂CD (t). Since ℓT0 = 1T and 1T rn = ℓ0 rn = δ0n = 0 for n ̸= 0, every column sum of L̂CD vanishes: 1T L̂CD = 0T . The corrected generator L̂eff := L̂ + L̂CD . Thus, it is a probability conserving quantity. We emphasize, however, that L̂eff is not a stochastic rate matrix: L̂CD = (∂πeq /∂t)1T is dense with sign-indefinite off-diagonal entries, so L̂eff generally has negative off-diagonals. C.

Magnus fourth-order integrator

Time integration of the modified Fokker–Planck equation dρ(t) = L̂(t) + L̂CD (t) ρ(t), (S37) dt requires a method that respects the structure of the generator at each step. The standard midpoint (Euler–exponential) scheme ρ(t + ∆t) = eG(t+∆t/2)∆t ρ(t), where G = L̂ + L̂CD , is accurate to O(∆t2 ) but becomes insufficient when the generator varies rapidly within a single step, as is the case near the midpoint of a fast protocol where L̂CD (t) is largest. We therefore adopt the fourth-order Magnus integrator [52, 53], which achieves O(∆t4 ) accuracy while preserving probability normalisation exactly at each step. The Magnus integrator advances ρ from t to t + ∆t via ρ(t + ∆t) = eΩ(t, t+∆t) ρ(t), where the Magnus exponent is approximated using the two Gauss–Legendre quadrature nodes √ ! 1 3 t1,2 = t + ∓ ∆t 2 6 as ∆t Ω(t, t + ∆t) = G1 + G2 + 2

√

3 ∆t2 G2 , G 1 , 12

(S38)

(S39)

(S40)

with Gk = L̂(tk ) + L̂CD (tk ) evaluated at the quadrature nodes, and [G2 , G1 ] = G2 G1 − G1 G2 the matrix commutator. The matrix exponential eΩ is computed at each step via the Padé approximation [54, 55], as implemented in scipy.linalg.expm. The commutator term in Eq. (S40) captures the leading non-commutativity of the generator between t1 and t2 , which is the dominant source of error in the midpoint scheme. For a slowly varying protocol the commutator is small and Eq. (S40) reduces to the midpoint scheme; for fast protocols, where L̂CD changes appreciably within a step, the commutator correction is essential.

9 D.

Analytic evolution for the harmonic trap

For the harmonic potential V (x, κ0 (t)) = 21 κ0 (t)x2 , the Fokker–Planck equation admits an exact reduction that bypasses the matrix-valued Liouville evolution entirely. Since a Gaussian initial condition remains Gaussian under harmonic confinement, we write r α(t) −α(t)x2 e , (S41) ρ(x, t) = π where α(t) > 0 is the time-dependent precision parameter. Substituting Eq. (S41) directly into the Fokker–Planck 2 ∂(κxρ) 1 ∂ ρ equation ∂ρ + βγ ∂t = ∂x ∂x2 and requiring the equality to hold for all x yields the scalar-valued Riccati equation α̇(t) = −

4α2 (t) 2κ(t) α(t) − , γ βγ

(S42)

with initial condition α(0) = βκi /2, corresponding to the equilibrium distribution for the initial stiffness κi . For the bare protocol κ(t) = κ0 (t), the solution α(t) lags behind the instantaneous equilibrium target αeq (t) = βκ0 (t)/2. For the analytic counterdiabatic protocol of Ref. [17, 56], replacing κ0 (t) with κCD (t) = κ0 (t) + γ κ̇0 (t)/(2κ0 (t)) enforces α(t) = αeq (t) exactly for all t ∈ [0, τ ]. Equation (S42) is integrated using a high-order adaptive solver (scipy.solve ivp, DOP853, rtol = 10−12 , atol = 10−14 ), with dense output evaluated at every reporting time. This approach carries no spatial discretisation error: ρ(x, t) is represented analytically via α(t) at every instant and projected onto the spatial grid {xi } only when evaluating scalar diagnostics such as the TVD, KL divergence, and dissipated work. The residual |α(t) − αeq (t)| reaches at most ∼ 10−11 , consistent with the ODE solver tolerance, and is entirely free of the O(∆x2 ) discretisation error that affects the matrix-valued spectral LCD on a finite grid. The harmonic trap therefore serves as a clean benchmark: the analytic reference isolates the O(∆x2 ) grid error in the spectral CD condition and the O(∆t2 ) heat-quadrature error in the dissipated work, neither of which is present in the scalar evolution. E.

Numerical evaluation of excess dissipated work

At each Magnus step tk (Appendix V C), the instantaneous power is recorded as Pk = ζ̇(tk )

N X ∂V (xi , ζt ) k

i=1

∂ζ

ρi (tk ),

(S43)

where ρi (tk ) is the i-th component of the actual Magnus4-evolved distribution at time tk . The cumulative work W(tk ) is obtained by integrating the discrete power series {Pk } using Simpson’s rule, W(t2m ) =

m−1 X

∆t P2j + 4P2j+1 + P2j+2 , 3 j=0

(S44)

which is accurate to O(∆t4 ). This is essential: a lower-order quadrature such as the midpoint or trapezoidal rule (O(∆t2 )) introduces a spurious residual in Wdiss that dominates the genuine tracking error and obscures the escorting equality. With the fourth-order Simpson integration, |Wdiss (t)| converges to much lower errors, consistent with the O(∆t4 ) Magnus integrator accuracy and the independently measured TVD, consistently in the range of ∼ 10−13 − 10−12 , confirming the trajectory-wise escorting equality at a level set by the time-discretisation of the dynamics rather than by any deficiency of the counterdiabatic construction. VI. A.

ADDITIONAL EXPERIMENTAL RESULTS

Free energy difference for the driven harmonic trap

For the harmonic potential V (x, t) = 12 κ(t)x2 , the canonical partition function at inverse temperature β is the Gaussian integral r Z ∞ 2 2π Z(κ) = e−βκx /2 dx = . (S45) βκ −∞

10 The Helmholtz free energy follows as F(κ) = −β

−1

βκ 1 ln . ln Z(κ) = 2β 2π

(S46)

The free energy difference between the final stiffness κf and the initial stiffness κi is therefore 1 κf ∆F = F(κf ) − F (κi ) = ln , 2β κi

(S47)

where the 2π factors cancel exactly. For the parameters used throughout this work, κi = 1, κf = 4, and β = 1, this gives the free-energy difference value to be: ∆F = 21 ln 4 = ln 2 ≈ 0.6931 .

(S48)

In the numerical simulations, the free energy is computed from the discrete partition function on the spatial grid {xi }N i=1 , Zd (κ) =

N X

e−βκxi /2 ,

(S49)

i=1

with Fd (κ) = −β −1 ln Zd (κ). For our grid of N = 80 points on q ∈ [−4, 4], the discrete free energy difference ∆Fd = Fd (κf ) − Fd (κi ) = 0.6931 agrees with the continuum result (S48) to ∼ 5 × 10−5 , confirming that the spatial discretisation introduces negligible error in the thermodynamic reference. TABLE S3: Comparison of terminal work and dissipation metrics across varying protocol durations τ . For the system driven by the harmonic potential Vh (x, ζt ), the metrics contrast the uncontrolled dynamics (bare) against the exact analytical counterdiabatic drive (analytic) and our numerical Liouvillian framework (LCD). All values are reported relative to the exact equilibrium free energy change ∆F = 0.6931. τ

W(τ )

W analytic (τ )

W LCD (τ )

∆F

Wdiss (τ )

analytic Wdiss (τ )

LCD Wdiss (τ )

0.001 0.003 0.010 0.030 0.100 0.300 1.000 3.000 10.00

1.4977 1.4958 1.4890 1.4702 1.4097 1.2727 1.0215 0.8274 0.7339

0.6931 0.6931 0.6931 0.6931 0.6931 0.6931 0.6931 0.6931 0.6931

0.8046 0.8027 0.7959 0.7771 0.7166 0.5796 0.3284 0.1343 0.0408

4.005 × 10−13 2.265 × 10−14 5.818 × 10−13 2.898 × 10−14 5.933 × 10−13 3.186 × 10−13 −5.635 × 10−13 5.285 × 10−14 −2.345 × 10−13

−6.867 × 10−12 −5.163 × 10−12 −5.701 × 10−12 −4.944 × 10−12 −5.149 × 10−12 −3.715 × 10−12 −2.268 × 10−12 −8.511 × 10−13 −2.666 × 10−13

Table S3 reports the terminal work W(τ ) and dissipated work Wdiss (τ ) = W(τ ) − ∆F for protocol durations spanning four decades, τ ∈ [10−3 , 10]. Under bare dynamics, the work W(τ ) exceeds ∆F at all finite τ , with Wdiss ranging from 0.80 (τ = 10−3 , strongly irreversible) to 0.04 (τ = 10, near quasi-static), consistent with the expected Wdiss ∼ 1/τ scaling in linear response. Both the analytical CD and the spectral LCD yield W(τ ) = ∆F to within ∼ O(10−12 ) across all protocol durations, with residuals attributable to the ODE solver tolerance (analytic) and the fourth-order Simpson work quadrature (LCD). This confirms that the escorting condition Wdiss (τ ) ≈ 0 holds irrespective of driving speed, eliminating the dissipative bias in the free energy estimate without statistical sampling.

VII.

STRESS TEST: THE QUARTIC COALESCENCE MODEL

The double-well and harmonic-trap systems studied above are both well-conditioned : the biorthogonal decomposition underlying the spectral LCD generator remains numerically stable throughout the protocol, and the resulting LCD driven density tracks the instantaneous equilibrium distribution to within machine precision. To probe the limits of this construction, we consider a deliberately adversarial model in which the slowest relaxation mode’s spectral gap closes exponentially over most of the protocol.

11 A.

Model and protocol

The quartic coalescence potential is V (x, ζ) = x4 − 16(1 − ζ) x2 ,

ζ(t) = min(t/τ, 1),

(S50)

annealed linearly from ζ = 0 to ζ = 1 over duration τ , at inverse temperature β = 1. The discrete generator L̂(t) is built on a grid of N = 80 points on x ∈ [−4.5, 4.5] via the same detailed-balance-preserving (Sasa–Tasaki-type) discretization [35] used throughout this work, Eq. (8) . √ At ζ = 0, Eq. (S50) is a symmetric double well with minima at x = ± 8 and a central barrier of height ∆Vb (ζ) = V (0, ζ) − V (xmin , ζ) = 64 (1 − ζ)2 ,

(S51)

so ∆Vb (0) = 64. As ζ → 1 the two minima coalesce (hence coalescence model ) into the single well V (x, 1) = x4 , with ∆Vb (1) = 0. Figure S3 shows the potential at six representative points spanning this transition, split explicitly at ζ = 0.5 (where ∆Vb = 16).

ζ = 0.00 τ

ζ = 0.20 τ

gap closed (metastable) 75 50 25 0 25

20 0 20 40

V(x, ζ)

gap closed (metastable)

100 50 0

ζ = 0.50 τ

ζ = 0.70 τ

gap open

ζ = 1.00 τ

gap open

150

200

100

150

200 100

0 2

gap open

300

100

50 4

ζ = 0.40 τ

gap closed (metastable) 150

x FIG. S3: The quartic coalescence potential: The potential V (x, ζ) depicted at six points along the anneal. For ζ < 0.5 (top row) the potential is a deep, symmetric double well; for ζ ≥ 0.5 (bottom row) the barrier has dropped below ∆Vb = 16 and continues to collapse smoothly into the single quartic well at ζ = 1.

Exponentially closing spectral gap

The slowest non-trivial eigenvalue λ1 (ζt ) of L̂(t) governs the inter-well (symmetric ↔ antisymmetric) relaxation rate of the double well. For an overdamped, symmetric double-well potential, Kramers/Eyring theory [44, 57] predicts this rate is exponentially suppressed by the barrier height, |λ1 (ζ)| ∼ C(ζ) exp − β ∆Vb (ζt ) = C(ζ) exp − 64(1 − ζ)2 , (S52) with C(ζ) a slowly-varying prefactor set by the local curvature at the well and barrier. We verified Eq. (S52) directly against the numerically diagonalized L̂(t): the exponential decay exp −64(1 − ζ)2 reproduces the computed |λ1 (ζ)|

12 to within a prefactor C(ζ) ∈ [2, 7] across the entire range ζ ∈ [0.5, 1] (Table S4). Below t/τ ≈ 0.33–0.4, ∆Vb ≳ 36 and exp(−∆Vb ) ≲ 10−16 falls below double-precision machine epsilon: λ1 (ζ) is no longer numerically resolvable and the dense eigensolver returns a value dominated by floating-point round-off (the noisy plateau at 1/|λ1 | ∼ 1012–13 visible in Fig. 8 for t/τ ≲ 0.33), rather than the true (exponentially smaller) physical gap. TABLE S4: Verification of the Kramers-type scaling. Numerical corroboration of Eq. (S52). The table monitors the first non-vanishing relaxation eigenvalue |λ1 (ζ)| against the classical barrier-height scaling prediction exp[−∆Vb (ζ)] across varying parameter configurations ζ. ζ 0.50 0.60 0.70 0.80 0.90 1.00

e−∆Vb

|λ1 (ζ)| (numeric)

∆Vb (ζ)

−7

7.30 × 10 1.89 × 10−4 1.23 × 10−2 2.00 × 10−1 1.04 2.72

16.11 10.31 5.80 2.58 0.64 0.00

ratio −7

1.01 × 10 3.33 × 10−5 3.03 × 10−3 7.60 × 10−2 5.25 × 10−1 1.00

7.2 5.7 4.1 2.6 2.0 2.7

Consequence for the counterdiabatic generator

The spectral LCD generator weights each mode by cn (t) = ζ̇t ℓn (∂ L̂/∂ζ)r0 /λn . Because λ1 (ζ) is exponentially small for ζ ≲ 0.5, c1 is exponentially large exactly where the counterdiabatic correction is most needed to suppress inter-well leakage – and the biorthogonal decomposition itself becomes numerically ill-conditioned in the same regime (cond(R) reaches ∼ 1012 near ζ = 0). This is a fundamentally different obstruction from the truncation-order question studied for the double well and harmonic trap (Section VII D): it is not that more modes are needed for convergence, but that the dominant mode’s own coefficient is unreliable. Table S5 confirms the practical consequence. Unlike the double well and harmonic trap, where the LCD generator suppresses Wdiss by 9–12 orders of magnitude relative to bare Fokker–Planck dynamics, here it achieves only a modest, τ -dependent improvement (∼33–1500×): Wdiss , rather than collapsing to the ∼ 10−12 floor seen elsewhere in this work. The exact equilibrium free-energy difference, ∆F = F (ζ=1) − F(ζ=0) = 62.9407...,

(S53)

computed directly from the discrete partition functions F(ζ) = −β −1 ln

e−βV (xi ,ζ) ,

(S54)

is in agreement with the reported reference value ∆F = 62.9407... [58]. TABLE S5: Quartic coalescence model. Work, dissipation at the terminal time step t/τ = 1.0. Comparison of bare vs. spectral LCD, for ∆F = 62.9407... (β = 1, N = 80; ∆t = 10−3 for τ ∈ [0.1, 0.5]. τ

W(τ )

W LCD (τ )

∆F

Wdiss (τ )

LCD Wdiss (τ )

improvement

0.10 0.20 0.50

86.48 78.87 71.48

63.66 63.03 62.95

62.94 62.94 62.94

23.54 15.92 8.54

0.71 0.09 0.006

33× 176× 1423×

We include this model precisely because it fails gracefully: it delineates the regime – deep, near-degenerate metastability, with a barrier large compared to kB T over a substantial fraction of the protocol – in which the spectral biorthogonal construction should not be expected to deliver near-machine-precision counterdiabatic driving, in contrast to the double-well and harmonic-trap results of Sections VII B–VII C. It is essential to distinguish the numerical breakdown of the spectral construction from a genuine divergence of counterdiabatic driving. For the symmetric coalescence potential, reflection symmetry P L̂P = L̂ [where (Pf )i = fN −1−i ] dictates that the stationary state r0 = πeq and the parametric perturbation ∂ L̂/∂ζ are even, whereas the

13 slowest relaxation mode r1 is odd. Consequently, the transition matrix element ℓT1 (∂ L̂/∂ζ)r0 vanishes identically by parity: escorting a symmetric πeq never requires transporting probability across the central barrier. Where the eigenvalues are well-separated [e.g., at ζ = 0.8 with |λ1 | ≈ 0.20], we confirm this parity protection numerically, observing an odd-to-even overlap suppression of |ℓT1 (∂ L̂/∂ζ)r0 |/|ℓT2 (∂ L̂/∂ζ)r0 | ∼ 10−11 . However, as the barrier height grows (ζ ≲ 0.4, ∆Vb ≳ 36), the Kramers gap |λ1 | ∼ e−β∆Vb drops below the double-precision machine floor (∼ 10−12 –10−13 ), rendering the eigenvectors severely ill-conditioned [cond(D) > 108 ]. In this regime, numerical eigensolvers fail to resolve near-degenerate eigenvalues into definite-parity eigenvectors; the protective zero-overlap is lost to floating-point round-off, and the assembled modal coefficient c1 = ℓT1 (∂ L̂/∂ζ)r0 /λ1 degenerates into an unphysical 0/0 numerical quotient that corrupts the spectral reconstruction of L̂CD [Eq. (34)]. By contrast, the numerics for the exact counterdiabatic operator remains completely bounded throughout the protocol. Because the closed-form rank-one representation L̂CD = (∂πeq /∂t)1T [Eq. (35)] bypasses eigenbasis inversion entirely, it never constructs c1 , maintains an O(1) operator norm, and tracks πeq (t) to a fidelity of TVD ∼ 10−6 (limited solely by spatial grid resolution and integrator tolerance). This confirms that the barrier-closure failure is an artifact of the spectral representation, not a physical limit on control. A genuine divergence of the classical counterdiabatic term—the true statistical analogue of a quantum level crossing—arises only when the target distribution itself must move rapidly across a closing gap, such as during an asymmetric barrier crossing or a first-order bifurcation of πeq (t). In those asymmetric geometries, ∥∂πeq /∂t∥ grows without bound. Thus, unlike the quantum adiabatic gauge potential, whose norm diverges universally whenever an energy gap closes, the classical Liouvillian control diverges only when probability must be actively transported through that closure.

Spectral truncation as imperfect escorting

The counterdiabatic correction L̂CD (t) [Eq. (34)] is a sum over all N − 1 non-zero relaxation modes. In practice, PM ℓT (t) (∂ L̂/∂t) r0 (t) (M ) the sum can be truncated to the M slowest modes, L̂CD (t) = − n=1 n rn (t) ℓT0 , yielding an imperfect λn (t) but systematically improvable escorting: the driven distribution no longer tracks πeq (t) exactly, but the tracking error decreases monotonically as M increases, interpolating between the uncontrolled bare dynamics (M = 0) and the full LCD (M = N − 1). The effectiveness of this truncation rests on the spectral structure of L̂(t). Since the relaxation rates |λn | grow with mode index (Figure S4), the coefficients cn ∼ |λn |−1 decay, and the dominant contribution comes from the modes whose relaxation timescales |λn |−1 are comparable to or longer than the driving timescale τ . The number of modes required for a given accuracy therefore depends on the protocol speed: for slow protocols fewer modes suffice to suppress the TVD by several orders of magnitude; for faster protocols modes with shorter relaxation timescales contribute more significantly, and M ∼ (N − 1)/2 is needed to achieve comparable improvement. TABLE S6: Truncated Liouvillian counterdiabatic driving. Maximum TVD, maximum DKL , and maximum absolute intermediate dissipated work max |Wdiss (t)| evaluated for several truncation thresholds M across the double-well and harmonic trap potentials at a fixed driving period τ = 0.1 (∆t = 10−6 ). Double Well Truncation (M ) 5 10 15 20 25 30 35 Full (M = 79)

max TVD −3

2.66 × 10 7.44 × 10−5 1.35 × 10−6 1.39 × 10−7 7.77 × 10−9 1.35 × 10−9 1.53 × 10−10 2.14 × 10−12

max DKL −5

3.33 × 10 2.71 × 10−8 1.10 × 10−11 1.16 × 10−13 8.54 × 10−16 4.53 × 10−16 4.54 × 10−16 3.88 × 10−16

Harmonic Trap max |Wdiss (t)| −4

1.34 × 10 1.97 × 10−7 1.96 × 10−10 5.42 × 10−12 2.86 × 10−12 2.85 × 10−12 2.84 × 10−12 2.85 × 10−12

Max TVD

Max DKL

max |Wdiss (t)|

−4

5.69 × 10−5 2.34 × 10−5 6.10 × 10−6 7.28 × 10−7 1.61 × 10−7 1.41 × 10−8 3.71 × 10−9 5.15 × 10−12

2.07 × 10 6.56 × 10−5 3.21 × 10−5 1.32 × 10−5 7.86 × 10−6 3.97 × 10−6 2.64 × 10−6 1.95 × 10−12

1.17 × 10 9.01 × 10−6 2.50 × 10−6 2.64 × 10−7 5.04 × 10−8 3.74 × 10−9 5.94 × 10−10 4.26 × 10−16

1.0

1.25

1/|λ1 (t)| 1/|λ2 (t)| 1/|λ3 (t)|

0.8

1.00 1/|λ1 (t)| 1/|λ2 (t)| 1/|λ3 (t)|

0.75 0.50

1/|λ4 (t)| 1/|λ5 (t)|

0.6

1/|λ4 (t)| 1/|λ5 (t)|

0.4 0.2

0.25 0.00 0.00

Inverse relaxation rate, 1/|λn (t)|

0.25

0.50

t/τ

0.75

1.00

(a) Inverse spectral gap magnitude of an overdamped FP equation driven by a double-well potential for τ = 0.1.

0.0

0.2

0.4

0.6

t/τ

0.8

1.0

(b) Inverse spectral gap magnitude of an overdamped FP equation driven by a harmonic-trap potential for τ = 0.1.

−1 FIG. S4: Temporal evolution of the inverse spectral gap. The inverse spectral gap ∆−1 n (t) ≡ |λ0 (t) − λn (t)| computed between the ground state (instantaneous equilibrium distribution) zero-mode (λ0 = 0) and the first few dominant biorthogonal eigenvalues λn (t) as a function of time t for (a) the double-well potential and (b) the harmonic trap landscape. The trajectories map the instantaneous relaxation timescales governing the non-symmetric transition rate matrix throughout the driving protocol.

TABLE S7: Spectral truncation convergence metrics for the time-varying DW potential. performance for fast (τ = 0.01) and slow (τ = 0.50) protocol durations. Errors drop systematically as the number of retained modes ntrunc approaches the full discrete grid size (N = 79). τ

Truncation (M )

Max TVD

Max DKL −4

Max |Wdiss (t)|

0.01 (dt = 10−6 )

5 10 15 20 25 30 35 Full (n = 79)

7.139 × 10 2.788 × 10−4 8.263 × 10−6 8.460 × 10−7 6.220 × 10−8 1.051 × 10−8 1.383 × 10−9 1.709 × 10−12

4.686 × 10 1.033 × 10−6 1.083 × 10−9 1.291 × 10−11 6.315 × 10−14 2.268 × 10−15 3.891 × 10−16 3.631 × 10−16

2.697 × 10−4 6.129 × 10−7 9.541 × 10−10 1.678 × 10−11 2.056 × 10−12 1.948 × 10−12 1.943 × 10−12 1.942 × 10−12

0.50 (Slow)

5 10 15 20 25 30 35 Full (n = 79)

6.605 × 10−4 1.594 × 10−5 2.784 × 10−7 2.809 × 10−8 1.566 × 10−9 2.709 × 10−10 3.059 × 10−11 1.077 × 10−12

1.594 × 10−6 1.069 × 10−9 4.111 × 10−13 4.860 × 10−15 4.409 × 10−16 3.959 × 10−16 4.004 × 10−16 4.125 × 10−16

3.490 × 10−5 4.347 × 10−8 4.035 × 10−11 1.344 × 10−12 1.021 × 10−12 1.018 × 10−12 1.019 × 10−12 1.020 × 10−12

−3

Efficient simulations in symmetric representation

Because the Fokker-Planck operator L̂(t) is non-symmetric due to the spatial discretization construction, computing its spectral decomposition directly requires dense, non-symmetric eigensolvers that scale poorly [O(N 3 )] and are prone to numerical instabilities, such as spurious complex eigenvalue pairs. Under conditions of detailed balance, this bottleneck can be systematically bypassed. By defining a time-dependent similarity transformation to the symmetric

15 TABLE S8: Spectral truncation convergence metrics for a harmonically confined Brownian particle. Trajectory-wise quantities converge rapidly with the truncation scale M , demonstrating highly predictable low-rank approximations. τ

Truncation (M )

Max TVD

Max DKL

−4

Max |Wdiss (t)|

0.01 (Fast)

5 10 15 20 25 30 35 Full (n = 79)

2.293 × 10 8.761 × 10−5 4.981 × 10−5 2.558 × 10−5 1.753 × 10−5 1.062 × 10−5 7.863 × 10−6 2.202 × 10−12

2.597 × 10 6.485 × 10−5 3.564 × 10−5 1.432 × 10−5 9.819 × 10−6 5.084 × 10−6 3.352 × 10−6 4.930 × 10−16

1.532 × 10−4 6.675 × 10−5 3.457 × 10−5 1.650 × 10−5 1.060 × 10−5 5.961 × 10−6 3.982 × 10−6 5.703 × 10−12

0.50 (Slow)

5 10 15 20 25 30 35 Full (n = 79)

1.516 × 10−4 3.211 × 10−5 1.255 × 10−5 4.063 × 10−6 2.165 × 10−6 9.791 × 10−7 6.216 × 10−7 1.299 × 10−12

6.234 × 10−6 4.117 × 10−8 4.537 × 10−9 4.555 × 10−10 1.372 × 10−10 3.107 × 10−11 1.335 × 10−11 4.558 × 10−16

5.381 × 10−6 2.684 × 10−7 5.349 × 10−8 8.275 × 10−9 3.123 × 10−9 9.346 × 10−10 4.667 × 10−10 3.308 × 10−12

gauge, −1/2 1/2 L̂sym (t) = πeq (t)L̂(t)πeq (t),

the generator is mapped to a real symmetric tridiagonal matrix. This structure allows the use of highly optimized tridiagonal algorithms (e.g., standard bisection and inverse iteration), which drastically reduce the computational complexity of eigenbasis reconstruction. To rigorously validate this acceleration scheme, we assess both its numerical fidelity and performance across a wide parameter space. We evaluate the framework at a spatial grid resolution of N = 80 across uniform temporal snapshots t/τ ∈ {0, 0.25, 0.50, 0.75, 1.00} for driving protocols varying by three orders of magnitude (τ ∈ {0.01, 0.1, 1.0}). Benchmarked execution times are reported as mean values over 100 iterations after a 10-cycle warm-up phase. The results for the double-well and driven harmonic landscapes are detailed in Table S9 and Table S10, respectively. Across all protocols, the symmetry of the gauge-transformed operator is maintained near the double-precision limit, with the non-symmetry metric strictly bounded by max ∥L̂sym − L̂Tsym ∥ ∼ O(10−14 − 10−13 ). Spectral preservation under Eq. (VII E) is corroborated by checking the maximum absolute discrepancy between the eigenvalues of the original non-symmetric operator and its symmetric-gauge counterpart. The numerical difference remains exceptionally small, varying from O(10−12 ) in the smoothly behaving harmonic potential to O(10−8 ) in the highly driven double-well potential, where sharp boundary layers can transiently enhance numerical gradients. Crucially, this symmetric-gauge formulation yields a dramatic performance improvement. For the double-well landscape, transitioning to the tridiagonal solver cuts the average diagonalisation time from 2.01 ms to 0.17 ms, representing a consistent ∼ 12× speedup. For the harmonic potential, the speedup is even more pronounced, consistently exceeding 15×. Since counterdiabatic driving requires evaluating the Liouvillian eigenbasis at every temporal step along the driving path, this order-of-magnitude reduction in execution time is critical for making continuous non-equilibrium tracking computationally feasible in complex systems.

Scalability: Sparsity of the Liouville Operator and Extension to Higher Dimensions

A natural question is whether our biorthogonalization framework to compute the L̂CD extends beyond one spatial dimension. We address this by analysing the structure of the discretised Liouville operator L̂(t) and showing that its sparsity, combined with the spectral truncation result of Section VII D, makes the construction tractable in two and three dimensions. In d spatial dimensions, the state space is discretized on a regular grid of Nx points per dimension, giving a total of N = Nxd grid points. The probability vector ρ(t) ∈ RN is obtained by flattening the d-dimensional grid into a single vector, and the Liouville operator becomes an N × N rate matrix. Table S11 summarises the matrix size, memory requirement, and full-diagonalisation cost for a spatial discretization of say Nx = 50 across dimensions.

Wall-clock time per call (ms)

100

10 1 50

100

Spatial discretization, N

Dense eig(L) -- Double well eigh_tridiagonal(Lsym ) -- Double well

110

120

Dense eig(L) -- Harmonic trap eigh_tridiagonal(Lsym ) -- Harmonic trap

FIG. S5: Diagonalisation wall-clock time versus spatial discretisation N . Dense non-symmetric eigensolver −1/2 1/2 (applied to L̂(t)) compared with the tridiagonal symmetric eigensolver applied to L̂symm = πeq L̂ πeq . The −8 −13 symmetrised form yields a consistent 12–15× speedup, with eigenvalues agreeing to 10 –10 . TABLE S9: Symmetric-gauge verification for double-well potential. Symmetry residual, eigenvalue agreement, and diagonalisation wall-clock time (averaged over 100 runs, N = 80) at five checkpoints along the protocol for three durations τ . The tridiagonal solver achieves a consistent 11–12× speedup. τ 0.01

0.10

1.00

t/τ 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00

max ∥L̂sym − L̂Tsym ∥ 1.42 × 10−13 1.42 × 10−13 1.42 × 10−13 1.42 × 10−13 1.14 × 10−13 1.42 × 10−13 1.42 × 10−13 1.42 × 10−13 1.14 × 10−13 1.14 × 10−13 1.42 × 10−13 1.42 × 10−13 1.42 × 10−13 1.42 × 10−13 1.14 × 10−13

max ∥∆λ∥ 5.60 × 10−8 4.25 × 10−8 1.55 × 10−8 7.02 × 10−9 5.74 × 10−9 5.60 × 10−8 4.25 × 10−8 1.55 × 10−8 1.20 × 10−8 5.74 × 10−9 5.60 × 10−8 4.25 × 10−8 1.55 × 10−8 7.02 × 10−9 5.74 × 10−9

tdense (ms) 2.016 1.906 2.014 2.073 1.953 2.005 1.904 1.991 2.056 1.941 2.005 1.894 1.991 2.059 1.951

ttri (ms) 0.166 0.167 0.165 0.167 0.167 0.166 0.166 0.166 0.167 0.168 0.166 0.166 0.165 0.167 0.167

Speedup 12.1× 11.4× 12.2× 12.4× 11.7× 12.1× 11.5× 12.0× 12.3× 11.6× 12.1× 11.4× 12.0× 12.3× 11.7×

Full diagonalisation is therefore already intractable in 3D for spatial discretization of say, Nx = 50. The key structural property that rescues the situation is that L̂(t) is extremely sparse, which can be seen from Figures (S1, S2). In d dimensions with nearest-neighbour coupling (the Sasa–Tasaki discretisation [42]), each grid point i = (i1 , . . . , id ) is connected only to its 2d nearest neighbours i ± ek (k = 1, . . . , d), where ek is the unit vector in direction k. The column of L̂ corresponding to site i therefore has exactly 2d off-diagonal non-zero entries (the transition rates to/from neighbours) plus one diagonal entry (the negative total outgoing rate). The total number of non-zero entries is thus: nnz(L̂) = (2d + 1) N = O(d N ),

(S55)

growing linearly in the number of states, not quadratically. The sparsity — fraction of non-zero entries — is (2d+1)/N ,

17 TABLE S10: Symmetric-gauge verification: harmonic trap. Same layout as Table S9. The tridiagonal solver achieves a 14–16× speedup, with eigenvalue agreement ranging from 10−12 to 10−8 along the protocol. τ 0.01

0.10

1.00

t/τ 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00

max ∥L̂sym − L̂Tsym ∥ 4.26 × 10−14 4.26 × 10−14 4.26 × 10−14 4.26 × 10−14 5.68 × 10−14 4.26 × 10−14 4.26 × 10−14 4.26 × 10−14 2.84 × 10−14 5.68 × 10−14 4.26 × 10−14 4.26 × 10−14 4.26 × 10−14 4.26 × 10−14 5.68 × 10−14

max ∥∆λ∥ 2.56 × 10−12 2.27 × 10−12 1.81 × 10−10 2.56 × 10−8 6.56 × 10−8 2.56 × 10−12 2.27 × 10−12 1.81 × 10−10 1.71 × 10−8 6.56 × 10−8 2.56 × 10−12 2.27 × 10−12 1.81 × 10−10 2.56 × 10−8 6.56 × 10−8

tdense (ms) 2.444 2.392 2.275 2.682 2.568 2.450 2.395 2.279 2.649 2.571 2.451 2.390 2.267 2.670 2.563

ttri (ms) 0.163 0.165 0.166 0.169 0.170 0.163 0.165 0.166 0.169 0.170 0.163 0.165 0.166 0.169 0.170

Speedup 15.0× 14.5× 13.7× 15.8× 15.1× 15.1× 14.5× 13.7× 15.7× 15.1× 15.0× 14.5× 13.7× 15.8× 15.1×

TABLE S11: Scaling of the Liouville operator L̂(t) with spatial dimension d for Nx = 50 grid points per dimension. Memory assumes dense storage (8 bytes [FLOAT64] per entry). Diagonalisation cost assumes O(N 3 ) complexity. d 1 2 3

N = Nxd 50 2,500 125,000

Matrix size 50 × 50 2500 × 2500 1.25×105 × 1.25×105

Memory < 1 MB 50 MB 125 GB

Diag. cost ∼ 105 ops ∼ 1010 ops ∼ 1015 ops

which decreases as N −1 for fixed d: sparsity =

(2d + 1) N 2d + 1 N →∞ = −−−−→ 0. N2 N

(S56)

TABLE S12: Non-zeros per column and sparsity of L̂(t) for Nx = 50 per dimension. d 1 2 3

nnz/column 3 5 7

Total nnz 150 12,500 875,000

Sparsity 6.0% 0.2% 0.006%

The sparsity arises directly from the local structure of the Fokker–Planck equation: probability can only flow between adjacent grid points, so the generator L̂(t) is a graph Laplacian on the d-dimensional lattice, with non-zeros only at positions corresponding to neighboring pairs. In memory, L̂(t) is therefore never stored as a dense N × N array but as a sparse matrix with O(dN ) entries — reducing the 3D memory requirement from 125 GB (dense) to ∼ 7 MB (sparse, N = 125,000, d = 3). As already noted in Section V B for the equilibrium problem the rank-one closed form is O(N ), and needs none of this; the sparse/truncated construction is relevant to the diagnostic analysis and to the NESS setting where r0 must be computed numerically.

Record · ID 405675 · SHA-256 9bb90cd7133f807b
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.