# Inferring entropy production from time-dependent moments

Prashant Singh and Karel Proesmans  Niels Bohr International Academy, Niels
Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen,
Denmark<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Measuring entropy production of a system directly from the experimental data
is highly desirable since it gives a quantifiable measure of the time-
irreversibility for non-equilibrium systems and can be used as a cost function
to optimize the performance of the system. Although numerous methods are
available to infer the entropy production of stationary systems, there are
only a limited number of methods that have been proposed for time-dependent
systems and, to the best of our knowledge, none of these methods have been
applied to experimental systems. Herein, we develop a general non-invasive
methodology to infer a lower bound on the mean total entropy production for
arbitrary time-dependent continuous-state Markov systems in terms of the
moments of the underlying state variables. The method gives surprisingly
accurate estimates for the entropy production, both for theoretical toy models
and for experimental bit erasure, even with a very limited amount of
experimental data.

##  1 Introduction

In the last two decades, the framework of stochastic thermodynamics has
enabled us to give thermodynamic interpretation to mesoscopic systems that are
driven arbitrarily far from equilibrium [1]. Using the assumption of local
detailed balance and a time scale separation between the system and its
environment, one can define thermodynamic quantities such as heat, work and
entropy production for general Markov systems [2, 3]. General results such as
fluctuation theorems and thermodynamic uncertainty relations have been
proposed and verified experimentally [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. A
fundamental quantity in stochastic thermodynamics is the entropy production
[1, 2]. Indeed, minimising entropy production is an important problem in
several applications, such as heat engines [15], information processing [16]
and biological systems [17, 18]. Furthermore, it provides information and
places bounds on the dynamics of the system in the form of the thermodynamic
speed limits [20, 22, 23, 19, 25, 21, 24, 26, 27] and the dissipation-time
uncertainty relation [28, 29, 30, 31]. Meanwhile measuring entropy production
experimentally is difficult since one needs extensive knowledge about the
details of the system, requiring a large amount of data.

Over the last couple of years, several methods have been developed to infer
entropy production and other thermodynamic quantities in a model-free way or
with limited information about the underlying system [32]. These techniques
involve estimating entropy production using some experimentally accessible
quantities. This includes methods based on thermodynamic bounds, such as the
thermodynamic uncertainty relation [33, 34, 35], inference schemes based on
the measured waiting-time distribution [36, 37, 38, 39, 40], and methods based
on ignoring higher-order interactions [41, 42]. However, most of these methods
focus only on steady-state systems and do not work for time-dependent
processes [43, 44, 45, 46, 47]. Nonetheless, time-dependent dynamics plays a
crucial role in many systems. One can for example think about nano-electronic
processes, such as bit erasure or AC-driven circuits, but also in biological
systems, where oscillations in, e.g., metabolic rates or transcription
factors, effectively lead to time-dependent driving [50, 51, 48, 49]. There
are few methods to infer entropy production for time-dependent systems [11,
52, 12, 53, 54, 55], and it is often challenging to apply them to experimental
data. For instance, thermodynamic bound-based methods require the measurement
of quantities such as response functions [11, 52] or time-inverted dynamics
[12], which require an invasive treatment of the system. Another interesting
method combines machine learning with a variational approach [54]. This method
does however require a large amount of experimental data and a high temporal
resolution. This might also be the reason why, to the best of our knowledge,
no experimental studies exist on the inference of entropy production under
time-dependent driving.

This paper aims to address these issues by deriving a general lower bound on
the total entropy production for continuous-state time-dependent Markov
systems in terms of the time-dependent moments of the underlying state
variables. The method gives an analytic lower bound for the entropy production
if one only focuses on the first two moments and reduces to a numerical scheme
if higher moments are taken into account. The scheme also gives the optimal
force field that minimizes the entropy production for the given time-dependent
moments. After presenting the theoretical calculations and applications to
analytically solvable models, we will test our method in an experimental
scenario involving the erasure of computational bits [50]. The method works
surprisingly well for all examples, even when only the first two moments are
taken into account and with a very limited amount of data ($<100$
trajectories).

In the following, we first recall some known results from stochastic
thermodynamics in Section 2, which will then be used to derive bounds first
for one-dimensional systems in Section 3 and subsequently in general
$d$-dimensional systems in Section 4. Subsequently, we apply our method to the
theoretical toy examples in Section 5 and to the experimental trajectories
associated with erasure operation of computational bits in Section 6. Finally,
we conclude in Section 7.

##  2 Stochastic thermodynamics

We consider a $d$-dimensional stochastic system whose state at time $t$ is
denoted by a variable $\mathbb{x}(t)=\\{x_{1}(t),x_{2}(t),....,x_{d}(t)\\}$.
The system experiences a time-dependent force field
$\mathbb{F}\left(\mathbb{x}(t),t\right)$ and is in contact with a thermal
reservoir at constant temperature $T$. The stochastic variable $\mathbb{x}(t)$
evolves according to an overdamped Langevin equation

$\displaystyle\frac{d\mathbb{x}}{dt}=\left(k_{B}T\right)^{-1}\mathbb{D}~{}\mathbb{F}\left(\mathbb{x}(t),t\right)+\boldsymbol{\zeta}(t),$
(1)

$k_{B}$ is the Boltzmann constant and $\boldsymbol{\zeta}(t)$ represents the
Gaussian white noise with $\left\langle\boldsymbol{\zeta}(t)\right\rangle=0$
and
$\left\langle\boldsymbol{\zeta}(t)\boldsymbol{\zeta}(t^{\prime})\right\rangle=2\mathbb{D}~{}\delta(t-t^{\prime})$,
with $\mathbb{D}$ being the diffusion matrix. One can also write the Fokker-
Planck equation associated with the probability distribution
$P\left(\mathbb{x},t\right)$ as

$\displaystyle\frac{\partial P\left(\mathbb{x},t\right)}{\partial
t}=-\mathbb{\nabla}\cdot\left[\mathbb{v}\left(\mathbb{x},t\right)P\left(\mathbb{x},t\right)\right],$
(2)

with $\mathbb{v}\left(\mathbb{x},t\right)$ being the probability flux given by

$\displaystyle\mathbb{v}\left(\mathbb{x},t\right)=\frac{\mathbb{D}}{k_{B}T}\left[\mathbb{F}\left(\mathbb{x},t\right)-k_{B}T~{}\boldsymbol{\nabla}\log
P\left(\mathbb{x},t\right)\right].$ (3)

This Fokker-Plank equation arises in a broad class of physical systems like
colloidal particles [56], electronic circuits [51] and spin systems [57].
Stochastic thermodynamics dictates that the average total entropy dissipated
when the system runs up to time $t_{f}$ is given by

$\displaystyle
S_{\text{tot}}(t_{f})=k_{B}\int_{0}^{t_{f}}dt~{}\int_{-\infty}^{\infty}d\mathbb{x}~{}P\left(\mathbb{x},t\right)\mathbb{v}\left(\mathbb{x},t\right)\mathbb{D}^{-1}\mathbb{v}\left(\mathbb{x},t\right).$
(4)

Although this formula, in principle, allows one to calculate the entropy
production of any system described by a Fokker-Planck equation, one would need
the force field and probability distribution at each point in space and time.
This is generally not feasible in experimental set-ups.

In this paper, our aim is to develop a technique that gives information about
the total entropy production in terms of experimentally accessible quantities.
Using methods from the optimal transport theory [59, 58, 20], we will derive a
non-trivial lower bound on $S_{\text{tot}}(t_{f})$ completely in terms of the
moments of $\mathbb{x}(t)$. This means, we can infer an estimate for
$S_{\text{tot}}(t_{f})$ simply by measuring the first few moments.

##  3 Bound for one-dimensional systems

Let us first look at the one-dimensional case $(d=1)$ and denote the $n$-th
moment of $x(t)$ by $X_{n}(t)=\langle x(t)^{n}\rangle$. We will first look at
the case where one only has access to the first two moments $X_{1}(t)$ and
$X_{2}(t)$, as the resulting estimate for the entropy production has a rather
simple and elegant expression. Subsequently we will turn to the more general
case.

###  3.1 First two moments

Our central goal is to obtain a lower bound on $S_{\text{tot}}(t_{f})$,
constrained on the known time-dependent moments. This can be done by
minimizing the following action:

$\displaystyle\mathbb{S}\left(x,v,t_{f}\right)$
$\displaystyle=\frac{D}{k_{B}}S_{\text{tot}}(t_{f})+\int_{0}^{t_{f}}dt~{}\left[\mu_{1}(t)X_{1}(t)+\mu_{2}(t)X_{2}(t)\right]+\zeta_{1}(t_{f})X_{1}(t_{f})+\zeta_{2}(t_{f})X_{2}(t_{f}),$
$\displaystyle=\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}dx~{}P(x,t)\left[v^{2}(x,t)+\mu_{1}(t)x+\mu_{2}(t)x^{2}\right]+\int_{-\infty}^{\infty}dx~{}P(x,t_{f})\left[\zeta_{1}(t_{f})x+\zeta_{2}(t_{f})x^{2}\right],$
(5)

with respect to the driving protocol, $F(x,t)$. We have introduced four
Lagrange multipliers, two of which, namely $\zeta_{1}(t_{f})$ and
$\zeta_{2}(t_{f})$ associated with fixing of first two moments at the final
time $t_{f}$ and the other two, $\mu_{1}(t)$ and $\mu_{2}(t)$ at the
intermediate times. In the action, we have included the contribution of
moments at the final time $t_{f}$ separately, as dictated by the Pontryagin’s
Maximum Principle [60]

Minimizing the action with respect to $F(x,t)$ seems very complicated, due to
the non-trivial dependence of $P(x,t)$ on $F(x,t)$ in Eq. (2). To circumvent
this problem, we use methods from the optimal transport theory and introduce a
coordinate transformation $y(x,t)$ as [58]

$\displaystyle\dot{y}(x,t)=v\left(y(x,t),t\right),~{}~{}~{}~{}~{}\text{with
}y(x,0)=x,$ (6)

and overdot indicating the time derivative. Mathematically, this coordinate
transformation is equivalent to changing from Eulerian to Lagrangian
description in fluid mechanics. It has also been used to derive other
thermodynamic bounds such as the thermodynamic speed limit [20]. Using this
equivalence, one can show that [61]

$\displaystyle\int_{-\infty}^{\infty}dx~{}P(x,t)~{}g(x,t)$
$\displaystyle=\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}g(y(x,t),t),$ (7)
$\displaystyle\int_{-\infty}^{\infty}dx~{}P(x,t)~{}v(x,t)~{}g(x,t)$
$\displaystyle=\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}\frac{\partial}{\partial
t}y(x,t)~{}g(y(x,t),t),$ (8)

for any test function $g(x,t)$ and initial probability distribution
$P_{0}(x)=P(x,0)$. Plugging $g(x,t)=\delta\left(x-y(x^{\prime},t)\right)$, we
get

$\displaystyle
P\left(y(x^{\prime},t),t\right)=\frac{P_{0}(x^{\prime})}{|\partial_{x^{\prime}}y(x^{\prime},t)|}.$
(9)

On the other hand, putting $g(x,t)=v(x,t)$ in Eq. (8), we obtain

$\displaystyle
S_{\text{tot}}(t_{f})=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt~{}\int_{-\infty}^{\infty}dx~{}P_{0}\left(x\right)~{}\dot{y}\left(x,t\right)^{2}.$
(10)

Similarly, the moments and their time-derivatives are given by

$\displaystyle X_{n}(t)$
$\displaystyle=\int_{-\infty}^{\infty}dx~{}x^{n}~{}P(x,t)=\int_{-\infty}^{\infty}dx~{}y(x,t)^{n}~{}P_{0}(x),$
(11) $\displaystyle\dot{X}_{n}(t)$
$\displaystyle=n\int_{-\infty}^{\infty}dx~{}x^{n-1}~{}P(x,t)~{}v(x,t)=n\int_{-\infty}^{\infty}dx~{}y(x,t)^{n-1}~{}P(x,t)~{}\dot{y}(x,t).$
(12)

We can also write the force-field $F(x,t)$ in terms of $y(x,t)$ by using Eqs.
(3) and (6). Hence we can calculate various quantities such as probability
distribution, mean total dissipation, and force field utilizing $y(x,t)$. In
fact, it turns out that $y(x,t)$ is uniquely determined for a given force-
field. Using the relations in Eqs. (7) and (8), we can reformulate the
optimisation of the action $\mathbb{S}\left(x,v,t_{f}\right)$ in Eq. (5) with
respect to $F(x,t)$ as an optimisation problem in $y(x,t)$. Meanwhile one can
rewrite $\mathbb{S}\left(x,v,t_{f}\right)$ as

$\displaystyle\mathbb{S}\left(y,\dot{y},t_{f}\right)=$
$\displaystyle\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}dxP_{0}(x)\left[\dot{y}(x,t)^{2}+\mu_{1}(t)y(x,t)+\mu_{2}(t)y(x,t)^{2}\right]$
$\displaystyle+\int_{-\infty}^{\infty}dxP_{0}(x)\left[\zeta_{1}(t_{f})y(x,t_{f})+\zeta_{2}(t_{f})y(x,t_{f})^{2}\right].$
(13)

For a small change $y(x,t)+\delta y(x,t)$, the total change in action is

$\displaystyle\delta\mathbb{S}=$
$\displaystyle\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}\left[-2\ddot{y}(x,t)+\mu_{1}(t)+2\mu_{2}(t)y(x,t)\right]~{}\delta
y(x,t)$
$\displaystyle+\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}\left[2\dot{y}(x,t_{f})+\zeta_{1}(t_{f})+2\zeta_{2}(t_{f})y(x,t_{f})\right]~{}\delta
y(x,t_{f})$
$\displaystyle-2\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}\dot{y}(x,0)~{}\delta
y(x,0).$ (14)

For the optimal protocol, this small change has to vanish. Vanishing of the
first line for arbitrary $\delta y(x,t)$ gives the Euler-Lagrangian equation

$\displaystyle 2\ddot{y}(x,t)=\mu_{1}(t)+2\mu_{2}(t)y(x,t).$ (15)

