Global thermodynamic manifold for conservative control of stochastic systems

In the context of conservative driving, ‘full control’ refers to the ability to manipulate the entire potential-energy landscape. For a system with $N$ discrete states, the state energies can be assembled into a vector $\bm{V}=(V^{0},V^{1},\dots,V^{N-1})$ with $V^{\nu}=V(x_{\nu})$. For a continuous system, one may imagine a continuum of control parameters $V(x)$ parameterized by $x$. The forces conjugate to the state energies are indicator functions

where the functional derivative accounts for continuous $x$. The average deviation of the conjugate forces $\omega$ from their equilibrium values is then precisely the (negative) deviation of the probability distribution from equilibrium:

Analogous to (18), we define the full-control friction tensor $\zeta(V)$—which depends on the entire potential-energy landscape $V$—by its integral kernel

such that

Interpreted geometrically, the full-control friction tensor is a metric on a (perhaps infinite-dimensional) manifold that we call the global thermodynamic manifold $\mathcal{M}$, with each possible landscape identified with a point on $\mathcal{M}$. The squared length of a tangent vector on the global manifold corresponds to the linear-response excess power:

where $||\cdot||_{\zeta}$ denotes a norm with respect to the metric $\zeta$.

We will see later that the metric $\zeta$ has a non-empty null space spanned by $1$ (respectively the constant function $1$ for continuous $\Omega$ or the constant vector $\bm{1}$ for discrete $\Omega$). We should expect this intuitively from the linear-response relation implied by Eq. (25):

Uniformly raising or lowering the energy landscape does not affect the system’s dynamics, so this action does not alter the system distribution: $\delta p=0$ for $\mathrm{d}V/\mathrm{d}t\propto 1$. An equivalent argument may be made from conservation of probability ³³3Taking the inner product of both sides of (27) with the constant function $1$ and applying $(1,\delta p)=0$ leads to the same conclusion. Since $\zeta$ is symmetric, in this context the arbitrariness of the potential-energy zero-point and conservation of probability are equivalent. This highlights a duality between energy and probability, manifest in the geometry: the linear-response relation (27) states that $\delta p$ is the vector dual to $\mathrm{d}V/\mathrm{d}t$. This is also reflected in the conjugate forces to the state energies (22)..

Since $\text{det}\,\zeta=0$, the metric is said to be degenerate and the full-control manifold $\mathcal{M}$ is thus a non-Riemannian manifold. This degeneracy does not meaningfully impede the usual geometric analysis, however, as the null vector (in the coordinates of $V$) is the same everywhere on the manifold. The degeneracy may be managed by projecting points $V\in\mathcal{M}$ onto some hypersurface on which the metric is non-degenerate, e.g., by fixing the potential zero-point, setting $V(y)=0$ for some $y\in\Omega$. We address the construction of Riemannian full-control submanifolds for discrete-state systems in Sec. IV.3.

In practical settings, the number of control parameters is significantly smaller than the number of system states. For instance, early studies of fluctuation theorems involved the manipulation of single RNA molecules with the molecular extension as the only controllable parameter [38, 39]. The success of these demonstrations and other nonequilibrium measurements of equilibrium free-energy landscapes relies on executing experimental protocols with low excess work [4]. Partial (or parametric) control is thus an area of clear practical importance.

For conservative control with a set of $m$ control parameters $\left\{u^{0},\dots,u^{i},\dots,u^{m-1}\right\}$, protocols are constrained by $V(t)=V(\bm{u}(t))$ to a lower-dimensional control-parameter space. Geometrically, this space is a submanifold $\tilde{\mathcal{M}}$ of the global manifold $\mathcal{M}$.

Since the control parameters represent a constrained control over the potential energy $V$, there is a map (unique up to specification of the energy zero point) between points $\bm{u}$ on $\tilde{\mathcal{M}}$ and points $V(\bm{u})$ on $\mathcal{M}$. In differential geometry, such a map is an immersion:

The partial-control manifold is called a submanifold of the global manifold with immersion $\mathcal{F}$.