Let us now analyse the other two terms. Recall that for all possible paths
$y(x,t)$, optimisation is being performed with the fixed initial condition
$y(x,0)=x$ [see Eq. (6)]. Consequently, we have $\delta y(x,0)=0$ which
ensures that the last line in Eq. (14) goes to zero. On the other hand, the
value of $y(x,t_{f})$ at the final time can be different for different paths
which gives $\delta y(x,t_{f})\neq 0$. For the second line to vanish, we must
then have the pre-factor associated with $\delta y(x,t_{f})$ equal to zero.
Thus, we get two boundary conditions in time

$\displaystyle 2\dot{y}(x,t_{f})$
$\displaystyle=-\zeta_{1}(t_{f})-2\zeta_{2}(t_{f})y(x,t_{f}),$ (16)
$\displaystyle y(x,0)$ $\displaystyle=x.$ (17)

Observe that the first equation clearly demonstrates the importance of
incorporating the terms at the final time $t_{f}$ in the action described in
Eq. (5). If these terms are omitted, we would obtain $\dot{y}(x,t_{f})=0$,
which according to Eq. (12) suggests that the time derivative of moments at
the final time is zero. However, this is not true. To get consistent boundary
condition, it is necessary for $\zeta_{1}(t_{f})$ and $\zeta_{2}(t_{f})$ to be
non-zero.

Proceeding to solve Eq. (15), one can verify that

$\displaystyle\dot{y}(x,t)=\lambda_{1}(t)+\lambda_{2}(t)~{}y(x,t),$ (18)

where $\lambda_{1}(t)$ and $\lambda_{2}(t)$ are functions that are related to
the Lagrange multipliers $\mu_{1}(t)$ and $\mu_{2}(t)$. To see this relation,
one can take the time derivative in Eq. (18) and compare it with the Euler-
Lagrange equation (15) to obtain

$\displaystyle\mu_{1}(t)=2\left(\dot{\lambda}_{1}(t)+\lambda_{1}(t)\lambda_{2}(t)\right)~{}~{}\text{and}~{}~{}\mu_{2}(t)=\dot{\lambda}_{2}(t)+\lambda_{2}(t)^{2}.$
(19)

Let us now compute these $\lambda(t)$-functions. For this, we use the solution
of $\dot{y}(x,t)$ in $\dot{X}_{n}(t)$ in Eq. (12) for $n=1$ and $n=2$. This
yields

$\displaystyle\lambda_{1}(t)=\frac{2X_{2}(t)\dot{X}_{1}(t)-X_{1}(t)\dot{X}_{2}(t)}{2(X_{2}(t)-X_{1}^{2}(t))}~{}~{}\text{and}~{}~{}\lambda_{2}(t)=\frac{\dot{X}_{2}(t)-2X_{1}(t)\dot{X}_{1}(t)}{2(X_{2}(t)-X_{1}^{2}(t))}.$
(20)

One can then use Eq. (10) to show that $S_{\text{tot}}(t_{f})$ satisfies the
bound

$\displaystyle
S_{\text{tot}}(t_{f})\geq~{}S_{\text{tot}}^{\text{12}}(t_{f})=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt~{}\left[\frac{\dot{A}_{2}(t)^{2}}{4A_{2}(t)}+\dot{X}_{1}(t)^{2}\right],$
(21)

where $A_{2}(t)=X_{2}(t)-X_{1}(t)^{2}$ is the variance of $x(t)$ and
$S_{\text{tot}}^{\text{12}}(t_{f})$ is the bound on the mean total entropy
dissipated with first and second moments fixed. Eq. (21) is one of the central
results of our paper. It gives a lower bound
$S_{\text{tot}}^{\text{12}}(t_{f})$ on the mean total entropy dissipated
completely in terms of the first moment and variance.

In our analysis, we have focused on minimizing the average total entropy
generated up to time $t_{f}$ with respect to the force field. However, an
alternative approach involves minimizing the entropy production rate at each
intermediate time, as opposed to the total entropy production, but, in A, we
show that this approach generally leads to sub-optimal protocols and therefore
no longer gives a bound on the total entropy production.

###  3.2 Optimal protocol with first two moments

Having obtained a lower bound for the entropy production, it is a natural
question to ask for which protocols this bound is saturated. To answer that
question, we first turn to the flux $v(x,t)$ that gives rise to this optimal
value. For this, we use Eqs. (6) and (18) and obtain

$\displaystyle v(x,t)=v^{\text{12}}(x,t)=\lambda_{1}(t)+\lambda_{2}(t)x,$ (22)

with $\lambda(t)$-functions given in Eq. (20). One can verify that this
equation is always satisfied for Gaussian processes. This is shown in B.

Next, we turn to the probability distribution associated with the optimal
dissipation. For this, one needs the form of $y(x,t)$ which follows from Eq.
(18)

$\displaystyle y(x,t)=\alpha(t)x+\beta(t),$ (23) with
$\displaystyle\alpha(t)=\exp\left(\int_{0}^{t}dt^{\prime}\lambda_{2}(t^{\prime})\right)~{}\text{and
}~{}~{}\beta(t)=\alpha(t)\int_{0}^{t}dt^{\prime}~{}\frac{\lambda_{1}(t^{\prime})}{\alpha(t^{\prime})}.$
(24)

Plugging this in Eq. (9), we obtain the distribution as

$\displaystyle
P(x,t)=P^{\text{12}}(x,t)=\frac{1}{\alpha(t)}P_{0}\left(\frac{x-\beta(t)}{\alpha(t)}\right),$
(25)

which then gives the optimal force-field as

$\displaystyle
F^{\text{12}}\left(x,t\right)=\frac{k_{B}T}{D}\left[\lambda_{1}(t)+\lambda_{2}(t)~{}x+D\frac{\partial}{\partial
x}\ln P_{0}\left(\frac{x-\beta(t)}{\alpha(t)}\right)\right].$ (26)

We see that the optimal protocol comprises solely of a conservative force
field with the associated energy landscape characterized by two components:
the first two terms correspond t a time-dependent harmonic oscillator, whereas
the last term is has the same shape as the equilibrium landscape associated
with the initial state. Furthermore, from Eq. (25), we also observe that this
optimal protocol preserves the shape of the initial probability distribution.
A similar result was also recently observed in the context of the precision-
dissipation trade-off relation for driven stochastic systems [62].

###  3.3 Fixing first $m$ moments

So far, we have derived a bound on $S_{\text{tot}}(t_{f})$ based only on the
knowledge of first two moments. We will now generalise this for arbitrary
number of moments. Consider the general situation where the first $m$ moments
of the variable $x(t)$ are given. Our goal is to optimise the action

$\displaystyle\mathbb{S}\left(x,v,t_{f}\right)$
$\displaystyle=\frac{D}{k_{B}}S_{\text{tot}}(t_{f})+\sum_{i=1}^{m}\int_{0}^{t_{f}}dt~{}\mu_{i}(t)~{}X_{i}(t)+\sum_{i=1}^{m}\zeta_{i}(t_{f})X_{i}(t_{f}),$
(27)
$\displaystyle=\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}dx~{}P(x,t)~{}\left[v^{2}(x,t)+\sum_{i=1}^{m}\mu_{i}(t)~{}x^{i}\right]+\sum_{i=1}^{m}\int_{-\infty}^{\infty}dx~{}P(x,t_{f})~{}\zeta_{i}(t)~{}x^{i},$
(28)

with respect to $v(x,t)$. Once again, $\mu(t)$\- functions are the Lagrange
multipliers for $m$ moments with contribution at the final time $t=t_{f}$
included separately. As done before in Eq. (6), we introduce the new
coordinates, $y(x,t)$, and rewrite the action in terms of this map

$\displaystyle\mathbb{S}\left(y,\dot{y},t_{f}\right)=$
$\displaystyle\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}\left[\dot{y}(x,t)^{2}+\sum_{i=1}^{m}\mu_{i}(t)~{}y(x,t)^{i}\right]+\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}\sum_{i=1}^{m}\zeta_{i}(t_{f})~{}y(x,t_{f})^{i}.$

For a small change $y(x,t)+\delta y(x,t)$, the total change in action is

$\displaystyle\delta\mathbb{S}=$
$\displaystyle\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}\left[-2\ddot{y}(x,t)+\sum_{i=1}^{m}i\mu_{i}(t)~{}y(x,t)^{i-1}\right]~{}\delta
y(x,t)$
$\displaystyle+\int_{-\infty}^{\infty}dx~{}P_{0}(x)\left[2\dot{y}(x,t_{f})+\sum_{i=1}^{m}i~{}\zeta_{i}(t_{f})~{}y(x,t_{f})^{i-1}\right]~{}\delta
y(x,t_{f})$
$\displaystyle~{}~{}~{}~{}~{}~{}-2\int_{-\infty}^{\infty}dx~{}P_{0}(x)\dot{y}(x,0)~{}\delta
y(x,0),$ (29)

which vanishes for the optimal path. Following the same line of reasoning as
Section 3.1, we then obtain the Euler-Lagrange equation

$\displaystyle 2\ddot{y}(x,t)=\sum_{i=1}^{m}i~{}\mu_{i}(t)~{}y(x,t)^{i-1},$
(30)

with two boundary conditions

$\displaystyle 2\dot{y}(x,t_{f})$
$\displaystyle=-\sum_{i=1}^{m}i~{}\zeta_{i}(t_{f})~{}y(x,t_{f})^{i-1},$ (31)
$\displaystyle y(x,0)$ $\displaystyle=x.$ (32)

However, we still need to compute functions $\mu_{i}(t)$ and $\zeta_{i}(t)$
with $i=1,2,..,m$ in order to completely specify the boundary conditions. To
calculate $\zeta$-functions, we use Eq. (12) and obtain

$\displaystyle-\frac{2}{i}\dot{X}_{i}(t_{f})$
$\displaystyle=\zeta_{1}(t_{f})X_{i-1}(t_{f})+2\zeta_{2}(t_{f})X_{i}(t_{f})+...+m\zeta_{m}(t_{f})X_{m-2+i}(t_{f}).$
(33)

On the other hand, for $\mu(t)$-functions, we use the second derivative of
moments and obtain

$\displaystyle\ddot{X}_{i}(t)$
$\displaystyle=\mathcal{B}_{i}(t)+\frac{i}{2}\left[\mu_{1}(t)X_{i-1}(t)+2\mu_{2}(t)X_{i}(t)+...+m\mu_{m}(t)X_{m-2+i}(t)\right],$
(34) $\displaystyle\text{with }\mathcal{B}_{i}(t)$
$\displaystyle=i(i-1)\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}\dot{y}(x,t)^{2}~{}y(x,t)^{i-2}.$
(35)

Both Eqs. (33) and (34) hold for all positive integer values of $i\leq m$.
Solving this gives all $\mu(t)$ and $\zeta(t)$ functions in terms of
$\mathcal{B}_{i}$, $i=1,...,n$ and $X_{i}$, $i=1,...,m+n-2$.

Once Eq. (30) is solved, the bound on average total entropy dissipated can
then be written as

$\displaystyle S_{\text{tot}}(t_{f})\geq
S_{\text{tot}}^{\text{12..m}}(t_{f})=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}\dot{y}(x,t)^{2},$
(36)

Unfortunately, we could not obtain an analytic expression for Eq. (30) due to
the presence of higher-order moments and $\mathcal{B}_{i}$, which are not a
priori known. Instead, we will focus on numerically integrating Eq. (30). One
can use the value of $y(x,t)$ to calculate the higher-order moments and
therefore the $\mu(t)$’s in Eq. (34), which in turn can be used to obtain
$y(x,t+\Delta t)$. Finally, one can make sure that the boundary condition, Eq.
(31), is satisfied through a shooting method.

##  4 Bound in higher dimensions

Having developed a general methodology in one dimension, we now extend these
ideas to higher dimensional systems with state variables
$\mathbb{x}(t)=\\{x_{1}(t),x_{2}(t),....,x_{d}(t)\\}$. Herein again, our aim
is to obtain a bound on $S_{\text{tot}}(t_{f})$ with information about
moments, including mixed moments.

###  4.1 First two moments

As before, let us first introduce the notation

$\displaystyle
X_{1}^{i}(t)=\int_{-\infty}^{\infty}d\mathbb{x}~{}x_{i}~{}P\left(\mathbb{x},t\right),~{}~{}~{}~{}~{}X_{2}^{i,j}(t)=\int_{-\infty}^{\infty}d\mathbb{x}~{}x_{i}x_{j}~{}P\left(\mathbb{x},t\right),$
(37)

for first two (mixed) moments of $x_{i}(t)$ variables. Once again, this is a
constrained optimisation problem and we introduce the following action
functional

$\displaystyle\mathbb{S}\left(\mathbb{x},\mathbb{v},t_{f}\right)$
$\displaystyle=\frac{D}{k_{B}}S_{\text{tot}}(t_{f})+\int_{0}^{t_{f}}dt~{}\left[\sum_{i=1}^{d}\mu_{i}(t)X_{1}^{i}(t)+\sum_{i,j=1}^{d}\lambda_{ij}(t)~{}X_{2}^{i,j}(t)\right]~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{i=1}^{d}\alpha_{i}(t_{f})X_{1}^{i}(t_{f})+\sum_{i,j=1}^{d}\beta_{ij}(t_{f})~{}X_{2}^{i,j}(t_{f}),$
(38)
$\displaystyle=\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}d\mathbb{x}~{}P\left(\mathbb{x},t\right)\left[\mathbb{v}\left(\mathbb{x},t\right)^{2}+\sum_{i=1}^{d}\mu_{i}(t)~{}x_{i}+\sum_{i,j=1}^{d}\lambda_{ij}(t)~{}x_{i}x_{j}\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\int_{-\infty}^{\infty}d\mathbb{x}~{}P\left(\mathbb{x},t_{f}\right)\left[\sum_{i=1}^{d}\alpha_{i}(t_{f})~{}x_{i}+\sum_{i,j=1}^{d}\beta_{ij}(t_{f})~{}x_{i}x_{j}\right],$
(39)

where functions $\mu_{i}(t)$, $\lambda_{ij}(t)$, $\alpha_{i}(t)$ and
$\beta_{ij}(t)$ all stand for Lagrange multipliers associated with given
constraints. Further carrying out this optimisation is challenging since
distribution $P\left(\mathbb{x},t\right)$ depends non-trivially on the force-
field $\mathbb{F}\left(\mathbb{x},t\right)$. We, therefore, introduce a
$d$-dimensional coordinates

$\displaystyle\dot{\mathbb{y}}(\mathbb{x},t)=\mathbb{v}\left(\mathbb{y}(\mathbb{x},t),t\right),~{}~{}~{}~{}~{}\text{with
}\mathbb{y}(\mathbb{x},0)=\mathbb{x},$ (40)

and rewrite $\mathbb{S}\left(\mathbb{x},\mathbb{v},t_{f}\right)$ completely in
terms of $\mathbb{y}(\mathbb{x},t)$. We then follow same steps as done in
Section 3.1 but in higher dimensional setting. Since this is a straightforward
generalisation, we relegate this calculation to C and present only the final
results here. For the optimal protocol, we find

$\displaystyle\dot{y}_{i}\left(\mathbb{x},t\right)=\dot{X}_{1}^{i}(t)+\sum_{j=1}^{d}\frac{\gamma_{ij}(t)}{2}~{}\left[y_{j}\left(\mathbb{x},t\right)-X_{1}^{j}(t)\right],$
(41)

with functions $\gamma_{ij}(t)$ related to the Lagrange multipliers
$\mu_{i}(t)$, $\lambda_{ij}(t)$ via Eq. (128). Furthermore, in C, we show that
$\gamma_{ij}(t)$ functions are the solutions of a set of linear equations:

$\displaystyle-2\dot{\mathbb{A}}(t)=\boldsymbol{\gamma}(t)\mathbb{A}(t)+\mathbb{A}(t)\boldsymbol{\gamma}(t).$
(42)

where $\boldsymbol{\gamma}(t)$ and $\mathbb{A}(t)$ are two matrices with
elements $\gamma_{ij}(t)$ and
$A_{ij}(t)=X_{2}^{i,j}(t)-X_{1}^{i}(t)X_{1}^{j}(t)$ respectively. Notice that
the elements of $\mathbb{A}(t)$ comprise only of variance and covariance which
are known. Therefore, the matrix $\mathbb{A}(t)$ is fully known at all times.

To solve Eq. (42), we note that $\mathbb{A}(t)$ is a symmetric matrix and
thus, there exists an orthogonal matrix $\mathbb{O}(t)$ such that

$\displaystyle\mathbb{A}(t)=\mathbb{O}(t)~{}\boldsymbol{\Lambda}(t)~{}\mathbb{O}^{T}(t),$
(43)

where $\boldsymbol{\Lambda}(t)$ is a diagonal matrix whose elements correspond
to the positive eigenvalues of $\mathbb{A}(t)$. $\mathbb{O}(t)$ can be
constructed using eigenvectors of $\mathbb{A}(t)$ and it satisfies the
orthogonality condition $\mathbb{O}^{-1}(t)=\mathbb{O}^{T}(t)$. Combining the
transformation Eq. (43) with Eq. (42), we obtain

$\displaystyle{\mathcal{G}}_{ij}(t)$
$\displaystyle=-\frac{2~{}\left[\mathbb{O}^{T}(t)~{}\dot{\mathbb{A}}~{}\mathbb{O}(t)\right]_{ij}}{{\Lambda}_{ii}(t)+{\Lambda}_{jj}(t)},~{}~{}~{}~{}\text{with
}\boldsymbol{\mathcal{G}}(t)=\mathbb{O}^{T}(t)~{}\boldsymbol{\gamma}(t)~{}\mathbb{O}(t).$
(44)

This relation enables us to write $\dot{y}_{i}\left(\mathbb{x},t\right)$ in
Eq. (41) as

$\displaystyle\dot{y}_{i}\left(\mathbb{x},t\right)=\dot{X}_{1}^{i}(t)+\sum_{j=1}^{d}\frac{\left[\mathbb{O}(t)~{}\boldsymbol{\mathcal{G}}(t)~{}\mathbb{O}^{T}(t)\right]_{ij}}{2}~{}\left[y_{j}\left(\mathbb{x},t\right)-X_{1}^{j}(t)\right].$
(45)

Finally, plugging this in the formula

$\displaystyle S_{\text{tot}}(t_{f})$
$\displaystyle=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt~{}\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}(\mathbb{x})~{}\dot{\mathbb{y}}\left(\mathbb{x},t\right)^{2},$
(46)

we obtain the bound in general $d$ dimensions as

$\displaystyle S_{\text{tot}}^{12}(t_{f})$
$\displaystyle=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt\left[\dot{\mathbb{X}}_{1}(t)^{2}+\frac{1}{4}\text{Trace}\Big{\\{}\mathbb{O}(t)~{}\boldsymbol{\mathcal{G}}^{2}(t)~{}\mathbb{O}^{T}(t)~{}\mathbb{A}(t)\Big{\\}}\right],$
(47)
$\displaystyle=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt\left[\dot{\mathbb{X}}_{1}(t)^{2}+\frac{1}{4}\text{Trace}\Big{\\{}\boldsymbol{\mathcal{G}}^{2}(t)~{}\boldsymbol{\Lambda}(t)~{}\Big{\\}}\right],$
(48)

where the elements of $\boldsymbol{\mathcal{G}}(t)$ matrix are fully given in
terms of $\mathbb{A}(t)$ in Eq. (44) which in turn depends on the variance and
covariance of $\mathbb{x}(t)$. Eq. (48) represents a general bound on the mean
entropy production in higher dimensions. For the two-dimensional case, we
obtain an explicit expression (see C.1)

$\displaystyle S_{\text{tot}}^{12}(t_{f})$
$\displaystyle=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt\left[\dot{X}_{1}^{1}(t)^{2}+\frac{\dot{A}_{11}(t)^{2}}{4A_{11}(t)}+\dot{X}_{1}^{2}(t)^{2}+\frac{\dot{A}_{22}(t)^{2}}{4A_{22}(t)}\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+\frac{\left\\{A_{12}(t)\frac{d}{dt}\left(A_{11}(t)A_{22}(t)\right)-2A_{11}(t)A_{22}(t)\dot{A}_{12}(t)\right\\}^{2}}{4A_{11}(t)A_{22}(t)\left\\{A_{11}(t)+A_{22}(t)\right\\}\left\\{A_{11}(t)A_{22}(t)-A_{12}(t)^{2}\right\\}}\right].$
(49)

The bound is equal to the sum of two one dimensional bounds in Eq. (21)
corresponding to two coordinates $x_{1}(t)$ and $x_{2}(t)$ plus a part that
arises due to the cross-correlation between them. In absence of this cross-
correlation $\left(A_{12}(t)=0\right)$, Eq. (49) correctly reduces to the sum
of two separate one dimensional bounds.

###  4.2 Optimal protocols with first two moments

Let us now look at the flux and the distribution associated with this optimal
dissipation. For this, we first rewrite Eq. (41) in vectorial form as

$\displaystyle\dot{\mathbb{y}}\left(\mathbb{x},t\right)=\dot{\mathbb{X}}_{1}(t)+\frac{\boldsymbol{\gamma}(t)}{2}\left[\mathbb{y}\left(\mathbb{x},t\right)-\mathbb{X}_{1}(t)\right],$
(50)

and solve it to obtain the solution $\mathbb{y}\left(\mathbb{x},t\right)$ as

$\displaystyle\mathbb{y}\left(\mathbb{x},t\right)=\mathbb{X}_{1}(t)+\boldsymbol{\Omega}(t)\left[\mathbb{x}-\mathbb{X}_{1}(t)\right].$
(51)

Here $\boldsymbol{\Omega}(t)$ stands for the time-ordered exponential and is
given by

$\displaystyle\boldsymbol{\Omega}(t)$
$\displaystyle=\mathcal{T}\Big{[}\exp\left(\int_{0}^{t}dt^{\prime}~{}\frac{\boldsymbol{\gamma}(t^{\prime})}{2}\right)\Big{]},$
$\displaystyle=\sum_{n=0}^{\infty}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\int_{0}^{t_{2}}dt_{3}...\int_{0}^{t_{n-1}}dt_{n}~{}\frac{\boldsymbol{\gamma}(t_{1})}{2}\frac{\boldsymbol{\gamma}(t_{2})}{2}....\frac{\boldsymbol{\gamma}(t_{n})}{2},$
(52)

with $\boldsymbol{\gamma}(t)$ defined in Eq. (44) and $\mathcal{T}$
corresponding to the time-ordering operator. Finally, substituting the
solution $\mathbb{y}\left(\mathbb{x},t\right)$ from Eq. (51) in Eqs. (40) and
(120), we obtain the optimal flux and optimal distribution as

$\displaystyle\mathbb{v}^{12}\left(\mathbb{x},t\right)$
$\displaystyle=\dot{\mathbb{X}}_{1}(t)+\frac{\boldsymbol{\gamma}(t)}{2}\left[\mathbb{x}-\mathbb{X}_{1}(t)\right],$
(53) $\displaystyle P^{12}\left(\mathbb{x},t\right)$
$\displaystyle=\frac{P_{0}\Big{(}\mathbb{X}_{1}(t)+\boldsymbol{\Omega}^{-1}(t)\left[\mathbb{x}-\mathbb{X}_{1}(t)\right]\Big{)}}{|\text{det}\left(\boldsymbol{\Omega}(t)\right)|}.$
(54)

Now the optimal force-field follows straightforwardly from Eq. (3) as

$\displaystyle\mathbb{F}\left(\mathbb{x},t\right)=\frac{k_{B}T}{D}\left[\dot{\mathbb{X}}_{1}(t)+\frac{\boldsymbol{\gamma}(t)}{2}\Bigg{(}\mathbb{x}-\mathbb{X}_{1}(t)\Bigg{)}+D~{}\boldsymbol{\nabla}\ln
P_{0}\Big{(}\mathbb{X}_{1}(t)+\boldsymbol{\Omega}^{-1}(t)\left[\mathbb{x}-\mathbb{X}_{1}(t)\right]\Big{)}\right].$
(55)

To sum up, we have derived a lower bound to the total dissipation
$S_{\text{tot}}(t_{f})$ in general $d$ dimensions just from the information
about first two moments. Together with this, we also calculated the force-
field that gives rise to this optimal value. In what follows, we consider the
most general case where we know the first-$m$ (mixed) moments of the state
variable $\mathbb{x}(t)=\\{x_{1}(t),x_{2}(t),...,x_{d}(t)\\}$.

###  4.3 Fixing first $m$ moments in general $d$ dimensions

Let us define a general mixed moment as

$\displaystyle
X_{m}^{i_{1},i_{2},..,i_{m}}(t)=\int_{-\infty}^{\infty}d\mathbb{x}~{}P\left(\mathbb{x},t\right)~{}x_{i_{1}}x_{i_{2}}....x_{i_{3}},$
(56)

where $i_{j}\in\mathbb{Z}^{+}$ and $i_{1}\geq i_{2}\geq i_{3}\geq..\geq i_{m}$
for all $1\leq i_{l}\leq d$. For $m=1$ and $m=2$, Eq. (56) reduces to first
two moments in Eq. (37). Here, we are interested in the general case for which
we consider the following action

$\displaystyle\mathbb{S}\left(\mathbb{x},\mathbb{v},t_{f}\right)$
$\displaystyle=\frac{D}{k_{B}}S_{\text{tot}}(t_{f})+\int_{0}^{t_{f}}dt\sum_{l=1}^{m}\sum_{~{}\\{i_{j}\\}=1}^{d}~{}\mu_{i_{1},i_{2},...,i_{l}}(t)~{}X_{l}^{i_{1},i_{2},..,i_{l}}(t)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{l=1}^{m}\sum_{~{}\\{i_{j}\\}=1}^{d}~{}\zeta_{i_{1},i_{2},...,i_{l}}(t_{f})~{}X_{l}^{i_{1},i_{2},..,i_{l}}(t_{f}).$
(57)

In conjunction to the previous sections, we again recast this action in terms
of $\mathbb{y}\left(\mathbb{x},t\right)$ defined in Eq. (40) as

$\displaystyle\mathbb{S}\left(\mathbb{y},\dot{\mathbb{y}},t_{f}\right)$
$\displaystyle=\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)\left[\dot{\mathbb{y}}\left(\mathbb{x},t\right)^{2}+\sum_{l=1}^{m}\sum_{~{}\\{i_{j}\\}=1}^{d}~{}\mu_{i_{1},i_{2},...,i_{l}}(t)\left\\{\prod_{q=1}^{l}y_{i_{q}}\left(\mathbb{x},t\right)\right\\}\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}+\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)\sum_{l=1}^{m}\sum_{~{}\\{i_{j}\\}=1}^{d}~{}\zeta_{i_{1},i_{2},...,i_{l}}(t_{f})\left\\{\prod_{q=1}^{l}y_{i_{q}}\left(\mathbb{x},t_{f}\right)\right\\}.$
(58)

Optimising this action in the same way as before, we obtain the Euler-Lagrange
equation

$\displaystyle 2\ddot{y}_{i}\left(\mathbb{x},t\right)$
$\displaystyle=\sum_{l=1}^{m}\sum_{~{}\\{i_{j}\\}=1}^{d}~{}\mu_{i,i_{2},...,i_{l}}(t)~{}y_{i_{2}}\left(\mathbb{x},t\right)y_{i_{3}}\left(\mathbb{x},t\right)...y_{i_{l}}\left(\mathbb{x},t\right),$
(59)

along with the boundary conditions

$\displaystyle 2\dot{y}_{i}\left(\mathbb{x},t_{f}\right)$
$\displaystyle=-\sum_{l=1}^{m}\sum_{~{}\\{i_{j}\\}=1}^{d}~{}\zeta_{i,i_{2},...,i_{l}}(t_{f})~{}y_{i_{2}}\left(\mathbb{x},t_{f}\right)y_{i_{3}}\left(\mathbb{x},t_{f}\right)...y_{i_{l}}\left(\mathbb{x},t_{f}\right),$
(60) $\displaystyle y_{i}\left(\mathbb{x},0\right)$ $\displaystyle=x_{i}.$
(61)