Associated with an immersion is its pushforward (or differential) $\mathcal{F}_{*}:\text{T}_{\bm{u}}\tilde{\mathcal{M}}\to\text{T}_{V(\bm{u})}\mathcal{M}$, mapping vectors on $\tilde{\mathcal{M}}$ (i.e., elements of the tangent space $\text{T}_{\bm{u}}\tilde{\mathcal{M}}$) to vectors on $\mathcal{M}$. In this setting, the matrix elements of the pushforward are simply obtained by the chain rule, and are identified with the forces $f_{i}(x)$ conjugate to respective parameters $u^{i}$ when the system is in state $x$:

The metric $\tilde{\zeta}(\bm{u})$ on $\tilde{\mathcal{M}}$, inherited from $\zeta(V)$, is then

(Substituting (24) into this form and applying (22) shows it to indeed be equivalent to the metric defined in Eq. (18).) This is constructed to ensure that geometric invariants (e.g., the lengths of tangent vectors) are equal when computed on $\mathcal{M}$ and on $\tilde{\mathcal{M}}$. In particular, since the squared length of a vector on a thermodynamic manifold corresponds to the linear-response excess power, the preservation of geometric invariants under coordinate transformations means that one computes the same excess power from $\tilde{\zeta}$ and $\bm{u}(t)$ as from $\zeta$ and $V(t)$. Figure 2 summarizes this picture.

For a discrete-state system with $|\Omega|=n$, full control is possible with a set of $n-1$ parameters due to the invariance of the dynamics under a uniform shift of the state energies (captured in the degeneracy of $\zeta$). We define a control-parameter set $\bm{u}\in\mathbb{R}^{n-1}$ as complete if the potential energy is fully determined (up to a global additive constant) by $\bm{u}$.

Define the Jacobian $\mathbb{J}=\mathcal{F}_{*}^{\mathsf{T}}$ of the transformation from the control parameters $\bm{u}$ to the state energies $\bm{V}$, i.e., $\left[\mathbb{J}\right]^{\phantom{i}\nu}_{i}=-f_{i}(x_{\nu})=\partial V^{\nu}/\partial u^{i}$. Then for probability vector $\bm{p}(t)$,

In App. A, we show that the completeness condition is equivalent to the invertibility of (31), i.e., the distribution $\bm{p}(t)$ is completely determined by the average conjugate forces $\left\langle\bm{f}\right\rangle_{\bm{p}(t)}$. The inverse of (31) is

for right inverse $\mathbb{J}^{\mathsf{R}}=\mathbb{J}^{\mathsf{T}}\left(\mathbb{J}\mathbb{J}^{\mathsf{T}}\right)^{-1}$. The time evolution of the average deviation $\left\langle\delta\bm{f}(t)\right\rangle$ of the conjugate force to a complete set of control parameters is thus self-contained in the sense that

where

is the generator of the dynamics of $\left\langle\bm{f}(t)\right\rangle$.

There is an equivalent geometric condition on the completeness of a control-parameter set, illustrated for a three-state system in Fig. 3. Let $\Delta_{\bm{f}}$ be the convex hull of the $n$ points $\left\{\bm{f}(x_{i})\right\}$, i.e., the set of all conjugate-force averages $\left\langle\bm{f}\right\rangle$. For a complete set of control parameters, $\Delta_{f}$ forms an $(n-1)$-dimensional simplex in $\mathbb{R}^{n-1}$. Since the space of probability vectors $\bm{p}$ for an $n$-state system defines an $(n-1)$-dimensional simplex (the probability simplex $\Delta_{\bm{p}}=\left\{\bm{p}\in\mathbb{R}^{n}\ :\ \bm{1}^{\mathsf{T}}\bm{p}=1\right\}$), there exists a bijective transformation between these spaces. Put simply, the space of conjugate-force averages $\left\langle\bm{f}\right\rangle$ for a complete set of control parameters is the same “size” as the space of probability vectors.