The Lagrange multipliers $\mu_{i_{1},i_{2},...,i_{l}}(t)$ and
$\zeta_{i_{1},i_{2},...,i_{l}}(t)$ can be computed using time-derivatives of
the (mixed) moments $X_{l}^{i_{1},i_{2},..,i_{l}}(t)$ as

$\displaystyle\dot{X}_{r}^{i_{1},i_{2},..,i_{r}}(t)$
$\displaystyle=-\frac{1}{2}\sum_{k=1}^{r}\sum_{l=1}^{m}\sum_{~{}\\{i_{j}^{\prime}\\}=1}^{d}~{}\zeta_{i_{k},i_{2}^{\prime},i_{3}^{\prime},...,i_{l}^{\prime}}(t)~{}X_{r+l-2}^{i_{1},i_{2},...,i_{k-1},i_{k+1}.,i_{r},i_{2}^{\prime}.,i_{3}^{\prime},...,i_{l}^{\prime}}(t)$
(62) $\displaystyle\ddot{X}_{r}^{i_{1},i_{2},..,i_{r}}(t)$
$\displaystyle=\frac{1}{2}\sum_{k=1}^{r}\sum_{l=1}^{m}\sum_{~{}\\{i_{j}^{\prime}\\}=1}^{d}~{}\mu_{i_{k},i_{2}^{\prime},i_{3}^{\prime},...,i_{l}^{\prime}}(t)~{}X_{r+l-2}^{i_{1},i_{2},...,i_{k-1},i_{k+1}.,i_{r},i_{2}^{\prime}.,i_{3}^{\prime},...,i_{l}^{\prime}}(t)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\mathcal{B}_{i_{1},i_{2},..,i_{r}}(t),$
(63) $\displaystyle\text{with }\mathcal{B}_{i_{1},i_{2},..,i_{r}}(t)$
$\displaystyle=\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)\sum_{k,l=1\atop
k\neq
l}^{r}~{}\frac{\dot{y}_{i_{l}}\left(\mathbb{x},t\right)~{}\dot{y}_{i_{k}}\left(\mathbb{x},t\right)}{y_{i_{l}}\left(\mathbb{x},t\right)~{}y_{i_{k}}\left(\mathbb{x},t\right)}\prod_{j=1}^{r}y_{i_{j}}\left(\mathbb{x},t\right).$
(64)

Once again, for $m>2$, we find dependence of the right hand side of Eqs. (62)
and (63) on (mixed) moments with order higher than $m$. As discussed before in
Section 3.3, one can tackle this by utilizing the solution
$y_{i}\left(\mathbb{x},t\right)$ recursively to obtain these higher moments
[as discussed after Eq. (35)].

Solving Eq. (59) numerically yields the solution
$\mathbb{y}_{*}\left(\mathbb{x},t\right)$ using which in Eq. (46) we obtain a
lower bound to $S_{\text{tot}}(t_{f})$. The associated optimal flux and
optimal distribution are calculated using Eqs. (40) and (120) respectively.

##  5 Theoretical examples

We now test our framework on two analytically solvable toy models. The first
one is consists of a one dimensional free diffusion model with initial
position drawn from the distribution $P_{0}(x_{0})\sim\exp(-x_{0}^{4})$. This
example will demonstrate how our method becomes more accurate when we
incorporate knowledge of higher moments. The second example involves two-
dimensional diffusion in a moving trap, serving as an illustration of our
method’s applicability in higher dimensions.

Figure 1: (a) Estimation of the total entropy production
$S_{\text{tot}}(t_{f})$ for model (65) by using the first two moments (in
green) and first four moments (in blue). The red curve is the exact result
whose expression is given in Eq. (67). (b) Probability distribution of the
total entropy dissipated $s_{\text{tot}}(t_{f})$ for the optimal protocol
obtained by fixing the first two moments (in green) and first four moments (in
blue) till time $t_{f}=1$. Here also, the red symbols correspond to the
simulation of Eq. (65). For both panels, we have used $D=1\text{ and }k_{B}=1$

###  5.1 Free diffusion with $P_{0}(x_{0})=\exp(-x_{0}^{4})/2\Gamma(5/4)$

We first consider a freely diffusing particle in one dimension

$\displaystyle\frac{dx}{dt}=\sqrt{2D}~{}\eta(t),~{}~{}~{}\text{with
}P_{0}(x_{0})=\exp(-x_{0}^{4})/2\Gamma(5/4).$ (65)

This example serves two purposes: First, we can obtain the exact expression of
$S_{\text{tot}}(t_{f})$ which enables us to rigorously compare the bounds
derived above with their exact counterpart. Second, we discuss the intricacies
that arise when we fix more than first two moments. To begin with, the
position distribution $P(x,t_{f})$ is given by

$\displaystyle P(x,t_{f})=\frac{1}{\sqrt{4\pi
Dt_{f}}}~{}\int_{-\infty}^{\infty}dx_{0}~{}P_{0}(x_{0})~{}\exp\left(-\frac{(x-x_{0})^{2}}{4Dt_{f}}\right).$
(66)

Since no external force acts on the system, the total entropy produced will be
equal to the total entropy change of the system. Therefore, the mean
$S_{\text{tot}}(t_{f})$ reads

$\displaystyle
S_{\text{tot}}(t_{f})=-k_{B}\int_{-\infty}^{\infty}dx~{}P(x,t_{f})~{}\log
P(x,t_{f})+k_{B}\int_{-\infty}^{\infty}dx_{0}~{}P_{0}(x_{0})~{}\log
P_{0}(x_{0}).$ (67)

We emphasize that Eqs. (66) and (67) are exact results and do not involve any
approximation.

Let us now analyse how our bound compares with this exact value. The first
four moments of the position are given by

$\displaystyle X_{1}(t)=X_{3}(t)=0,$ (68) $\displaystyle
X_{2}(t)=2Dt+\frac{\Gamma(3/4)}{\Gamma(1/4)},$ (69) $\displaystyle
X_{4}(t)=\frac{1}{4}+12~{}\frac{\Gamma(3/4)}{\Gamma(1/4)}Dt+12D^{2}t^{2}.$
(70)

Combining this with Eq. (36), shows that the dissipation bound with first two
moments reads

$\displaystyle
S_{\text{tot}}^{\text{12}}(t_{f})=\int_{0}^{t_{f}}dt\frac{k_{B}D}{2Dt+\frac{\Gamma(3/4)}{\Gamma(1/4)}}=\frac{k_{B}}{2}\log\left(\frac{2Dt_{f}+\frac{\Gamma(3/4)}{\Gamma(1/4)}}{\frac{\Gamma(3/4)}{\Gamma(1/4)}}\right),$
(71)

with associated optimal protocols and time-dependent distribution given by

$\displaystyle P^{12}(x,t)$
$\displaystyle=\frac{1}{2\Gamma\left(5/4\right)~{}\alpha(t)}~{}\exp\left[-x^{4}/\alpha(t)^{4}\right],~{}~{}~{}\text{with
}\alpha(t)=\sqrt{\frac{\frac{\Gamma\left(3/4\right)}{\Gamma\left(1/4\right)}}{2Dt+\frac{\Gamma\left(3/4\right)}{\Gamma\left(1/4\right)}}}$
(72) $\displaystyle v^{12}(x,t)$
$\displaystyle=\frac{Dx}{2Dt+\frac{\Gamma\left(3/4\right)}{\Gamma\left(1/4\right)}}.$
(73) $\displaystyle F^{12}(x,t)$
$\displaystyle=k_{B}T\left[\frac{x}{2Dt+\frac{\Gamma\left(3/4\right)}{\Gamma\left(1/4\right)}}-\frac{4x^{3}}{\alpha(t)^{4}}\right],$
(74)

As illustrated in Figure 1(a), the bound $S_{\text{tot}}^{\text{12}}(t_{f})$,
although not exact, is still remarkably close to the exact value in Eq. (67).
We have also plotted the full probability distribution of total dissipation in
Figure 1(b) using the optimal force-field in Eq. (74). Here, we see deviation
from the exact result indicating that fixing first two moments only does not
give a very good approximation of dissipation beyond its average value. To
improve this, we instead fix the first four moments of position and carry out
this analysis.

One can wonder whether it is possible to get an even better estimate by taking
into account the first four moments. To check this, we will now proceed by
solving the Euler-Lagrange equation (30) numerically for $m=4$. Since, this is
a second-order differential equation, we need two boundary conditions in time.
One of them is Eq. (32) which gives $y(x,0)=x$. However, the other one in Eq.
(31) gives a condition at final time $t_{f}$ which is difficult to implement
numerically. For this, we use the following trick. We expand the initial
$\dot{y}(x,0)$ as

$\displaystyle\dot{y}(x,0)=\theta_{0}+\theta_{1}x+....+\theta_{s}x^{s},$ (75)

where we typically choose $s=9$. The coefficients $\theta_{i}$ are free to
vary but must give the correct time derivative of first four moments in Eqs.
(68)-(70). We then assign a cost function as

$\displaystyle\mathcal{C}(t_{f})=\int_{-\infty}^{\infty}dx~{}P_{0}(x)~{}\left[2\dot{y}(x,t_{f})+\sum_{i=1}^{4}i~{}\zeta_{i}(t_{f})~{}y(x,t_{f})^{i-1}\right]^{2}.$
(76)

Observe that the cost function vanishes if the boundary condition (31) is
satisfied. Starting from $y(x,0)=x$ and $\dot{y}(x,0)$ in Eq. (75), we solve
the Euler-Lagrange equation (30) numerically and measure the cost function at
the final time. We then perform a multi-dimensional gradient descent in
$\theta$-parameters to minimise $\mathcal{C}(t_{f})$. Eventually, we obtain
the numerical form of the optimal map $y_{*}(x,t)$.

Using this map, we then obtain a refined bound on the dissipation,
$S_{\text{tot}}^{1234}(t_{f})$. This is illustrated in Figure 1(a). Clearly,
this bound is better than $S_{\text{tot}}^{12}(t_{f})$ and in fact, matches
very well with the exact one. Furthermore, we can also obtain the associated
optimal protocols by numerically inverting Eqs. (6) and (9). With these
protocols, we then obtain the distribution of total entropy dissipated which
is shown in Figure 1(b). Compared to the previous case, we again find that the
estimated distribution is now closer to the exact one.

###  5.2 Two dimensional Brownian motion

The previous example illustrated our method in a one-dimensional setting. In
this section, we look at a two-dimensional system with position
$\mathbb{x}(t)=\left(x_{1}(t),x_{2}(t)\right)$. The system undergoes motion in
presence of a moving two-dimensional harmonic oscillator
$U(x_{1},x_{2},t)=K\left(x_{1}-vt\right)\left(x_{2}-vt\right)$ and its
dynamics is governed by

$\displaystyle\frac{dx_{1}}{dt}=-K(x_{2}-vt)+\sqrt{2D}~{}\eta_{1}(t),$ (77)
$\displaystyle\frac{dx_{2}}{dt}=-K(x_{1}-vt)+\sqrt{2D}~{}\eta_{2}(t).$ (78)

We also assume that the initial position is drawn from Gaussian distribution
as

$\displaystyle P_{0}(x_{1}(0),x_{2}(0))=\frac{1}{2\pi
b_{0}^{2}}~{}\exp\left(-\frac{x_{1}(0)^{2}+x_{2}(0)^{2}}{2b_{0}^{2}}\right).$
(79)

For this example, our aim is to compare the bound derived in Eq. (49) with the
exact expression. Using the Langevin equation, we obtain the first two
cumulants and covariance to be

$\displaystyle
X_{1}^{1}(t)=X_{1}^{2}(t)=X_{1}(t)=\frac{v}{K}\left(Kt-1+e^{-Kt}\right),$ (80)
$\displaystyle A_{11}(t)=A_{22}(t)=b_{0}^{2}\cosh(2Kt)+\frac{D}{K}\sinh(2Kt),$
(81) $\displaystyle
A_{12}(t)=A_{21}(t)=-b_{0}^{2}\sinh(2Kt)+\frac{D}{K}\left[1-\cosh(2Kt)\right].$
(82)

Plugging them in Eq. (49) gives

$\displaystyle S^{12}_{\text{tot}}\left(t_{f}\right)$
$\displaystyle=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt\left[2\dot{X}_{1}(t)^{2}+\frac{\Big{(}\dot{A}_{11}(t)+\dot{A}_{12}(t)\Big{)}^{2}}{4\Big{(}A_{11}(t)+A_{12}(t)\Big{)}}+\frac{\Big{(}\dot{A}_{11}(t)-\dot{A}_{12}(t)\Big{)}^{2}}{4\Big{(}A_{11}(t)-A_{12}(t)\Big{)}}\right],$
(83)
$\displaystyle=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt\left[2v^{2}(1-e^{-Kt})^{2}+\frac{\left(Kb_{0}^{2}-D\right)^{2}e^{-4Kt}}{b_{0}^{2}e^{-2Kt}+\frac{D}{K}\left(1-e^{-2Kt}\right)}+\frac{\left(Kb_{0}^{2}+D\right)^{2}e^{4Kt}}{b_{0}^{2}e^{2Kt}+\frac{D}{K}\left(e^{2Kt}-1\right)}\right].$
(84)

In fact, one can carry out exact calculations for this model. To see this, let
us introduce the change of variable

$\displaystyle
R_{1}=\frac{x_{1}+x_{2}}{\sqrt{2}},~{}~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}~{}R_{2}=\frac{x_{1}-x_{2}}{\sqrt{2}}$
(85)

and rewrite the Langevin equations as

$\displaystyle\frac{dR_{1}}{dt}=-K\left(R_{1}-\sqrt{2}vt\right)+\sqrt{2D}~{}\eta_{R_{1}}(t),$
(86) $\displaystyle\frac{dR_{2}}{dt}=KR_{2}+\sqrt{2D}~{}\eta_{R_{2}}(t),$ (87)

where $\eta_{R_{1}}(t)$ and $\eta_{R_{2}}(t)$ are two independent Gaussian
white noises with zero mean and delta correlation. Clearly these two equations
are independent of each other which means that we can treat $R_{1}(t)$ and
$R_{2}(t)$ as two independent one dimensional Brownian motions. Also, observe
that both of them are Gaussian processes. Following B, we can now express
average dissipation completely in terms of first two cumulants. The first two
cumulants for our model are

$\displaystyle\langle
R_{1}(t)\rangle=\frac{\sqrt{2}v}{K}\left(Kt-1+e^{-Kt}\right),~{}~{}~{}~{}~{}\langle
R_{2}(t)\rangle=0,$ (88) $\displaystyle A_{R_{1}}(t)=\langle
R_{1}^{2}(t)\rangle-\langle
R_{1}(t)\rangle^{2}=b_{0}^{2}~{}e^{-2Kt}+\frac{D}{K}\left(1-e^{-2Kt}\right),$
(89) $\displaystyle A_{R_{2}}(t)=\langle R_{2}^{2}(t)\rangle-\langle
R_{2}(t)\rangle^{2}=b_{0}^{2}~{}e^{2Kt}+\frac{D}{K}\left(e^{2Kt}-1\right).$
(90)

Finally inserting them in Eq. (21), we obtain

$\displaystyle S_{\text{tot}}\left(t_{f}\right)$
$\displaystyle=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt\left[\left(\frac{d\langle
R_{1}(t)\rangle}{dt}\right)^{2}+\frac{\dot{A}_{R_{1}}(t)^{2}}{4A_{R_{1}}(t)}+\frac{\dot{A}_{R_{2}}(t)^{2}}{4A_{R_{2}}(t)}\right],$
$\displaystyle=\frac{k_{B}}{D}\int_{0}^{t_{f}}\left[2v^{2}(1-e^{-Kt})^{2}+\frac{\left(Kb_{0}^{2}-D\right)^{2}e^{-4Kt}}{b_{0}^{2}e^{-2Kt}+\frac{D}{K}\left(1-e^{-2Kt}\right)}+\frac{\left(Kb_{0}^{2}+D\right)^{2}e^{4Kt}}{b_{0}^{2}e^{2Kt}+\frac{D}{K}\left(e^{2Kt}-1\right)}\right],$
(91)

which matches with the estimated entropy production in Eq. (84) exactly. In
other words, our method becomes exact in this case.

##  6 Bit erasure

Figure 2: Left panel shows the control functions, $g(t)$, $f(t)$ and $m(t)$,
that are modulated experimentally to carry out the erasure operation. The
trajectories thus obtained are shown in the right panel. The grey trajectory
starts from the left well initially and ends up in the right well at the end
of the cycle. On the other hand, the red trajectory starts from the right well
and also ends up in the right well.

Motivated by the very good agreement between the bound and the actual entropy
production in simple theoretical models, we will now turn to experimental data
of a more complicated system, namely bit erasure [63]. The state of a bit can
often be described by a one-dimensional Fokker-Planck equation, with a double-
well potential energy landscape. This landscape can for example be produced by
an optical tweezer [64], a feed-back trap [50], or a magnetic field [65, 66].
Initially, the colloidal particle can either be in the left well (state $0$)
or in the right well (state $1$) of the double-well potential. By modulating
this potential, one can assure that the particle always ends up the right well
(state $1$). During this process, the expected amount of work put into the
system is always bounded by the Landauer’s limit, $\langle W\rangle\geq
k_{B}T\ln 2$ [50]. Over the last few years, Landauer’s principle has also been
tested in a more complex situation where symmetry between two states is broken
[67, 68]. To add an extra layer of complexity, we consider here the scenario
where this symmetry is broken by deploying an asymmetric double-well potential
and we will use the data from [67] to test our bound. The system satisfies the
Fokker-Planck equation

$\displaystyle\frac{\partial P(x,t)}{\partial
t}=\frac{D}{k_{B}T}\frac{\partial}{\partial x}\left[P(x,t)~{}\frac{\partial
U(x,t)}{\partial x}+k_{B}T~{}\frac{\partial P(x,t)}{\partial x}\right],$ (92)

where $U(x,t)$ is the time-dependent asymmetric double-well potential given by

$\displaystyle U(x,t)$
$\displaystyle=4E_{b}\left[-\frac{g(t)}{2}\tilde{x}^{2}+\frac{\tilde{x}^{4}}{4}-Af(t)~{}\tilde{x}\right],$
(93) $\displaystyle\text{with }\tilde{x}=$
$\displaystyle\begin{cases}\frac{x}{x_{m}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if
}x<0,\\\
\frac{x}{x_{m}}\left[1+m(t)(\eta-1)\right]^{-1}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if
}x\geq 0.\end{cases}$ (94)

Notice that initially, the potential $U(x,t)$ has its minimum located at
$-x_{m}$ and $+\eta x_{m}$, where $\eta\geq 1$ is the asymmetry factor. For
$\eta=1$, the potential is completely symmetric initially. On the other hand,
for $\eta\neq 0$, the two wells are asymmetric and the system is out of
equilibrium initially. Furthermore, the functions $m(t)$, $g(t)$ and $f(t)$
represent the experimental protocols to modulate $U(x,t)$ during erasure
operation. They are defined in Figure 2 and in Table 1. The experimental
protocol consists of symmetrizing $U(x,t)$ by changing the function $m(t)$,
while the other two functions $g(t)$ and $f(t)$ are kept fixed. Then, $g(t)$
and $f(t)$ are changed during which the barrier is lowered, the potential is
tilted and again brought back to its symmetric form at time $t=3t_{f}/4$. By
this time, the particle is always at the right well. Finally, the right well
is expanded to its original size with minimum located at $+\eta x_{m}$ at the
end of the cycle time $t_{f}$. Experimental trajectories thus obtained are
shown in Figure 2 (right panel).

Figure 3: First two moments associated with the erasure operation with the
cycle time $t_{f}/t_{m}=14.89$. The measured moments are fitted by the
Chebyshev polynomials as indicated in Eqs. (95) and (96). Figure 4: Estimated
mean total entropy produced $S_{\text{tot}}$ as a function of the rescaled
cycle time $t_{f}/t_{m}$. We have also performed a comparison with the exact
values of $S_{\text{tot}}$ obtained by simulating the Langevin equation with
potential $U(x,t)$ given in Eq. (93). The error bars are calculated through
least-squares fitting to moments $X_{1}(t)$ and $X_{2}(t)$.

Here, we are interested in estimating the mean entropy produced during the
erasure operation by using the bound in Eq. (21). For this, we take $20-50$
trajectories for different cycle times $t_{f}$ and measure the first two
moments. The precise values of cycle time and number of trajectories
considered are given in Table 2. We fit the measured moments using the basis
of Chebyshev polynomials as

$\displaystyle X_{1}(t)=\langle
x(t)\rangle=\chi_{0}T_{0}(t)+\chi_{1}T_{1}(t)+....+\chi_{s}T_{s}(t),$ (95)
$\displaystyle X_{2}(t)=\langle
x^{2}(t)\rangle=\psi_{0}T_{0}(t)+\psi_{1}T_{1}(t)+....+\psi_{s}T_{s}(t),$ (96)

where we typically choose $s=20$. For our analysis, we have rescaled the
position by $x\to x/x_{m}$ and $t\to t/t_{m}$ with $t_{m}=x_{m}^{2}/D$ being
the diffusive time scale. We have illustrated the resulting experimental plots
in Figure 3. In the experiment, $E_{b}/K_{B}T=13$, $A=0.2$ and
$t_{m}=2.52~{}\text{ms}$. With the form of $X_{1}(t)$ and $X_{2}(t)$ in Eqs.
(95) and (96), we use the formula in Eq. (21) to estimate the mean entropy
produced as a function of the cycle time.

Figure 4 illustrates the primary outcome of this section. The estimated mean
entropy production $S_{\text{tot}}$ has been plotted for different
dimensionless cycle times $t_{f}$. The results are essentially in perfect
agreement with simulation results despite the very low amount of data and the
non-trivial dynamics of the system. This shows that our method is well-suited
to study the entropy production of general time-dependent experimental
systems.

##  7 Conclusion

In this paper, we have constructed a general framework to derive a lower bound
on the mean value of the total entropy production for arbitrary time-dependent
systems described by overdamped Langevin equation. The bound, thus obtained,
is given in terms of the experimentally accessible quantities, namely the
(mixed) moments of the observable. If one only includes the first two moments,
one gets a simple analytical expression for the lower bound, whereas including
higher moments leads to a numerical scheme. We tested our results both on two
analytically solvable toy models and on experimental data for bit erasure,
taken from [67]. The lower bound is surprisingly close to the real entropy
production, even if one only takes into account the first two moments. This
makes it a perfect method to infer the entropy production and opens up the
possibility to infer entropy production in more complicated processes in
future research. One can for example think about biological systems, where
upstream processes such as glycolytic oscillations [48] or oscillatory
dynamics in transcription factors [49], can effectively lead to time-dependent
thermodynamic descriptions.

We end by stressing that, although we view thermodynamic inference as the main
application of our results in sections 3 and 4, the thermodynamic bounds Eqs.
(21) and (48) are also interesting in their own rights and might lead to a
broad range of applications similar to the thermodynamic speed limit and the
thermodynamic uncertainty relations [6, 19].

## Acknowledgement

We would like to thank Prithviraj Basak and John Bechhoefer for fruitful
discussions on the paper and for providing the data used in Section 6. We also
thank Jonas Berx for carefully reading the manuscript. This project has
received funding from the European Union’s Horizon 2020 research and
innovation program under the Marie Sklodowska-Curie grant agreement No. 847523
‘INTERACTIONS’ and grant agreement No. 101064626 ’TSBC’ and from the Novo
Nordisk Foundation (grant No. NNF18SA0035142 and NNF21OC0071284).

##  Appendix A Optimizing the entropy production rate

Throughout our paper, we have considered a method that gives us the optimal
value of total entropy produced $S_{\text{tot}}(t_{f})$ within fixed duration
$[0,t_{f}]$ based on the information about moments. Here, we showcase another
method that optimises the entropy production rate
$\sigma(t)=\dot{S}_{\text{tot}}(t)$ instead of the total dissipation
$S_{\text{tot}}(t)$. For simplicity, we focus on one dimensional case since
the extension to higher dimension is quite straightforward.

For systems obeying Langevin equation (1), the mean entropy production rate
$\sigma(t)$ is given by [2]

$\displaystyle\sigma(t)=\frac{k_{B}}{D}\int_{-\infty}^{\infty}dx~{}P(x,t)~{}v(x,t)^{2}.$
(97)

In order to optimize this with given moments

$\displaystyle X_{n}(t)$ $\displaystyle=\langle
x^{n}(t)\rangle=\int_{-\infty}^{\infty}dx~{}x^{n}P(x,t),$ (98)
$\displaystyle\dot{X}_{n}(t)$ $\displaystyle=\frac{d}{dt}\langle
x^{n}(t)\rangle=n\int_{-\infty}^{\infty}dx~{}x^{n-1}~{}P(x,t)~{}v(x,t),$ (99)

we consider the following objective function

$\displaystyle\mathbb{L}(v)$
$\displaystyle=\frac{D}{k_{B}}\sigma(t)+\sum_{i=1}^{m}\xi_{i}(t)\dot{X}_{i}(t),$
(100)
$\displaystyle=\int_{-\infty}^{\infty}dx~{}P(x,t)\left[v(x,t)^{2}+\sum_{i=1}^{m}i~{}\xi_{i}(t)~{}v(x,t)\right],$
(101)

where $\xi_{i}(t)$ for $1\leq i\leq m$ are the Lagrange multipliers
corresponding to the first $m$-moments. Now performing the minimisation
$\frac{\delta~{}\mathbb{L}(v)}{\delta v(x,t)}\Big{|}_{v^{12..m}(x,t)}=0$, we
obtain

$\displaystyle
v^{12..m}(x,t)=-\frac{1}{2}\sum_{i=1}^{m}~{}i~{}\xi_{i}(t)~{}x^{i-1},$ (102)

such that the Lagrange multipliers can be computed by inserting this solution
in Eq. (99). This gives a set of $m$ linear equations for
$\xi_{i}(t),~{}\xi_{2}(t),...,\xi_{m}(t)$ as

$\displaystyle-\frac{2}{i}\dot{X}_{i}(t)$
$\displaystyle=\xi_{1}(t)X_{i-1}(t)+2\xi_{2}(t)X_{i}(t)+...+m\xi_{m}(t)X_{m-2+i}(t),~{}~{}~{}(\text{for
}1\leq i\leq m).$ (103)

Figure 5: Comparison of two methods for free diffusion in Eq. (65) with given
first four moments. Method I minimises the rate $\sigma(t)$ in Eq. (97)
whereas Method II minimises the total dissipation $S_{\text{tot}}(t)$ as done
in Section 5.1. The two results are clearly different. Chosen parameters are
$k_{B}=1$ and $D=0.2$.

Solving these equations completely specifies $v(x,t)$ and plugging this in Eq.
(97), we find the optimal entropy production rate

$\displaystyle\sigma(t)\geq~{}\sigma^{12..m}(t)=\frac{k_{B}}{4D}\sum_{i,j=1}^{m}\xi_{i}(t)~{}\xi_{j}(t)~{}X_{i+j-2}(t).$
(104)

In this case also, we can get the closed form solutions for $m=2$

$\displaystyle v^{12}(x,t)=\lambda_{1}(t)+\lambda_{2}(t)x,$ (105)
$\displaystyle\sigma^{12}(t)=\frac{k_{B}}{D}~{}\left[\frac{\dot{A}_{2}(t)^{2}}{4A_{2}(t)}+\dot{X}_{1}(t)^{2}\right].$
(106)

where $\lambda$-functions are given in Eq. (20) and
$A_{2}(t)=X_{2}(t)-X_{1}(t)^{2}$ is the variance of $x(t)$. Both these results
match with Eqs. (21) and (22) based on our previous method (where we optimise
the total dissipation and not the rate). This means that for $m=2$, the
optimal value of $S_{\text{tot}}(t_{f})$ corresponds to optimising the rate
$\sigma(t)$ at every time instant $0\leq t\leq t_{f}$.

In the remaining part of this section, we prove that this does not remain true
for general $m$ and the two methods yield different results. To see this, let
us construct $y(x,t)$ corresponding to $v(x,t)$ in Eq. (102) by using Eq. (6)