Concretely, consider an $N$-spin system, for which there are $n=2^{N}$ states. We show in Sec. VI.3 that the set of all $k$-spin interactions for $k=1,\dots,N$ from the cluster expansion is a complete control set that spans a $(2^{N}-1)$-dimensional Riemannian full-control submanifold.

We begin with a brief review of the Drazin inverse, a generalized inverse for linear operators which plays a central role in the derivations that follow. The Drazin inverse of the rate matrix has also appeared in previous works on optimal control [40, 41, 42, 43, 23] and has applications in first-passage problems [44, 45].

Intuitively, the need for an appropriately defined inverse of the generator arises due to the time integral in the definition (44) of the full-control friction tensor. In the same way as the generator $\mathcal{G}_{t}$ plays the role of a time derivative on the space of probability distributions, its Drazin inverse $\mathcal{G}^{\mathcal{D}}_{t}$ plays the role of a time integral (over $[0,\infty)$) on the same space.

Let $\mathcal{A}$ be a linear operator on some Hilbert space. The Drazin inverse of $\mathcal{A}$ is the unique operator $\mathcal{A}^{\mathcal{D}}$ satisfying

for some $k\in\mathbb{Z}$ called the index of $\mathcal{A}$. For operators on Hilbert spaces of finite dimension (i.e., matrices), $k$ is the smallest non-negative integer such that $\text{rank}\,\mathcal{A}^{k}=\text{rank}\,\mathcal{A}^{\ell}$ for all $\ell>k$ [46]. We consider only the case of $k=1$ (such as for irreducible rate matrices [47]). In this case, the Drazin inverse of $\mathcal{A}$ may be constructed from

for projection matrix $P_{0}$ onto the null space $\mathcal{N}(\mathcal{A})$ of $\mathcal{A}$.

From properties (35a-c), one can show that

for identity operator $I$. That is, $\mathcal{A}\mathcal{A}^{\mathcal{D}}$ projects vectors in the Hilbert space onto the complement of $\mathcal{N}(\mathcal{A})$.

Concretely, suppose that $\mathbb{W}$ is an irreducible rate matrix with equilibrium distribution $\bm{\pi}$. Then the product $\mathbb{W}\mathbb{W}^{\mathcal{D}}$ (and $\mathbb{W}^{\mathcal{D}}\mathbb{W}$ by the commutativity property (35c)) maps a probability distribution to its deviation from equilibrium, i.e., $\mathbb{W}\mathbb{W}^{\mathcal{D}}\bm{p}=(I-\Pi)\bm{p}=\delta\bm{p}$ with $\Pi=\bm{\pi}\bm{1}^{\mathsf{T}}=\left[\bm{\pi}\,\cdots\,\bm{\pi}\right]$. As a corollary,

so that $\mathbb{W}^{\mathcal{D}}$ is a proper inverse of $\mathbb{W}$ when the domain is restricted to vectors in the hyperplane $\left\{\bm{v}\ :\ \bm{1}^{T}\bm{v}=0\right\}$ – including deviations $\delta\bm{p}$ from equilibrium and derivatives $\mathrm{d}\bm{p}/\mathrm{d}t$ of probability distributions.

From properties (35b,c), it follows that ⁴⁴4Property (35b) can be re-written as $(\mathcal{A}^{\mathcal{D}})^{2}\mathcal{A}=\mathcal{A}^{\mathcal{D}}$ by the commutative property (35c). Thus $\mathcal{A}x=0\implies\mathcal{A}^{\mathcal{D}}x=0$. An identical argument holds for right-multiplication.

The Drazin inverse may thus be interpreted as “inverting the invertible part” of a matrix [40], leaving the null space of $\mathbb{W}^{\mathcal{D}}$ coincident with the null space of $\mathbb{W}$.

Consider some continuous- or discrete-state Markov process with generator $\mathcal{G}$ satisfying the conditions established in Sec. II. We would like to manipulate the expression (24) for the full-control friction tensor into a more manageable form. We take advantage of the fact that the integrand $\Sigma_{\omega}(t)=\left\langle\delta\omega(t)\delta\omega(0)\right\rangle_{\text{eq}}$ is calculated for a fixed energy landscape so that $\mathrm{d}\pi/\mathrm{d}t=0$ and thus $\mathrm{d}\Sigma_{\omega}/\mathrm{d}t=\mathcal{G}\Sigma_{\omega}$. Then from (38), the full-control friction tensor is