$\displaystyle\dot{y}(x,t)$
$\displaystyle=-\frac{1}{2}\sum_{i=1}^{m}~{}i~{}\xi_{i}(t)~{}y\left(x,t\right)^{i-1}.$
(107)

Taking a time derivative, we find

$\displaystyle\ddot{y}(x,t)=-\frac{1}{2}\sum_{i=1}^{m}~{}i~{}\dot{\xi}_{i}(t)~{}y\left(x,t\right)^{i-1}+\frac{1}{4}\sum_{i,j=1}^{m}ij(i-1)~{}\xi_{i}(t)~{}\xi_{j}(t)~{}y(x,t)^{i+j-3}.$
(108)

For general $m$, this does not match with Eq. (30). For example, we get terms
like $y(x,t)^{m},~{}y(x,t)^{m+1}$ and so on in Eq. (108) which do not appear
in Eq. (30). Therefore, two methods are inequivalent for general $m>2$. This
in-equivalence has also been illustrated in Figure (5) for the free diffusion
model considered in Eq. (65) by fixing the first four moments in Eqs.
(68)-(70). In general, the method discussed in the main part of the paper will
give a tighter bound for the total dissipation.

##  Appendix B Bound $S_{\text{tot}}^{12}(t_{f})$ converges to exact
$S_{\text{tot}}(t_{f})$ for Gaussian distributions

In this appendix, we prove that the bound $S_{\text{tot}}^{12}(t_{f})$
obtained in Eq. (84) converges with the exact average dissipation for Gaussian
distributions. Let us begin with a general process $x(t)$ with time-dependent
potential $V(x,t)=\frac{K(t)}{2}x^{2}-x\mathcal{F}(t)$ where $K(t)$ and
$\mathcal{F}(t)$ are some arbitrary functions of time. This process satisfies
the Langevin equation

$\displaystyle\frac{dx}{dt}=-K(t)x+\mathcal{F}(t)+\sqrt{2D}~{}\eta(t).$ (109)

If the initial position is drawn from the Gaussian distribution, then the
process $x(t)$ also admits a Gaussian distribution given by

$\displaystyle P(x,t)=\frac{1}{\sqrt{2\pi
A_{2}(t)}}\exp\left[-\frac{\left(x-X_{1}(t)\right)^{2}}{2A_{2}(t)}\right].$
(110)

Combining this with Eq. (3), we can write the probability flux as

$\displaystyle
v(x,t)=\left(\frac{D}{A_{2}(t)}-K(t)\right)x+\left(\mathcal{F}(t)-\frac{DX_{1}(t)}{A_{2}(t)}\right).$
(111)

Since our aim is to express $v(x,t)$ completely in terms of moments and their
time derivatives, we proceed to find both $K(t)$ and $\mathcal{F}(t)$ in terms
of moments. Using Eq. (111), it is easy to show that the variance $A_{2}(t)$
turns out to be

$\displaystyle\frac{A_{2}(t)}{\mathcal{Z}^{2}(t)}=A_{2}(0)+2D\int_{0}^{t}\frac{dw}{\mathcal{Z}(w)^{2}},~{}~{}~{}~{}\text{with
}\mathcal{Z}(t)=e^{-\int^{t}_{0}K(w)~{}dw}.$ (112)

Taking time derivative on both sides, we get

$\displaystyle K(t)-\frac{D}{A_{2}(t)}=\frac{\dot{A_{2}}(t)}{2A_{2}(t)}.$
(113)

This expresses the first term in Eq. (111) in terms of the variance. Next we
look at the second term. For this, we first express the mean as

$\displaystyle\frac{X_{1}(t)}{\mathcal{Z}(t)}=X_{1}(0)+\int_{0}^{t}dw~{}\frac{\mathcal{F}(w)}{\mathcal{Z}(w)},$
(114)

and then take its time derivative to yield

$\displaystyle\mathcal{F}(t)-\frac{DX_{1}(t)}{A_{2}(t)}=\dot{X}_{1}(t)+\frac{\dot{A}_{2}(t)~{}X_{1}(t)}{2A_{2}(t)}.$
(115)

Substituting Eqs. (113) and (115) in Eq. (111), we obtain

$\displaystyle
v(x,t)=-\frac{\dot{A_{2}}(t)}{2A_{2}(t)}~{}x+\dot{X}_{1}(t)+\frac{\dot{A}_{2}(t)~{}X_{1}(t)}{2A_{2}(t)}.$
(116)

Now the expression of flux involves only first two cumulants and their time
derivatives. Finally, we use Eq. (4) to get the total entropy production as

$\displaystyle
S_{\text{tot}}\left(t_{f}\right)=\frac{k_{B}}{D}\int_{0}^{t_{f}}\left[\dot{X}_{1}(t)^{2}+\frac{\dot{A}_{2}(t)^{2}}{4A_{2}(t)}\right].$
(117)

But this is also the bound $S_{\text{tot}}^{12}(t_{f})$ derived in Eq. (84).
Hence, we conclude that $S_{\text{tot}}^{12}(t_{f})$ converges to its exact
counterpart $S_{\text{tot}}\left(t_{f}\right)$ for Gaussian distributions.

##  Appendix C Derivation of Eq. (41) in $d$ dimension

This appendix presents a derivation of the Euler-Lagrange equation (41) in $d$
dimensions. For this, we first consider $\mathbb{y}\left(\mathbb{x},t\right)$
in Eq. (40) which is mathematically equivalent to switching from Eulerian
description to the Lagrangian description in fluid mechanics [58]. Using the
equivalence of two descriptions, we have

$\displaystyle\int_{-\infty}^{\infty}d\mathbb{x}~{}P\left(\mathbb{x},t\right)~{}g\left(\mathbb{x},t\right)$
$\displaystyle=\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)~{}g\left(\mathbb{y}(\mathbb{x},t),t\right),$
(118)
$\displaystyle\int_{-\infty}^{\infty}d\mathbb{x}~{}P\left(\mathbb{x},t\right)~{}\mathbb{v}(\mathbb{x},t)~{}g\left(\mathbb{x},t\right)$
$\displaystyle=\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)~{}\dot{\mathbb{y}}(\mathbb{x},t)~{}g\left(\mathbb{y}(\mathbb{x},t),t\right).$
(119)

Plugging
$g\left(\mathbb{x},t\right)=\delta\left(\mathbb{x}-\mathbb{y}(\mathbb{x^{\prime}},t)\right)$,
we get

$\displaystyle
P\left(\mathbb{y}(\mathbb{x^{\prime}},t),t\right)=\frac{P_{0}(\mathbb{x^{\prime}})}{\text{det}\left(\mathbb{\nabla}\mathbb{y}(\mathbb{x^{\prime}},t)\right)}.$
(120)

On the other hand, putting
$g\left(\mathbb{x},t\right)=\mathbb{v}(\mathbb{x},t)^{2}$ in Eq. (118), we
obtain

$\displaystyle
S_{\text{tot}}(t_{f})=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt~{}\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}(\mathbb{x})~{}\dot{\mathbb{y}}\left(\mathbb{x},t\right)^{2}.$
(121)

In fact, one can use $\mathbb{y}\left(\mathbb{x},t\right)$ to recast the
optimisation of action
$\mathbb{S}\left(\mathbb{x},\dot{\mathbb{v}},t_{f}\right)$ in Eq. (39) as an
optimisation problem in $\mathbb{y}\left(\mathbb{x},t\right)$. To see this, we
first rewrite $\mathbb{S}\left(\mathbb{x},\dot{\mathbb{v}},t_{f}\right)$ as

$\displaystyle\mathbb{S}\left(\mathbb{y},\dot{\mathbb{y}},t_{f}\right)$
$\displaystyle=\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)\left[\dot{\mathbb{y}}\left(\mathbb{x},t\right)^{2}+\sum_{i=1}^{d}\mu_{i}(t)~{}y_{i}\left(\mathbb{x},t\right)+\sum_{i,j=1}^{d}\lambda_{ij}(t)~{}y_{i}\left(\mathbb{x},t\right)y_{j}\left(\mathbb{x},t\right)\right]$
$\displaystyle+\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)\left[\sum_{i=1}^{d}\alpha_{i}(t_{f})~{}y_{i}\left(\mathbb{x},t_{f}\right)+\sum_{i,j=1}^{d}\beta_{ij}(t_{f})~{}y_{i}\left(\mathbb{x},t_{f}\right)y_{j}\left(\mathbb{x},t_{f}\right)\right].$
(122)

For a small change in path
$\mathbb{y}\left(\mathbb{x},t\right)\to\mathbb{y}\left(\mathbb{x},t\right)+\delta\mathbb{y}\left(\mathbb{x},t\right)$,
the total change in action is given by

$\displaystyle\delta\mathbb{S}=$
$\displaystyle\int_{0}^{t_{f}}dt\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)\sum_{i=1}^{d}\left[-2\ddot{y}_{i}(\mathbb{x},t)+\mu_{i}(t)+\sum_{j=1}^{d}\tilde{\lambda}_{ij}(t)~{}y_{j}(\mathbb{x},t)\right]\delta
y_{i}(\mathbb{x},t)$
$\displaystyle+\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)~{}\sum_{i=1}^{d}\left[2\dot{y}_{i}\left(\mathbb{x},t_{f}\right)+\alpha_{i}(t_{f})+\sum_{j=1}^{d}\tilde{\beta}_{ij}(t_{f})~{}y_{j}\left(\mathbb{x},t_{f}\right)\right]\delta
y_{i}(\mathbb{x},t_{f})$
$\displaystyle-2\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)~{}\sum_{i=1}^{d}\dot{y}_{i}\left(\mathbb{x},0\right)\delta
y_{i}(\mathbb{x},0).$ (123)

where $\tilde{\lambda}_{ij}(t)=\lambda_{ij}(t)+\lambda_{ji}(t)$ and
$\tilde{\beta}_{ij}(t)=\beta_{ij}(t)+\beta_{ji}(t)$. For optimal path, this
change in action should vanish. Vanishing of the first line for arbitrary
$\delta y_{i}\left(\mathbb{x},t\right)$ gives the Euler-Lagrange equation

$\displaystyle
2\ddot{y}_{i}(\mathbb{x},t)=\mu_{i}(t)+\sum_{j=1}^{d}\tilde{\lambda}_{ij}(t)~{}y_{j}(\mathbb{x},t).$
(124)

On the other hand, we get appropriate boundary conditions by demanding
vanishing of the other two lines. Since, we begin with fixed initial condition
$\mathbb{y}(\mathbb{x},0)=\mathbb{x}$, the third line in Eq. (123) trivially
goes to zero. However, the value of $\mathbb{y}(\mathbb{x},t_{f})$ at the
final time is not fixed which means
$\delta\mathbb{y}\left(\mathbb{x},t_{f}\right)\neq 0$. Hence for the change
$\delta\mathbb{S}$ to become zero, we must have the pre-factor associated with
$\delta y_{i}\left(\mathbb{x},t_{f}\right)$ in the second line equal to zero.
Thus, we derive the following boundary conditions in time

$\displaystyle
2\dot{y}_{i}\left(\mathbb{x},t_{f}\right)+\alpha_{i}(t_{f})+\sum_{j=1}^{d}\tilde{\beta}_{ij}(t_{f})~{}y_{j}\left(\mathbb{x},t_{f}\right)=0,$
(125) $\displaystyle y_{i}(\mathbb{x},0)=x_{i}.$ (126)

Due to the linearity of Eq. (124) in $\mathbb{y}(\mathbb{x},t)$, the general
solution satisfies

$\displaystyle\dot{y}_{i}\left(\mathbb{x},t\right)=-\frac{\nu_{i}(t)}{2}-\sum_{j=1}^{d}\frac{\gamma_{ij}(t)}{2}~{}y_{j}\left(\mathbb{x},t\right),~{}~{}~{}~{}\text{for
}0\leq t\leq t_{f}.$ (127)

In order to see that this is consistent with Eq. (124), we take its time
derivative and compare the resulting $\ddot{y}_{i}\left(\mathbb{x},t\right)$
with Eq. (124). This gives rise to the following relations

$\displaystyle\mu_{i}(t)=-\dot{\nu}_{i}(t)+\frac{1}{2}\sum_{l=1}^{d}\frac{\gamma_{il}(t)~{}\nu_{l}(t)}{2},~{}~{}\text{and
}~{}~{}\tilde{\lambda}_{ij}(t)=-\dot{\gamma}_{ij}(t)+\sum_{l=1}^{d}\frac{\gamma_{il}(t)~{}\gamma_{lj}(t)}{2}.$
(128)

These relations ensure that Eqs. (124) and (127) are consistent with each
other. To summarize, we have derived the solution for the optimal path. We now
have to simply insert this solution in Eq. (121) to get the bound on total
dissipation. However, before that we need to specify $\nu(t)$ and $\gamma(t)$
functions. For this, we use the definition of moments

$\displaystyle
X_{1}^{i}(t)=\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)y_{i}\left(\mathbb{x},t\right),~{}~{}~{}\text{and
}X_{2}^{i,j}(t)=\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}\left(\mathbb{x}\right)y_{i}\left(\mathbb{x},t\right)y_{j}\left(\mathbb{x},t\right).$
(129)

and take their time derivative to yield

$\displaystyle\nu_{i}(t)$
$\displaystyle=-2\dot{X}_{1}^{i}(t)-\sum_{l=1}^{d}\gamma_{il}(t)~{}X_{1}^{l}(t),$
(130) $\displaystyle-2\dot{A}_{ij}(t)$
$\displaystyle=\sum_{l=1}^{d}\left[\gamma_{il}(t)A_{jl}(t)+\gamma_{jl}(t)A_{il}(t)\right],$
(131)

where $A_{ij}(t)=X_{2}^{i,j}(t)-X_{1}^{i}(t)X_{1}^{j}(t)$. We have quoted this
result in Eq. (42) of the main text.

###  C.1 Simplification in two dimensions

For two dimensional case, we can explicitly solve Eqs. (131) and obtain

$\displaystyle\gamma_{11}(t)$
$\displaystyle=\frac{-A_{22}(t)\dot{A}_{11}(t)\left[A_{11}(t)+A_{22}(t)\right]+2A_{12}(t)A_{22}(t)\dot{A}_{12}(t)+A_{12}(t)^{2}\left[\dot{A}_{11}(t)-\dot{A}_{22}(t)\right]}{\left[A_{11}(t)+A_{22}(t)\right]\left[A_{11}(t)A_{22}(t)-A_{12}(t)^{2}\right]},$
$\displaystyle\gamma_{22}(t)$
$\displaystyle=\frac{-A_{11}(t)\dot{A}_{22}(t)\left[A_{11}(t)+A_{22}(t)\right]+2A_{12}(t)A_{11}(t)\dot{A}_{12}(t)+A_{12}(t)^{2}\left[-\dot{A}_{11}(t)+\dot{A}_{22}(t)\right]}{\left[A_{11}(t)+A_{22}(t)\right]\left[A_{11}(t)A_{22}(t)-A_{12}(t)^{2}\right]},$
$\displaystyle\gamma_{12}(t)$
$\displaystyle=\frac{\dot{A}_{11}(t)A_{12}(t)A_{22}(t)+A_{11}(t)A_{12}(t)\dot{A}_{22}(t)-2A_{11}(t)\dot{A}_{12}(t)A_{22}(t)}{\left[A_{11}(t)+A_{22}(t)\right]\left[A_{11}(t)A_{22}(t)-A_{12}(t)^{2}\right]}.$

Plugging them in formula

$\displaystyle S_{\text{tot}}(t_{f})$
$\displaystyle=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt~{}\int_{-\infty}^{\infty}d\mathbb{x}~{}P_{0}(\mathbb{x})~{}\dot{\mathbb{y}}\left(\mathbb{x},t\right)^{2},$
(132)

we obtain the total dissipation bound as

$\displaystyle S_{\text{tot}}^{12}(t_{f})$
$\displaystyle=\frac{k_{B}}{D}\int_{0}^{t_{f}}dt\left[\dot{X}_{1}^{1}(t)^{2}+\frac{\dot{A}_{11}(t)^{2}}{4A_{11}(t)}+\dot{X}_{1}^{2}(t)^{2}+\frac{\dot{A}_{22}(t)^{2}}{4A_{22}(t)}\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+\frac{\left\\{A_{12}(t)\frac{d}{dt}\left(A_{11}(t)A_{22}(t)\right)-2A_{11}(t)A_{22}(t)\dot{A}_{12}(t)\right\\}^{2}}{4A_{11}(t)A_{22}(t)\left\\{A_{11}(t)+A_{22}(t)\right\\}\left\\{A_{11}(t)A_{22}(t)-A_{12}(t)^{2}\right\\}}\right].$
(133)

This result has been written in Eq. (49) in the main text.

##  Appendix D Experimental protocols

cycle time | m(t) | g(t) | f(t) | action
---|---|---|---|---
$t\in[0,0.25t_{f}]$ | $\left(1-\frac{t}{0.25t_{f}}\right)^{2}$ | $1$ | $0$ | symmetrize potential
$t\in[0.25t_{f},0.4t_{f}]$ | $0$ | $\left(2-\frac{t}{0.25t_{f}}\right)^{2}$ | $0$ | lower barrier
$t\in[0.4t_{f},0.6t_{f}]$ | $0$ | $\left(2-\frac{t}{0.25t_{f}}\right)^{2}$ | $\left(\frac{t/t_{f}-0.4}{0.2}\right)$ | tilt
$t\in[0.6t_{f},0.65t_{f}]$ | $0$ | $\left(2-\frac{t}{0.25t_{f}}\right)^{2}$ | $1$ | raise barrier
$t\in[0.65t_{f},0.75t_{f}]$ | $0$ | $\left(2-\frac{t}{0.25t_{f}}\right)^{2}$ | $\left(1-\frac{t/t_{f}-0.65}{0.1}\right)$ | Untilt
$t\in[0.75t_{f},t_{f}]$ | $\left(3-\frac{t}{0.25t_{f}}\right)^{2}$ | $1$ | $0$ | asymmetrize potential
Table 1: Expressions of the dimensionless functions $m(t)$, $g(t)$ and $f(t)$ that are used in order to change the double-well potential in the bit erasure experiment. Cycle time | No. of trajectories
---|---
13.36 | 40
14.89 | 53
17.30 | 40
19.79 | 28
21.19 | 20
28.59 | 22
40.78 | 16
Table 2: Number of trajectories considered for different cycle times
(rescaled) in the experiment.

## References

## References

  * [1] Seifert U 2012 Stochastic thermodynamics, fluctuation theorems and molecular machines Rep. Prog. Phys. 75 126001
  * [2] Peliti L and Pigolotti S 2021 Stochastic Thermodynamics: An Introduction (Princeton University Press)
  * [3] Sekimoto K 1998 Langevin equation and thermodynamics Prog. Th. Phys. Supp. 130 17–27
  * [4] Seifert U 2005 Entropy Production along a Stochastic Trajectory and an Integral Fluctuation Theorem Phys. Rev. Lett. 95 040602
  * [5] Hayashi K, Ueno H, Iino R and Noji H 2010 Fluctuation Theorem Applied to $\textbf{F}_{1}$-ATPase Phys. Rev. Lett. 104 218103
  * [6] Barato A C and Seifert U 2015 Thermodynamic uncertainty relation for biomolecular processes Phys. Rev. Lett. 114 158101
  * [7] Gingrich T R, Horowitz J M, Perunov N and England J L 2016 Dissipation bounds all steady-state current fluctuations Phys. Rev. Lett. 116 120601
  * [8] Proesmans K and Van den Broeck C 2017 Discrete-time thermodynamic uncertainty relation Europhys. Lett. 119 20001
  * [9] Hasegawa Y and Van Vu T 2019 Fluctuation theorem uncertainty relation Phys. Rev. Lett. 123 110602
  * [10] Timpanaro A M, Guarnieri G, Goold J and Landi G T 2019 Thermodynamic uncertainty relations from exchange fluctuation theorems Phys. Rev. Lett. 123 090604
  * [11] Koyuk T and Seifert U 2019 Operationally accessible bounds on fluctuations and entropy production in periodically driven systems Phys. Rev. Lett. 122 230601
  * [12] Proesmans K and Horowitz J M 2019 Hysteretic thermodynamic uncertainty relation for systems with broken time-reversal symmetry J. Stat. Mech.: Theory Exp. 054005
  * [13] Harunari P E, Fiore C E and Proesmans K 2020 Exact statistics and thermodynamic uncertainty relations for a periodically driven electron pump J. Phys. A: Math. Theor. 53 374001
  * [14] Pal S, Saryal S, Segal D, Mahesh T S and Agarwalla B K 2020 Experimental study of the thermodynamic uncertainty relation Phys. Rev. Res. 2 022044
  * [15] Pietzonka P and Seifert U 2018 Universal Trade-Off between Power, Efficiency, and Constancy in Steady-State Heat Engines Phys. Rev. Lett. 120 190602
  * [16] Proesmans K, Ehrich J and Bechhoefer J 2020 Finite-Time Landauer Principle Phys. Rev. Lett. 125 100602
  * [17] Ilker E, Güngör O, Kuznets-Speck B, Chiel J, Deffner S and Hinczewski M 2022 Shortcuts in Stochastic Systems and Control of Biophysical Processes Phys. Rev. X 12 021048
  * [18] Murugan A, Huse D A and Leibler S 2012 Speed, dissipation, and error in kinetic proofreading Proc. Natl. Acad. Sci. 109 12034-39
  * [19] Shiraishi N, Funo K and Saito K 2018 Speed limit for classical stochastic processes Phys. Rev. Lett. 121 070601
  * [20] Aurell E, Mejía-Monasterio C and Muratore-Ginanneschi P 2011 Optimal protocols and optimal transport in stochastic thermodynamics Phys. Rev. Lett. 106 250601
  * [21] Ito S and Dechant A 2020 Stochastic time evolution, information geometry, and the cramér-rao bound Phys. Rev. X 10 021056
  * [22] Aurell E, Gawedzki K, Mejía-Monasterio C, Mohayaee R and Muratore-Ginanneschi P 2012 Refined second law of thermodynamics for fast random processes J. Stat. Phys. 147 487–505
  * [23] Sivak D A and Crooks G E 2012 Thermodynamic metrics and optimal paths Phys. Rev. Lett. 108 190602
  * [24] Zhen Y, Egloff D, Modi K and Dahlsten O 2021 Universal bound on energy cost of bit reset in finite time Phys. Rev. Lett. 127 190602
  * [25] Proesmans K, Ehrich J and Bechhoefer J 2020 Optimal finite-time bit erasure under full control Phys. Rev. E 102 032105
  * [26] Van Vu T and Saito K 2022 Finite-Time Quantum Landauer Principle and Quantum Coherencee Phys. Rev. Lett. 128 010602
  * [27] Dechant A 2022 Minimum entropy production, detailed balance and Wasserstein distance for continuous-time Markov processes J. Phys. A: Math.Theor. 55 094001
  * [28] Falasco G and Esposito M 2020 Dissipation-time uncertainty relation Phys. Rev. Lett. 125 120604
  * [29] Neri I 2022 Universal tradeoff relation between speed, uncertainty, and dissipation in nonequilibrium stationary states SciPost Physics  12 139
  * [30] Kuznets-Speck B and Limmer,D T 2021 Dissipation bounds the amplification of transition rates far from equilibrium Proc. Natl. Acad. Sci. 118 e2020863118
  * [31] Yan L-L et al 2022 Experimental verification of dissipation-time uncertainty relation Phys. Rev. Lett. 128 050603
  * [32] Seifert U 2019 From Stochastic Thermodynamics to Thermodynamic Inference Ann. Rev. Cond. Mat. Phys. 10 171-192
  * [33] Gingrich T, Rotskoff G, and Horowitz J 2017 Inferring dissipation from current fluctuations J. Phys. A 50 184004
  * [34] Li J, Horowitz J M, Gingrich T R and Fakhri N 2019 Quantifying dissipation using fluctuating currents Nat. Commun. 10 1666
  * [35] Manikandan S K, Gupta D and Krishnamurthy S 2020 Inferring Entropy Production from Short Experiments Phys. Rev. Lett. 124 120603
  * [36] Martínez I A, Bisker G, Horowitz J M and Parrondo J M R 2019 Inferring broken detailed balance in the absence of observable currents Nat. Commun. 10 1-10
  * [37] Skinner D, Dunkel J 2021 Estimating entropy production from waiting time distributions Phys. Rev. Lett. 127 198101
  * [38] Harunari P E, Dutta A, Polettini M and Roldán E 2022 What to Learn from a Few Visible Transitions’ Statistics? Phys. Rev. X 12 041026
  * [39] Van der Meer J, Ertel B and Seifert U 2022 Thermodynamic Inference in Partially Accessible Markov Networks: A Unifying Perspective from Transition-Based Waiting Time Distributions Phys. Rev. X 12 031025
  * [40] Pietzonka P, Coghi F 2023 Thermodynamic cost for precision of general counting observables arXiv 2305.15392
  * [41] Lynn C, Holmes C, Bialek W, Schwab D 2022 Decomposing the local arrow of time in interacting systems Phys. Rev. Lett. 129 118101
  * [42] Lynn C, Holmes C, Bialek W, Schwab D 2022 Emergence of local irreversibility in complex interacting systems Phys. Rev. E 106 034102
  * [43] Roldán E and Parrondo J M R 2010 Estimating dissipation from single stationary trajectories Phys. Rev. Lett. 105 150607
  * [44] Lander B, Mehl J, Blickle V, Bechinger C and Seifert U 2012 Noninvasive measurement of dissipation in colloidal systems Phys. Rev. E. 86 030401
  * [45] Otsubo S, Ito S, Dechant A and Sagawa T 2020 Estimating entropy production by machine learning of short-time fluctuating currents Phys. Rev. E. 101 062106
  * [46] Van Vu T, Tuan Vo V and Hasegawa Y 2020 Entropy production estimation with optimal current Phys. Rev. E. 101 042138
  * [47] Kim D-K, Bae Y, Lee S and Jeong H 2020 Learning Entropy Production via Neural Networks Phys. Rev. Lett. 125 140604
  * [48] Chandra F, Buzi G and Doyle J 2011 Glycolytic oscillations and limits on robust efficiency Science 333 187
  * [49] Heltberg M, Krishna S and Jensen M 2019 On chaotic dynamics in transcription factors and the associated effects in differential gene regulation Nat. Comm. 10 71
  * [50] Jun Y, Gavrilov M and Bechhoefer J 2014 High-Precision Test of Landauer’s Principle in a Feedback Trap Phys. Rev. Lett. 113 190601
  * [51] Freitas N, Delvenne J-C, and Esposito M 2020 Stochastic and Quantum Thermodynamics of Driven RLC Networks Phys. Rev. X 10 031005
  * [52] Koyuk T and Seifert U 2020 Thermodynamic uncertainty relation for time-dependent driving Phys. Rev. Lett. 125 260604
  * [53] Dechant A and Sakurai Y 2023 Thermodynamic interpretation of Wasserstein distance arXiv:1912.08405
  * [54] Otsubo S, Manikandan S K, Sagawa T and Krishnamurthy S 2022 Estimating time-dependent entropy production from non-equilibrium trajectories Commun. Phys. 5 11
  * [55] Lee S, Kim D-K, Park J-M, Kim W K, Park H and Lee J S 2023 Multidimensional entropic bound: Estimator of entropy production for Langevin dynamics with an arbitrary time-dependent protocol Phys. Rev. Research 5 013194
  * [56] Blickle V, Speck T, Helden L, Seifert U and Bechinger, C 2006 Thermodynamics of a Colloidal Particle in a Time-Dependent Nonharmonic Potential Phys. Rev. Lett. 96 070603
  * [57] Garanin D A 1997 Fokker-Planck and Landau-Lifshitz-Bloch equations for classical ferromagnets Phys. Rev. B 55 3050
  * [58] Benamou J-D and Brenier Y 2000 A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem Numer. Math. 84 375-393
  * [59] Villani C 2003 Topics in Optimal Transportation (American Mathematical Society, Rhode Island)
  * [60] Kopp R 1962 Pontryagin maximum principle Mathematics in Science and Engineering 5 255-279
  * [61] Batchelor C K and Batchelor G 2000 An introduction to fluid dynamics (Cambridge University Press)
  * [62] Proesmans K 2023 Precision-dissipation trade-off for driven stochastic systems Commun. Phys. 6 226
  * [63] Landauer R 1961 Irreversibility and heat generation in the computing process IBM J. Res. Develop. 5 183-191
  * [64] Bérut A, Arakelyan A, Petrosyan A, Ciliberto S, Dillenschneider R and Lutz E 2012 Experimental verification of Landauer’s principle linking information and thermodynamics Nature 483 187-189
  * [65] Hong J, Lambson B, Dhuey S and Bokor J 2016 Experimental test of Landauer’s principle in single-bit operations on nanomagnetic memory bits Sci. Adv. 2 e1501492
  * [66] Martini L et al 2016 Experimental and theoretical analysis of Landauer erasure in nanomagnetic switches of different sizes Nano Energy 19 108-116
  * [67] Gavrilov M and Bechhoefer J 2016 Erasure without Work in an Asymmetric Double-Well Potential Phys. Rev. Lett. 117 200601
  * [68] Gavrilov M, Chétrite R and Bechhoefer J 2017 Direct measurement of weakly nonequilibrium system entropy is consistent with Gibbs–Shannon form Proc. Natl. Acad. Sci. 114 11097-11102