for $\Sigma_{\omega}^{0}\equiv\Sigma_{\omega}(0)$. The final line follows from the fact that $\lim_{t\to\infty}\Sigma_{\omega}(t)=0$ for ergodic dynamics.

Equation (40c) is a main result: the full-control friction tensor is decomposable into a product of the Drazin inverse of the generator and the equal-time equilibrium covariance matrix of the indicator functions $\omega$. Roughly, this separates the friction tensor into a static part $\Sigma_{\omega}^{0}$ and a dynamic part $\beta\mathcal{G}^{\mathcal{D}}$.

This expression simplifies further. Appendix B shows that the equal-time equilibrium covariance matrix of the indicator functions $\omega$ satisfies

Since $\pi$ is in the null-space of $\mathcal{G}^{\mathcal{D}}$, we obtain the decomposition

For discrete systems, this reduces the infinite-time integral (24) to a matrix equation

Together with the formula (36) for the Drazin inverse of a matrix, this expression allows for machine-precision numerical calculation of the full-control friction tensor $\zeta(V)$ – and thus any partial-control tensor $\tilde{\zeta}(\bm{u})$ via the analysis of Sec. IV.2.

In App. C, we show that for a complete control set $\bm{u}$, the friction tensor is

for generator $\mathbb{W}_{\bm{f}}$ of the dynamics of $\left\langle\bm{f}(t)\right\rangle$ as defined in Eq. (34) and equal-time correlation matrix $\Sigma_{\bm{f}}^{0}=\left\langle\delta\bm{f}\delta\bm{f}^{\mathsf{T}}\right\rangle_{\bm{\pi}(\bm{u})}$.

Lastly, since a protocol $V(t)$ can be completely characterized up to an (irrelevant) additive constant by the evolution of the equilibrium distribution $\pi(t)$, we use the operator decomposition to derive the relation (App. E)

The linear-response relation (27) is then equivalent to the approximation

This expression was derived previously using an operator expansion [41].

As noted in [11], the friction tensor $\tilde{\zeta}$ (for any set of control parameters) may be written as

the Hadamard product of the integral relaxation time matrix

and the equilibrium Fisher-information matrix

When the conjugate forces $\bm{f}$ are the indicator functions $\bm{\omega}$, equating (47) with the operator decomposition (43) gives a simple correspondence between the Drazin inverse of the rate matrix and the integral relaxation time:

Seen another way, defining the transition probability matrix with elements $p_{\mu\nu}(t)=\text{Pr}\left[\mu,t\ |\ \nu,0\right]$, the Drazin inverse of the rate matrix is [45]

Thus, much like $\mathcal{T}$, the Drazin inverse of the rate matrix describes the timescale of system relaxation. The magnitude of $\mathbb{W}^{\mathcal{D}}_{\mu\nu}$ characterizes how long it takes for $p_{\mu\nu}(t)$, the probability $p_{\mu}(t)$ conditioned on initialization in state $\nu$, to reach the equilibrium probability $\pi_{\mu}$.

When a generator $\mathcal{G}$ admits the spectral decomposition (9), the integral kernel of the Drazin inverse is

for relaxation times

See App. D for a concise proof. Substituting this expression into the decomposition (42) and using the fact that $\psi_{\alpha}(x)=\phi_{\alpha}(x)\pi(x)$ for systems obeying detailed balance [49], we obtain a spectral decomposition of the full-control friction tensor:

The degeneracy of $\zeta$ as inferred from Eq. (27) is manifest in this form: from the biorthogonality of the eigenfunctions, we have $\phi_{0}\zeta=\zeta\phi_{0}=0$, with $\phi_{0}=1$.

This spectral decomposition can be used to express the linear-response excess power (25) as a sum of independent contributions associated with the fixed-parameter relaxation modes of the system:

Under the change of variables $s\equiv t/\tau$ for total protocol time $\tau$, the linear-response excess work is

Equations (54)-(56) constitute another main result. The spectral decomposition of the excess work provides insight into the dynamical origin of dissipation: roughly, each term captures the degree to which the protocol traverses a given mode, scaled by that mode’s relaxation time. As a design principle, dissipation may be reduced by driving a system such that $\mathrm{d}V/\mathrm{d}t$ is maximally orthogonal to modes with long relaxation times.

A similar expression holds for partial-control submanifolds. From the transformation law for metrics,

for conjugate relaxation modes

The zeroth conjugate relaxation mode $\tilde{\bm{\psi}}_{0}$ is the equilibrium average $-\left\langle\bm{f}\right\rangle_{\pi}$ of the conjugate-force vector, and for $\alpha\geq 1$, $\tilde{\bm{\psi}}_{\alpha}$ captures the fixed-parameter relaxation of $-\left\langle\delta\bm{f}(t)\right\rangle$ corresponding to system relaxation along mode $\alpha$. The linear-response excess power in terms of a limited set of control parameters is then

Equation (57) is particularly powerful for continuous systems for which the spectrum of $\mathcal{L}$ is analytically known. While the full-control spectral decomposition (54) is a manifestly infinite sum for continuous state spaces, the partial-control spectral decomposition (57) may be a finite sum if the conjugate forces $f_{i}(x)$ are orthogonal to all but a finite set of modes $\psi_{\alpha}(x)$; Section (VI.2) demonstrates a concrete example of this.

We begin with a very simple model system: two states with arbitrary transition rates $w_{10}$ and $w_{01}$ (Fig. 4a, left). From Eq. (36), the Drazin inverse of the rate matrix is

From the decomposition (43) and instantaneous equilibrium distribution $\left(\pi_{0}(t),\pi_{1}(t)\right)^{\mathsf{T}}=\left(\frac{w_{01}}{w_{01}+w_{10}},\frac{w_{10}}{w_{01}+w_{10}}\right)^{\mathsf{T}}$ the full-control friction tensor is

for relaxation time $\tau_{\text{rel}}(t)=(w_{01}(t)+w_{10}(t))^{-1}$. Note that the two-state system has only one non-zero eigenvalue (and thus only one relaxation time $\tau_{\text{rel}}\equiv\tau_{1}$). It may be readily verified that the spectral decomposition $\zeta(\bm{V})=\beta\tau_{\text{rel}}\bm{\psi}_{1}\bm{\psi}_{1}^{\mathsf{T}}$ yields the same result.

The control parameters $\bm{V}$ may be freely defined subject to

Thus the energy difference $\Delta\equiv V^{1}-V^{0}$ is a natural parameter to consider. The linear-response excess power is

and any LR minimum-work protocol $\Delta(t)$ satisfies

A physical instantiation of a two-state system is a single-level quantum dot interacting with a metallic reservoir (Fig. 4a, right). In the wide-band approximation (i.e., the metallic reservoir’s density of states is assumed constant [50]), the energy $\varepsilon$ of the quantum dot and the chemical potential $\mu$ of the reservoir set the respective tunneling rates $w_{10}=2R_{0}[1+\text{e}^{\beta(\mu-\varepsilon)}]^{-1}$ into the reservoir and $w_{01}=2R_{0}[1+\text{e}^{\beta(\varepsilon-\mu)}]^{-1}$ out of the reservoir, where $R_{0}$ is the tunneling rate when $\mu=\varepsilon$ [5]. Then the friction tensor (for one control parameter, the friction coefficient) is

for $\Delta=\varepsilon-\mu$.

Consider the overdamped diffusion $X_{t}$ of a particle in a one-dimensional harmonic trap $V(x,t)=\frac{1}{2}\kappa(t)(x-x_{0}(t))^{2}$ with stiffness $\kappa(t)$ and trap center $x_{0}(t)$ (Fig. 4b). The Fokker-Planck operator governing the evolution of the particle’s position distribution $p(x,t)$ is