# Redactable and Sanitizable Signature Schemes: Applications and Limitations
for use in Decentralized Digital Identity Systems

Bryan Kumara  The Alan Turing Institute, London, U.K. Email: bkumara, mhooper,
cmaple, thobson, jcrowcroft @turing.ac.uk  Mark Hooper  The Alan Turing
Institute, London, U.K. Email: bkumara, mhooper, cmaple, thobson, jcrowcroft
@turing.ac.uk  Carsten Maple  The Alan Turing Institute, London, U.K. Email:
bkumara, mhooper, cmaple, thobson, jcrowcroft @turing.ac.uk  Timothy Hobson
The Alan Turing Institute, London, U.K. Email: bkumara, mhooper, cmaple,
thobson, jcrowcroft @turing.ac.uk  Jon Crowcroft  The Alan Turing Institute,
London, U.K. Email: bkumara, mhooper, cmaple, thobson, jcrowcroft
@turing.ac.uk

(18 April 2023)

###### Abstract

Redactable signature schemes and sanitizable signature schemes are methods
that permit modification of a given digital message and retain a valid
signature. This can be applied to decentralized identity systems for
delegating identity issuance and redacting sensitive information for privacy-
preserving verification of identity. We propose implementing these protocols
on a digital credential and compare them against other privacy-enhancing
techniques to assess their suitability.

##  1 Introduction

A critical aspect of economic growth is the integration of the internet as a
medium for innovation [11]. It is pivotal that that digital identities, which
refer to information used by computers to represent external agents, address
stakeholder concerns [11] [4] [2]. Every digital identity system comes with
its own set of trade-offs that considers security, reliability, resilience,
accessibility, and scalability amongst many other factors [12].
A key issue is the degree on which a single party may exert control over the
operation of the system. Centralized digital identity systems usually contain
a dedicated registration process, unique user identifiers, and the storage of
information in a central database [7]. However, the centralization of
authority and storage carries the risk of abuse such as exclusion and
censorship alongside individual profiling [14].
In contrast, decentralized identity systems typically approach the issue using
transparent, open networks and protocols [15]. They usually involve voluntary
participation and attempt to minimize data collection or formal enrolment. In
this approach, standards such as Decentralized Identifiers and Verifiable
Credentials help prove set membership to establish that individuals possess
certain attributes without needing a central authority [3].
This extended abstract examines how Redactable and Sanitizable Signature
Schemes can enable privacy preserving computation in the context of
decentralized identity systems and consider mitigating solutions to some of
their limitations. Redactable and Sanitizable Signature Scheme can be applied
to anonymous disclosure and delegating identity issuance respectively. We
propose implementing these techniques for a decentralized identity system to
better understand their performance against other privacy-enhancing
techniques.

##  2 Background

Digital signature allow users to verify the authenticity of a message [10]. A
valid signature for a given message proves the message was sent by a known
user and has not been altered during transit [10]. Traditional digital
signatures become invalid if even a single part of the associated message is
modified [10]. However, in some cases users may want to allow a third party to
alter a signed message and still retain a valid signature such as removing
personal details of patient data for a medical study [1]. To support this,
Sanitizable Signature Schemes (SSS) allow a designated third party, known as
the sanitizer, to change permissible parts of a signed message. Alternatively,
Redactable Signature Schemes (RSS) allow parties to redact parts of a signed
message [1]. RSS and SSS are not recent, however research into these
constructions are relatively limited compared to other privacy-enhancing
techniques such as multiparty computation and zero-knowledge proofs. For the
multitude constructions of RSS and SSS, readers should refer to [1]. In
general, these schemes are based on pairing-based cryptography which maps
elements between pairs of cryptographic groups. RSS and SSS gather the altered
parts of a message and commit them cryptographically, generating a derived
signature from the original signature. A verification process would confirm
that this derived signature is well-formed.
Both RSS and SSS share standard security properties including being
unforgeable, the prevention of deriving information of sanitized or redacted
parts, the immutability of inadmissible parts of a message, and accountability
of who modifies a message. Various constructions offering different properties
can be found in [1].

##  3 Applications to Decentralized Digital Identity

Though RSS and SSS have a broad range of applications, they can be
particularly useful in decentralized identity systems. The former can be
applied for Verifiable Credentials (VC) which is an open standard for digital
credentials that can represent information such as bank account ownership and
licenses by enabling selective disclosure [8]. For instance, the issuer of a
bank account ownership would be a bank, and this information can be held by
the bank account owner who becomes the holder of a credential. An associated
digital signature attests to this issuance. The holder manages the credentials
and can present it to verifiers such as mortgage companies as needed.
Verifiers can request the VC from the holders as seen in the figure below
[18]. A key step is to check if the signature of the issuer is valid, which
indicates message authenticity and integrity. However, holders may not want to
disclose the entire VC content on every request, as it may contain sensitive
information depending on the format of the digital identity.

Figure 1: DID Model proposed by W3C [18]

One option to limit VC disclosure would be using zero-knowledge proofs to
allow a party to show that a given statement is true, such as being in
possession of a certain amount of cash, while not revealing any additional
information beyond the fact that the said condition is met [13]. However,
zero-knowledge proofs are complicated techniques with multiple variants. Some
require trusted setups which may not be suitable for decentralized identity
systems with various representations and proof techniques [16]. RSS allows a
holder to redact segments of a VC but leave other components unredacted for
verification. Current theoretical performance speed of RSS indicates that this
is the more efficient solution due to their speed and smaller signature size,
while proof of knowledge protocols may have sizes larger than the signature of
an unredacted VC.
From a credential issuer’s perspective, decentralization may require some form
of delegation for scalability. Decentralized Identifiers (DIDs) are a type of
unique identifier that is persistent, verifiable, and without the need of a
centralized registry which enables a new model of decentralized digital
identity. The World Wide Web Consortium (W3C), the main international
standards organization for the World Wide Web, maintains an open standard for
DIDs. Often, DID issuers will want to delegate authority by allowing other
nodes to issue DIDs with the delegator’s approval. However, to avoid malicious
behaviour from delegates, such generating dud DIDs which may lead to Sybil
attacks, there needs to be accountability. One approach is to use a multi-
signature scheme which requires the signature of multiple parties to jointly
sign a proposed DID [9]. Using this method the delegate node would gather the
necessary information but cannot immediately issue a DID as it requires a
superior’s signature before being issued.
However, one weakness of this approach is the back-and-forth communication
required between a delegate and the delegating node. Every proposed DID needs
to be checked by the delegating node that it has been formatted correctly
before obtaining their signature to create the DID. In areas of low
connectivity and or down-time of either party, there will be a delay in
issuing DIDs from a delegate node. SSS allows a delegating node to give a
delegate blank templates that can be modified and some SSS allow sanitization
in a constrained fashion, this will ensure that a potential DID will follow
the correct format. If the downstream entity formats the DID correctly, a
valid signature can be automatically derived for the DID without another round
of communication. Like RSS, the performance time of SSS is reasonably quick
and is helped by the relatively small number fields in a W3C DID standards. Of
course, SSS can also be applied to other hierarchical models where delegation
of authority is needed, such as in Public Key Infrastructure, where root
certificate authorities delegate authorities to subordinate certificate
authorities.

##  4 Limitations and Mitigations

SSS and RSS possess limitations that need to be addressed. Using SSS, a
malicious delegate may generate multiple dud instances that are correctly
filled. Neither SSS nor multi-signatures address the veracity of the input,
and this issue of a single-source oracle for decentralized system is a problem
outside the scope of this abstract. However, the potential to generate
multiple dud instances is somewhat limited in a multi-signature scheme. If the
issuer receives an unusually high number of requests to approve DIDs, it may
indicate malicious activity and could refuse the multi-signature request. SSS
does not have traffic monitoring as sanitization requires no communication
between the two parties outside of the original issuance of sanitizable
message. One solution is to use threshold cryptography, where each DID a
delegate creates reveals a share of their secret key. Crossing the threshold
by generating too many DIDs remove the autonomy of the delegate. Though this
approach seems promising, depending on the threshold it can incur a high
computational cost [17].

Another general issue of the two is the lack of post-quantum research into the
field. Whilst other techniques such as zero-knowledge proofs and multi-party
computations have multiple instances and analysis into post-quantum variants,
there is a distinct lack of them for RSS and SSS. Only very recently, RSS has
had a lattice-based [6] and a hash-based construction [5] proposed; more
research will be needed to understand security and performance implications.
SSS has no post-quantum variant that the author is aware of, however the use
of the post-quantum SPHINCS+ scheme in RSS [5] is promising and may be used in
constructions of SSS that uses chameleon hashes to permit sanitization.
Further work is required to examine SSS and RSS benefits for decentralized
identity systems. We aim to use RSS and SSS on DIDs based on the W3C standard
to examine the size and speed of these schemes compared to alternatives such
as zero-knowledge proofs and multi-signatures.

##  5 Conclusion

Decentralized identity systems can implement selective disclosure and delegate
identity issuance using RSS and SSS respectively. These related schemes permit
the modification of a message whilst retaining a valid signature, providing an
alternative to techniques such as zero-knowledge proofs and threshold
cryptography. However, a lack of research into the two needs to be remedied
before widespread adoption is possible. We propose future work examining the
performance of these schemes against other privacy-enhancing techniques.

##  6 Acknowledgments

This work was supported, in whole or in part, by the Bill and Melinda Gates
Foundation [INV-001309]. Under the grant conditions of the Foundation, a
Creative Commons Attribution 4.0 Generic License has already been assigned to
the Author Accepted Manuscript version that might arise from this submission.

## References

  * [1] Kai Samelin Arne Bilzhause  “Position Paper: The Past, Present, and Future of Sanitizable and Redactable Signatures”  In _International Conference on Availability, Reliability and Security_ Association for Computing Machinery, 2017
  * [2] Joseph J. Atick  “Digital identity: The Essential Guide”  ID4Africa, 2013
  * [3] Zoltan Andras Lux Dirk Thatmann Sebastian Zickau Felix Beierle  “Distributed-Ledger-based Authentication with Decentralized Identifiers and Verifiable Credentials”  In _2nd Conference on Blockchain Research and Applications for Innovative Networks and Services_ IEEE, 2020
  * [4] J.L. Camp  “Digital Identity”  In _IEEE Technology and Society_ 24 IEEE, 2004
  * [5] Fei Zhu Xun Yi Alsharif Abuadbba Junwei Lup Surya Nepal Xinyi Huang  “Efficient Hash-Based Redactable Signature for Smart Grid Applications”  In _European Symposium on Research in Computer Security_ 27 Association for Computing Machinery, 2022
  * [6] Yong Zhao Shaojun Yang Wei Wu Xinyi Huang  “A Lattice-Based Redactable Signature Scheme using Cryptographic Accumulators for Trees”  In _The Computer Journal_ Oxford University Press, 2022
  * [7] Chen Yan Ji Fang  “Centralized Identity Authentication Research Based on Management Application Platform”  In _International Conference on Information Science and Engineering_ 1 IEEE, 2009
  * [8] Alexander Rieger Johannes Sedlmeir  “Digital Identities and Verifiable Credentials”  In _Business and Information Systems Engineering_ 63 Springer, 2021
  * [9] Yannick Seurin Jonas Nick  “MuSig2: Simple Two-Round Schnorr Multi-Signatures”  In _CRYPTO 2021_ IEEE, 2021
  * [10] Ravneet Kaur  “Digital Signature”  In _International Conference on Computing Sciences_ IEEE, 2012
  * [11] Ali M Al-Khouri  “Digital identity: Transforming GCC economies”  In _Innovation: Organization and Management_ 16.2 Association for Computing Machinery, 2013
  * [12] Samia Bouzefrane Maryline Laurent  “Digital Identity Management”  Elsevier, 2015
  * [13] Ueli Maurer  “Unifying Zero-Knowledge Proofs of Knowledge”  In _Progress in Cryptology - AFRIACRYPT_ Springer, 2009
  * [14] Bertjan Doose Naomi Ellemers  “Self and Social Identity”  In _Annual Review of Psychology_ 53 Annual Reviews Inc., 2002, pp. 161–186
  * [15] Khalifa Toumi Omar Dib  “Decentralized Identity Systems: Architecture, Challenges, Solutions and Future Directions”  In _Annals of Emerging Technologies in Computing_ 4.5 nternational Association of EducatorsResearchers, 2020
  * [16] Xiaoqiang Sun F. Yu Peng Zhang Zhiwei Sun Weixin Xie Xiang Peng  “A Survey on Zero-Knowledge Proof in Blockchain”  In _IEEE Network_ 35.4 IEEE, 2021
  * [17] Amandine Jambert Sebastien Canard  “On Extended Sanitizable Signature Schemes”  In _Topics in Cryptology - CT - RSA_ Springer, 2010
  * [18] Young B. Yunhee Kang  “An Empirical Study of a Trustworthy Cloud Common Data Model Using Decentralized Identifiers”  In _The Development and Application of Fuzzy Logic_ 11 MDPI, 2022