# Optimizing the Energetics of the Finite-time Driving of Field Theories

Atul Tanaji Mohite<EMAIL_ADDRESS>Department of Theoretical
Physics and Center for Biophysics, Saarland University, Saarbrücken, Germany

###### Abstract

The phase transitions for many-body systems have been understood using field
theories. A few canonical physical model classes encapsulate the underlying
physical properties of a large number of systems. The finite-time driving of
such systems and associated optimal energetic costs have not been investigated
yet. We consider two universality classes Model A and Model B, that describe
the dynamics for the non-conserved and conserved scalar order parameters
respectively. Here, using the recent developments in stochastic thermodynamics
and optimal transport theory, we analytically compute the optimal driving
protocols by minimizing the mean stochastic work required for finite-time
driving. Further, we numerically optimize the mean and variance of the
stochastic work simultaneously. Such a multi-objective optimization is called
a Pareto optimization problem and its optimal solution is a Pareto front. We
discover a first-order Pareto phase transition in the Pareto front.
Physically, it corresponds to the coexistence of two classes of optimal
driving protocols analogous to the liquid-gas coexistence for the equilibrium
phase transition. Our framework sheds light on the finite-time optimal driving
of the fields and the trade-off between the mean and fluctuations of the
optimal work.

Introduction \- Understanding finite-time optimal processes is a frontier of
non-equilibrium statistical physics. Optimal processes are related to the
Thermodynamic length, a metric distance between thermodynamic states from an
Information Geometrical viewpoint [1]. In the last decade, information
geometric interpretation of the thermodynamics emerged [2, 3, 4, 5, 6, 7, 8,
8, 9], revealing the connection between optimal paths for physical processes
and the underlying thermodynamic metric structure [10]. It was deployed on,
driven barrier crossing [11], finite time erasure of a classical bit [12],
driven open quantum system [13], optimum non-equilibrium control of
nanomagnetic spin systems [14] and other physical problems [15, 16, 17, 18,
19, 20]. This framework holds true in a linear response regime close to the
equilibrium and single particle processes, thus limiting its application on
many-body far from equilibrium collective phenomena in biophysics.

A few canonical models termed ‘universality classes’ are characterized by
similar critical properties that reveal the underlying identical symmetries of
apparently different many-body physical systems. A ‘phase’ distinguishes the
collective physical property of an interacting many-body system. The phase
transition is the change of the universality class, which is accompanied by
critical phenomena, a divergence or anomaly in the physical properties of the
system [21]. Phase transitions are ubiquitous in nature and have been well
understood using field theoretical techniques which are derived using the
underlying phenomenological physical symmetries [22, 21, 23, 24]. Those are
justified due to the existence of a systematic coarse-graining from
fluctuating microscopic dynamics of an interacting many-body system to
fluctuating field dynamics [25, 26, 27]. Even though quasistatic phase
transitions are well understood, an important question has never been
addressed, how to energetically optimize a finite time driving of phase and
what is the corresponding optimal process/protocol?

Figure 1:  The free energy landscape is denoted by the solid black line.
Changing control parameters exchanges the stability of the landscape minima,
thus changing the phase or state of the system. In a small noise limit, the
initial and final Gaussian probability distribution around the stable minima
is denoted by the dotted red line.

The framework of stochastic thermodynamics (ST) [28, 29, 30, 31] and finite-
time thermodynamics [32, 33] deals with the average of fluctuating
thermodynamic quantities, which are observed for mesoscopic systems undergoing
finite time processes. In ST, a non-vanishing entropy production rate (EPR)
quantifies the excessive energy cost for a non-equilibrium process (i.e. a
finite time process here). Optimal transport theory (OTT) as a tool to solve
nonlinear differential equations is used in Maths [34, 35, 36] and it’s
connection to ST is derived [37, 38]. It relies on the equivalence between the
EPR for finite time processes and an information-geometric distance measure
called ‘Wasserstein distance’ ($\mathcal{W}_{2}$), thus leading to a
refinement of the second law of thermodynamics [39, 40]. This framework was
implemented on multiple physical problems with a motif of finding the finite-
time optimal driving between the states. For instance, a single particle
Langevin dynamics [37, 38], consistent with [41, 42], a many-body system [40],
a thermodynamic equilibrium state [43], stochastic thermodynamic machines
[44], finite time bit erasure [45, 46, 47], optimization of stochastic heat
engines [9, 48, 49, 50, 51].

Here, we propose a framework to minimize irreversible work to drive in a
finite time an interacting many-body system described by its Ginzburg-Landau
(GL) functional from a state in one phase to a state in the other phase. The
driving dynamics are enforced by generic protocols in the probability
distribution space under the assumption of full control of the driven
dynamics. To compute EPR minimizing optimal driving protocols, we utilize the
analogy between ST [28, 29, 30, 31] and OTT [34, 35, 36]. Further, we optimize
a linear combination of the mean and fluctuations of work by modelling it as a
multi-objective optimization problem [52]. It exhibits a first-order phase
transition in the Pareto front (PF). The PF is defined as a limiting manifold
beyond which multiple quantities (mean and variance of work here) can not be
simultaneously optimized [53, 54].

Optimizing mean work for a change of phase: We consider here a standard GL
functional describing the phase for an order parameter $\phi$:

$\mathcal{H}[\phi]=\int\bigg{[}\frac{a}{2}\phi^{2}-\frac{c}{3}\phi^{3}+\frac{b}{4}\phi^{4}+\frac{\kappa}{2}(\nabla\phi)^{2}\bigg{]}d^{d}\mathbf{r}.$
(1)

$\\{a,b,c,\kappa\\}$ are phenomenological control parameters determining the
phase of the system [55]. A Langevin equation for the dynamics of $\phi$ is
given by eq. 2, so that $\beta=\\{0,1\\}$ corresponds to model A (non-
conserved $\phi$) and model B (conserved $\phi$) dynamics respectively [21].

$\dot{\phi}(\mathbf{r},t)=-\mu(i\nabla)^{2\beta}\bigg{(}\frac{\delta\mathcal{H}}{\delta\phi}\bigg{)}+\Lambda(\mathbf{r},t),$
(2)

where $\mu$ is the mobility, $T$ is the temperature of heat bath and
$\Lambda(\mathbf{r},t)$ is the white Gaussian noise satisfying the Einstein
relation
$\langle\Lambda(\mathbf{r},t)\Lambda(\mathbf{r}^{\prime},t^{\prime})\rangle=2\mu
T(i\nabla)^{2\beta}\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(t-t^{\prime})$.
The probability distribution $\mathcal{P}([\phi],t)$ for $\phi$ is obtained by
converting the Langevin equation to the Fokker-Planck equation. Hence-fourth
we switch to Fourier space, $\phi_{q}$ denotes the $q^{th}$ Fourier mode of
$\phi$.

Assuming a weak noise limit, we compute initial and final equilibrium Gaussian
probability distributions $\mathcal{P}_{i}([\phi_{q}])$ and
$\mathcal{P}_{f}([\phi_{q}])$ as a function of initial and final control
parameters of GL functional [55]. The small noise limit ensures that the
dynamics of each Fourier mode $\phi_{q}$ follows an independent Langevin
equation with $\mathcal{P}_{i}([\phi_{q}])$ and $\mathcal{P}_{f}([\phi_{q}])$
being Gaussian. Using ST, the stochastic work for driving $q^{th}$ Fourier
mode in a finite driving time $\tau$ is defined as [31]:

$W_{q}=\Delta_{0}K\frac{\phi_{q}(0)^{2}}{2}+\Delta_{\tau}K\frac{\phi_{q}(\tau)^{2}}{2}-\int_{0}^{\tau}\partial_{t}K(q,t)\frac{\phi_{q}^{2}(t)}{2},$
(3)

where $K(q,t)$ is the driving protocol in the Langevin equation for $q^{th}$
Fourier mode:

$\partial_{t}\phi_{q}=\mu(q)\Big{(}K(q,t)\mathbf{\phi}_{q}+a(q,t)\Big{)}+\Lambda_{q},$
(4)

where $\Lambda_{q}$ is the Gaussian white noise in Fourier space satisfying
$\langle\Lambda_{q}(t)\Lambda_{q^{\prime}}(t^{\prime})\rangle=2\mu
T(-qq^{\prime})^{\beta}\delta(q+q^{\prime})\delta(t-t^{\prime})$. $\Delta
K_{0}=K_{i}-K_{0}^{+}$ and $\Delta K_{\tau}=K_{\tau}^{-}-K_{f}$ are energy
cost associated with the discontinuity of $K(q,t)$ at $t=0$ and $t=\tau$
respectively. $K_{i}(q)=-(\alpha_{i}+\kappa q^{2})$ and
$K_{f}(q)=-(\alpha_{f}+\kappa q^{2})$ are the initial and final trap stiffness
obtained by linearizing eq. 2 [55].

Summing eq. 64 over all noise realization of the trajectory leads to an
expression for the mean work $\langle W\rangle=\Delta\mathcal{F}+T\Sigma$.
Here $\langle.\rangle$ denotes the noise average.
$\Delta\mathcal{F}=\ln(\mathcal{Z}_{i}/\mathcal{Z}_{f})$ is the free energy
difference between the initial and final state, where
$\mathcal{Z}=\int\mathcal{P}([\phi_{q}])d^{d}\phi_{q}$ is the partition
function. $\Delta\mathcal{F}$ is the quasi-static contribution to work.
$\Sigma$ (EPR) is the mean irreversible work, a finite time contribution to
the mean stochastic work.

For the driving dynamics of $\phi$, we map the work optimization problem onto
OTT using previously established techniques [40, 37, 38, 43]. A lower bound
for the mean irreversible work (EPR) in terms of $\mathcal{W}_{2}$ reads:

$\Sigma\geq\frac{\mathcal{W}_{2}\big{(}\mathcal{P}_{i}[\phi_{q}],\mathcal{P}_{f}[\phi_{q}]\big{)}^{2}}{\tau\>\mu(q)\>T}.$
(5)

In eq. 5, the lower bound saturates when $\mathcal{P}_{i}[\phi_{q}]$ and
$\mathcal{P}_{f}[\phi_{q}]$ are Gaussian which is ensured by the small noise
limit. OTT gives the Optimal Transport Map (OTM) as a solution for time-
dependent interpolation of the mean and the correlation function from the
initial to the final state. The OTM for each independent Fourier mode $q$ is:

$\begin{split}\phi^{*}(q,t)&=s\bar{\phi}_{f}\>\delta_{q,0}+(1-s)\bar{\phi}_{i}\>\delta_{q,0},\\\
C^{*}(q,t)&=\Big{(}s\sqrt{C_{f}(q)}+(1-s)\sqrt{C_{i}(q)}\Big{)}^{2},\end{split}$
(6)

where $s=t/\tau$ is re-scaled time. $\bar{\phi}_{i}$ and $\bar{\phi}_{f}$ are
mean of the initial and final probability distribution.
$\textit{C}_{i}(q)=1/(\alpha_{i}+\kappa q^{2})$ and
$\textit{C}_{f}(q)=1/(\alpha_{f}+\kappa q^{2})$ are initial and final
correlation functions of $q^{th}$ Fourier mode. The OTM ensures that the
probability distribution at an intermediate time is Gaussian if
$\mathcal{P}_{i}[\phi_{q}]$ and $\mathcal{P}_{f}[\phi_{q}]$ are Gaussian.
Hence in eq. 79, the linear ansatz for the driving protocol is justified.
Driven dynamics of $\phi_{q}$ in eq. 79 should be consistent with OTM. The
consistency condition yields a solution for the optimal driving protocol
$K^{*}(q,t)$ [55]:

$\begin{split}K^{*}(q,t)=\frac{\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}}{\mu(q)\>\tau\>\sqrt{C^{*}(q,t)}}-\frac{T}{C^{*}(q,t)},\end{split}$
(7)

where $C_{f}(q)\leq C_{i}(q)$ ensures the negativity of $K^{*}(q,t)$. For
$C_{f}(q)\geq C_{i}(q)$, $\exists\>\tau_{c}$ for
$\\{K^{*}(q,t)<0,\>\tau\geq\tau_{c}\\}$, which corresponds to a driving limit
for the change of phase. For $\tau\to\infty$, $C^{*}(q,t)=T/K^{*}(q,t)$
recovers the initial and final equilibrium correlations. The finite value of
$\tau$ gives the $\mu(q)$ dependent discontinuity in $K^{*}(q,t)$ at $t=0$ and
$t=\tau$, which generalizes Ref. [41] to finite-time driving of the field
dynamics. The first term in eq. 7 is the discontinuity in $K^{*}(q,t)$ due to
the finite time driving which depends on $\mu$ and $q$. The optimal advective
term $a^{*}(q,t)=\delta_{q,0}(\bar{\phi}_{f}-\bar{\phi}_{i})/(\tau\mu(q))$
vanishes for $q\neq 0$, thus the advective energetic cost is attributed to
changing the mean of $\phi$ only. $\langle W\rangle$ is simplified by plugging
eq. 7 in eq. 64 [55]:

$\begin{split}\langle
W\rangle&=\sum_{q}\Bigg{\\{}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\mu(q)\tau}+\frac{\big{(}\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}\big{)}^{2}}{\mu(q)\tau}\\\
&\hskip
56.9055pt+\frac{T}{2}\log\bigg{(}\frac{C_{i}(q)}{C_{f}(q)}\bigg{)}\Bigg{\\}}\end{split}$
(8)

The optimum mean work in eq. 8 has the structure $\langle
W\rangle=\Delta\mathcal{F}+T\Sigma$, thus satisfying the finite time second
law of thermodynamics [39]. The first two terms scale inversely with $\tau$ so
that
$T\Sigma=\mathcal{W}_{2}(\mathcal{P}_{i}([\phi_{q}])\mathcal{P}_{f}([\phi_{q}]))^{2}/(\tau\mu(q))$
holds, recovering the equivalence utilized between ST and OTT. The third term
is the quasistatic contribution $\Delta\mathcal{F}$, which can be proven by
calculating $\mathcal{Z}$. We compute closed-form analytical expressions for
$\langle W\rangle$ by integrating out $q$ [55]. The set of eqs. 5, 7 and 8
sheds light on an important distinction between conserved and non-conserved
dynamics. In particular, the frequency dependence of $\mu(q)$ for the
conserved (Model B) dynamics implies that more work is needed to drive the
low-frequency modes. Physically, it is the consequence of the conservation law
imposed for the order parameter using the Model B dynamics.

Figure 2: (a) For model A, the numerically computed irreversible work (points)
and analytical expression (solid line) are plotted for different temperatures
as a function of driving time $\tau$. (b) The cost function
$\mathcal{J}^{\lambda}_{0}$ for the $0^{th}$ Fourier mode given by eq. 9 is
plotted, $\lambda$ is varied from $1$ to $0$. The inset shows at higher values
of $\mu$, a sudden jump in $\mathcal{J}^{\lambda}_{0}$ is observed. (c) Cost
function $\mathcal{J}_{\lambda}^{0}$ as a function of $\lambda$ for forward
(blue line) and reverse (orange line) change of control parameter $\lambda$
(inset). Different slopes of the $\mathcal{J}^{\lambda}_{0}$ correspond to two
minima. The minina exchange global stability at a critical value of
$\lambda=\lambda_{c}$. Hysteresis loop for
${d\mathcal{J}^{\lambda}_{0}}/{d\lambda}$ characterises first order phase
transition. (d) The phase diagram for a Pareto optimization of single Fourier
mode. For $\mu<\mu_{c}$, $\mathcal{J}^{\lambda}_{0}$ has a single minimum,
hence no phase transition. $\mu>\mu_{c}$, two minima corresponding to the mean
and fluctuation-dominated protocols are separated by a first-order phase
transition at $\lambda=\lambda_{c}$. (e) Pareto front plotted in mean $\langle
W_{0}\rangle$ and SD $\sigma_{0}$ plane has a discontinuity at a critical
value of $\lambda$ denoted by $\lambda_{c}$. (f-i) Pareto front for multiple
Fourier modes using the cost function eq. 11 plotted in $\langle
W\rangle-\sigma^{2}$ plane (bottom left) exhibits a discontinuity which is
attributed to a first-order phase transition observed for the low-frequency
modes (top right $\&$ left). On the contrary higher frequency mode does not
exhibit a phase transition (lower left). See for the control parameters used
in the plots [55].

We further determined the work numerically by simulating the driven dynamics
in eq. 79 with $K^{*}(q,t)$ using the spectral method on a two-dimensional
discrete lattice. For model A, the analytical expression for work is in
agreement with the numerical work for different initial and final control
parameters of the GL functional fig. 2(a). We observe deviation from the
analytical expression at higher temperatures which is attributed to deviations
of $\mathcal{P}_{i}([\phi_{q}])$ and $\mathcal{P}_{f}([\phi_{q}])$ from
Gaussianity. For model B, we compute $\langle W\rangle$ in eq. 8 by a
summation of the discrete Fourier modes to avoid the infrared divergence of
$\mathcal{W}_{2}$ on the RHS.

Pareto optimization of the mean and standard deviation of work: The framework
of OTT yields exact results for optimizing the mean work. Two work
distributions with the same mean but different variances (fluctuations) lead
to the same solution of the optimization problem because the work fluctuations
are not neglected. The framework of the Pareto optimization problem (POP)
deals with the multi-objective optimization problem, i.e. simultaneously
optimizing multiple physical quantities [56]. Multi-objective optimization has
been connected to different classes of non-equilibrium phase transitions [53,
54]. Optimization of more than one quantity is observed in evolutionary trade-
offs [57, 58], vascular flow networks [59, 60], spatial complex networks [61],
multi-layer transport networks [62] and betting strategies [63].

We formulate a POP, to simultaneously optimize the mean and standard deviation
(SD) of the work for a Fourier mode ($q$). We minimize the cost function
$\mathcal{J}^{\lambda}_{q}$ in eq. 9,

$\mathcal{J}^{\lambda}_{q}=\lambda\langle W_{q}\rangle+(1-\lambda)\sigma_{q},$
(9)

where $\mathcal{J}^{\lambda}_{q}$ is a linear combination of the mean work
$\langle W_{q}\rangle$ and the SD of the work $\sigma_{q}=\sqrt{\langle
W_{q}^{2}\rangle-\langle W_{q}\rangle^{2}}$. $\lambda$ characterises the
relative weight for the mean and the SD of the work. $\lambda=1$ minimizes
average work and $\lambda=0$ work fluctuations. For a single Fourier mode, POP
for B with mobility $\mu^{\prime}q^{2}$ is equivalent to the optimization
problem for model A with a mobility $\mu=\mu^{\prime}q^{2}$. A method to
numerically compute the protocol minimizing $\mathcal{J}^{\lambda}_{q}$ exists
[52] for a harmonic trap. For any driving protocol, this method relies on:
$(1)$ To solve a set of ordinary differential equations (ODEs) for the
linearized (Gaussian) driving dynamics in eq. 79. $(2)$ the solution of ODEs
is used to compute $\langle W_{q}\rangle$ and $\sigma_{q}$, equivalently
leading to $\mathcal{J}^{\lambda}_{q}$. $(3)$ run a Stochastic gradient
descent (SGD) algorithm to minimize $\mathcal{J}^{\lambda}_{q}$ recursively.
We generalise this method for any value of $\mu$ [55]. We perform the SGD
algorithm with $3.10^{5}$ steps which numerically computes the optimal driving
protocol $K^{*}(q,t)$ that minimizes $\mathcal{J}^{\lambda}_{q}$. For $\mu=1$,
as a cross check, we recover results form Ref.[52].

We find a first-order phase transition for POP similar to liquid-gas phase
transition [53], where $\mu$ is a relevant control parameter. The existence of
the phase transition is attributed to the non-trivial $\mu$ dependence of the
$K^{*}(q,t)$ in eq. 7, thus $\mathcal{J}_{q}^{\lambda}$ is a nonlinear
function of $\mu$. The first order phase transition is characterized by a kink
in $\mathcal{J}^{\lambda}_{q}$, hence $\mathcal{J}^{\lambda}_{q}$ has
different slopes below and above a critical value of $\lambda$,
$\lambda_{c}=0.43$, fig. 2(b). Slow change of $\mu$ builds up the kink, fig.
2(b). Different slopes of $\mathcal{J}^{\lambda}_{q}$ correspond to two
different minima of $\mathcal{J}^{\lambda}_{q}$ and the universality classes
of the optimal driving protocols [53, 54]. The minima exchange global
stability at $\lambda=\lambda_{c}$. Optimal protocols for
$\lambda\geq\lambda_{c}$ belong to the universality class of mean-dominated
optimal protocols which are characterized by a discontinuity at the $t=0$ and
$t=\tau$ and a smooth interpolation for the driving in between. On the
contrary, $\lambda\leq\lambda_{c}$ belongs to the universality class of
fluctuations-dominated optimal protocols which are constant over the driving
time and exhibit discontinuity at the $t=0$ and $t=\tau$. Refer to fig 3b in
[55].

We observe a discontinuity in the PF for the critical value of
$\lambda=\lambda_{c}$, such that the non-convex part of the PF for
$\lambda=\lambda_{c}$ cannot be accessed with $\mathcal{J}^{\lambda}_{q}$
given by eq. 9, fig. 2(e). To quantify the first-order phase transition using
the hysteresis, we start with $\lambda=0,1$ and slowly change the control
parameter $\lambda$, in the forward (blue line) and the backward direction
(orange line) respectively. In fig. 2(c), the hysteresis loop for
$d\mathcal{J}^{\lambda}_{q}/d\lambda$ and $\mathcal{J}^{\lambda}_{q}$ (in
inset) is plotted. The optimal protocol of the previous value of $\lambda$ is
used to minimize $\mathcal{J}^{\lambda}_{q}$ after varying $\lambda$. $10^{5}$
hysteresis steps were used for the convergence of SGD. For different Fourier
modes, the emergence of the $\lambda_{c}$ is summarized by a phase diagram in
fig. 2(f-h).

For $\mu=1$, at higher values of $q$, the ratio
$K_{i}(q)/K_{f}(q)=(\alpha_{i}+q^{2})/(\alpha_{f}+q^{2})\rightarrow 1$
approaches the quasi-static slow driving limit, hence a unique protocol
minimizes mean and fluctuation of work so that it obeys the linear relation
($\langle W_{q}\rangle=\langle\Delta\mathcal{F}\rangle+\sigma_{q}^{2}/\mu T$)
[64, 65, 66, 67]. In this limit, PF is a single point in $\langle
W_{q}\rangle-\sigma_{q}$ space. For $\mu\geq 1$, we observe a gradual shift
from the slow driving limit imposed by $K_{i}(q)/K_{f}(q)\rightarrow 1$, such
that the mean-dominated and fluctuation-dominated protocols are different. At
a critical value $\mu_{c}$, the first order phase transition emerges, fig.
2(d). The phase diagram for different $q$ is summarized in fig 3c [55].

POP for a single Fourier mode does not encapsulate the effect of all Fourier
modes. $N$ independent Fourier modes are required for a discrete lattice with
$N$ sites such that $l^{th}$ Fourier mode reads $q=l\pi/N$. The total work
$W=\sum_{q}W_{q}$ is a sum of work needed to drive each Fourier mode $W_{q}$
111The linear (Gaussian) driving dynamics of each Fourier mode are decoupled
and hence evolve independently and justify this decomposition. We propose a
POP for $W$ composed of work needed to drive all Fourier modes, such that
$\langle W\rangle=\sum_{q}\langle W_{q}\rangle$ and $\sigma=\sqrt{\langle
W^{2}\rangle-\langle W\rangle^{2}}=\sqrt{\sum_{q}\sigma_{q}^{2}}$. The total
cost function $\mathcal{J}_{\lambda}=\lambda\langle
W\rangle+(1-\lambda)\sigma$ for all Fourier modes can not be expressed as a
sum of $\mathcal{J}^{\lambda}_{q}$. Hence $\mathcal{J}_{\lambda}$ does not
satisfy the physical constraint that the dynamics of each Fourier mode are
decoupled and independent. To address this physical constraint, we redefine
the cost function eq. 9 to

$\mathcal{J}^{\lambda}_{q}=\lambda\frac{\langle
W_{q}\rangle}{T}+(1-\lambda)\frac{\sigma_{q}^{2}}{T^{2}}$ (10)

In the eq. 10 the SD of the work is replaced with the variance of the work.
The temperature $T$ is used to re-scale $\mathcal{J}^{\lambda}_{q}$ to a
dimensionless quantity. For the POP of a single Fourier mode with eq. 10, we
observe a first-order phase transition. The qualitative behaviour is the same
as for POP with eq. 9. The plots fig. 2(b-h) are reproduced with modified
critical values of $\lambda_{c}$ and a modified phase diagram for
$\mu-\lambda_{c}$ which qualitatively preserves the phase transition for the
POP using eq. 10. The total cost function $\mathcal{J}^{\lambda}$ for the
optimization of mean and variance of all Fourier modes reads:

$\mathcal{J}^{\lambda}=\lambda\sum_{q}\frac{\langle
W_{q}\rangle}{T}+(1-\lambda)\sum_{q}\frac{\sigma_{q}^{2}}{T^{2}}=\sum_{q}\mathcal{J}_{q}^{\lambda}$
(11)

$\mathcal{J}^{\lambda}$ becomes the sum of $\mathcal{J}^{\lambda}_{q}$ for the
$q^{th}$ Fourier mode. Hence physical effects of N-independent Fourier modes
are summarized by analysing a harmonic trap for a single Fourier mode as
previously done. Further, we run a SGD algorithm to solve POP using the total
cost function, eq. 11. The driven dynamics of each Fourier mode evolve
independently as given by eq. 79, hence $\langle W_{q}\rangle$ and
$\sigma_{q}$ are computed independently for each Fourier mode as previously
done for the POP of a single Fourier mode. To optimize $\mathcal{J}^{\lambda}$
with the SGD, for each Monte-Carlo step, from the set of protocols {K(q,t)} we
randomly choose a protocol of $l^{th}$ Fourier mode and introduce a stochastic
perturbation. The new protocol is accepted if it minimizes the cost function
$\mathcal{J}^{\lambda}$. The POP with eq. 11 is summarized in fig. 2(i), the
PF computed corresponding to the POP of the eq. 11 is a linear combination of
the PF found for each independent Fourier mode that minimizes the eq. 10.

Conclusions and outlook We propose a generic framework to optimize the average
work finite-time driving of an order parameter field of an interacting multi-
body system modelled. Here, driving is defined as driving from one state to
another state. Using the methods from OTT and ST, we compute optimal driving
protocols that minimize mean irreversible work exhibiting discontinuities
[41]. We use entropy production as a finite time correction term for the
irreversible work which was shown to be proportional to the $\mathcal{W}_{2}$.
Further, we propose a multi-objective optimisation problem for the mean and
variance of the stochastic work. It reveals a first-order phase transition for
Pareto optimization and the mobility of the dynamics is a relevant control
parameter characterizing the phase transition. Moreover, a POP to optimize the
mean and variance of work for the field dynamics is reduced to the set of the
individual POP for each Fourier mode.

We propose four extensions of this work. $(i)$ To optimize the mean work for a
change of phase in finite time with multiple order parameters. In such cases,
the cross-correlation between the order parameters plays an important role
[69]. Without a cross-correlation between the order parameters, the dynamics
of the order parameter can be decoupled and the framework proposed in this
letter can be applied independently. $(ii)$ To relax the assumption of the
Gaussian small white noise used in this letter, hence applicable to a broader
class of biophysical systems that exhibit memory effects. $(iii)$ A recent
experimental development of the optimal control of the finite-time single
particle can be extended to many-body experimental systems 222The theoretical
analysis of various many-body systems done using the field theories has
experimental relevance. For instance, liquid crystals, liquid-liquid phase
separation, superconducting alloys and thin films etc., which could
experimentally verify our results [71, 72]. $(iv)$ Building upon [73, 74], one
can design more sophisticated irreversible thermodynamic cycles for the field-
theoretical description of the many-body systems 333Our framework can be seen
as one stroke of the traditional four-stroke engine. It is unclear how one can
optimize the efficiency over a cycle which is expected to have more
sophisticated optimal driving protocols.

## References

  * Crooks [2007] G. E. Crooks, Phys. Rev. Lett. 99, 100602 (2007).
  * Feng and Crooks [2008] E. H. Feng and G. E. Crooks, Phys. Rev. Lett. 101, 090602 (2008).
  * Machta [2015] B. B. Machta, Phys. Rev. Lett. 115, 260603 (2015).
  * Ito [2018] S. Ito, Phys. Rev. Lett. 121, 030605 (2018).
  * Brandner and Saito [2020] K. Brandner and K. Saito, Phys. Rev. Lett. 124, 040602 (2020).
  * Van Vu and Hasegawa [2021] T. Van Vu and Y. Hasegawa, Phys. Rev. Lett. 126, 010601 (2021).
  * Ito and Dechant [2020] S. Ito and A. Dechant, Phys. Rev. X 10, 021056 (2020).
  * Yoshimura and Ito [2021] K. Yoshimura and S. Ito, Phys. Rev. Research 3, 013175 (2021).
  * Nakazato and Ito [2021] M. Nakazato and S. Ito, Geometrical aspects of entropy production in stochastic thermodynamics based on wasserstein distance (2021), arXiv:2103.00503 [cond-mat.stat-mech] .
  * Sivak and Crooks [2012] D. A. Sivak and G. E. Crooks, Phys. Rev. Lett. 108, 190602 (2012).
  * Sivak and Crooks [2016] D. A. Sivak and G. E. Crooks, Phys. Rev. E 94, 052106 (2016).
  * Zulkowski and DeWeese [2014] P. R. Zulkowski and M. R. DeWeese, Phys. Rev. E 89, 052140 (2014).
  * Zulkowski and DeWeese [2015] P. R. Zulkowski and M. R. DeWeese, Phys. Rev. E 92, 032113 (2015).
  * Rotskoff  _et al._ [2017] G. M. Rotskoff, G. E. Crooks, and E. Vanden-Eijnden, Phys. Rev. E 95, 012148 (2017).
  * Zulkowski  _et al._ [2012] P. R. Zulkowski, D. A. Sivak, G. E. Crooks, and M. R. DeWeese, Phys. Rev. E 86, 041148 (2012).
  * Zulkowski  _et al._ [2013] P. R. Zulkowski, D. A. Sivak, and M. R. DeWeese, PLOS ONE 8, 1 (2013).
  * Bonança and Deffner [2014] M. V. S. Bonança and S. Deffner, The Journal of Chemical Physics 140, 244119 (2014), https://doi.org/10.1063/1.4885277 .
  * Deffner [2014] S. Deffner, Journal of Physics B: Atomic, Molecular and Optical Physics 47, 145502 (2014).
  * Large and Sivak [2019] S. J. Large and D. A. Sivak, Journal of Statistical Mechanics: Theory and Experiment 2019, 083212 (2019).
  * Deffner and Bonança [2020] S. Deffner and M. V. S. Bonança, EPL (Europhysics Letters) 131, 20001 (2020).
  * Hohenberg and Halperin [1977] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
  * Bray [2002] A. J. Bray, Advances in Physics 51, 481 (2002), https://doi.org/10.1080/00018730110117433 .
  * Chaikin and Lubensky [1995] P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995) p. 202.
  * Cates [2012] M. E. Cates, Complex fluids: The physics of emulsions (2012), arXiv:1209.2290 [cond-mat.soft] .
  * Kawasaki [1994] K. Kawasaki, Physica A: Statistical Mechanics and its Applications 208, 35 (1994).
  * Dean [1996] D. S. Dean, Journal of Physics A: Mathematical and General 29, L613 (1996).
  * Leonard  _et al._ [2013] T. Leonard, B. Lander, U. Seifert, and T. Speck, The Journal of Chemical Physics 139, 204109 (2013), https://doi.org/10.1063/1.4833136 .
  * Seifert [2005] U. Seifert, Phys. Rev. Lett. 95, 040602 (2005).
  * Seifert [2012] U. Seifert, Reports on Progress in Physics 75, 126001 (2012).
  * Van den Broeck and Esposito [2015] C. Van den Broeck and M. Esposito, Physica A: Statistical Mechanics and its Applications 418, 6 (2015), proceedings of the 13th International Summer School on Fundamental Problems in Statistical Physics.
  * Sekimoto [2010] K. Sekimoto,  _Stochastic Energetics_ (Courier Corporation, 2010).
  * Andresen [2011] B. Andresen, Angewandte Chemie International Edition 50, 2690 (2011).
  * Esposito and den Broeck [2011] M. Esposito and C. V. den Broeck, EPL (Europhysics Letters) 95, 40004 (2011).
  * Villani [2003] C. Villani,  _Topics in Optimal Transportation (American Mathematical Society)_ (American Mathematical Society, Providence, 2003).
  * Benamou and Brenier [2000] J.-D. Benamou and Y. Brenier, Numerische Mathematik 84, 375 (2000).
  * Santambrogio [2015] F. Santambrogio,  _Optimal Transport for Applied Mathematicians(Progress in Nonlinear Differential Equations and Their Applications Volume 87)_ (Springer, 2015).
  * Aurell  _et al._ [2011] E. Aurell, C. Mejía-Monasterio, and P. Muratore-Ginanneschi, Physical review letters 106, 250601 (2011).
  * Aurell  _et al._ [2012a] E. Aurell, C. Mejía-Monasterio, and P. Muratore-Ginanneschi, Phys. Rev. E 85, 020103 (2012a).
  * Aurell  _et al._ [2012b] E. Aurell, K. Gawedzki, C. Mejía-Monasterio, R. Mohayaee, and P. Muratore-Ginanneschi, Journal of Statistical Physics 147, 487 (2012b).
  * Dechant and Sakurai [2019] A. Dechant and Y. Sakurai, Thermodynamic interpretation of wasserstein distance (2019), arXiv:1912.08405 [cond-mat.stat-mech] .
  * Schmiedl and Seifert [2007] T. Schmiedl and U. Seifert, Phys. Rev. Lett. 98, 108301 (2007).
  * Gomez-Marin  _et al._ [2008] A. Gomez-Marin, T. Schmiedl, and U. Seifert, The Journal of Chemical Physics 129, 024114 (2008), https://doi.org/10.1063/1.2948948 .
  * Zhang [2020a] Y. Zhang, EPL (Europhysics Letters) 128, 30002 (2020a).
  * Zhang [2020b] Y. Zhang, Journal of Statistical Physics 178, 1336 (2020b).
  * Proesmans  _et al._ [2020a] K. Proesmans, J. Ehrich, and J. Bechhoefer, Phys. Rev. Lett. 125, 100602 (2020a).
  * Proesmans  _et al._ [2020b] K. Proesmans, J. Ehrich, and J. Bechhoefer, Phys. Rev. E 102, 032105 (2020b).
  * Proesmans and Bechhoefer [2021] K. Proesmans and J. Bechhoefer, EPL (Europhysics Letters) 133, 30002 (2021).
  * Movilla Miangolarra  _et al._ [2021] O. Movilla Miangolarra, R. Fu, A. Taghvaei, Y. Chen, and T. T. Georgiou, Phys. Rev. E 103, 062103 (2021).
  * Fu  _et al._ [2020] R. Fu, O. M. Miangolarra, A. Taghvaei, Y. Chen, and T. T. Georgiou, 2020 59th IEEE Conference on Decision and Control (CDC) , 3034 (2020).
  * Taghvaei  _et al._ [2022] A. Taghvaei, O. M. Miangolarra, R. Fu, Y. Chen, and T. T. Georgiou, IEEE Control Systems Letters 6, 434 (2022).
  * Chen  _et al._ [2020] Y. Chen, T. T. Georgiou, and A. Tannenbaum, IEEE Transactions on Automatic Control 65, 2979 (2020).
  * Solon and Horowitz [2018] A. P. Solon and J. M. Horowitz, Phys. Rev. Lett. 120, 180605 (2018).
  * Seoane and Solé [2015] L. F. Seoane and R. Solé, Phys. Rev. E 92, 032807 (2015).
  * Seoane and Solé [2015] L. F. Seoane and R. V. Solé, A multiobjective optimization approach to statistical mechanics (2015), arXiv:1310.6372 [cond-mat.stat-mech] .
  * si [2024] Supplemental information (2024).
  * Miettinen [2009] K. Miettinen,  _Nonlinear Multiobjective Optimization_ (Springer, 2009).
  * Shoval  _et al._ [2012] O. Shoval, H. Sheftel, G. Shinar, Y. Hart, O. Ramote, A. Mayo, E. Dekel, K. Kavanagh, and U. Alon, Science 336, 1157 (2012), https://www.science.org/doi/pdf/10.1126/science.1217405 .
  * Noor and Milo [2012] E. Noor and R. Milo, Science 336, 1114 (2012), https://www.science.org/doi/pdf/10.1126/science.1223193 .
  * Ronellenfitsch and Katifori [2019] H. Ronellenfitsch and E. Katifori, Phys. Rev. Lett. 123, 248101 (2019).
  * Ronellenfitsch and Katifori [2016] H. Ronellenfitsch and E. Katifori, Phys. Rev. Lett. 117, 138301 (2016).
  * Guillier  _et al._ [2017] S. Guillier, V. Muñoz, J. Rogan, R. Zarama, and J. Valdivia, Physica A: Statistical Mechanics and its Applications 467, 465 (2017).
  * Santoro  _et al._ [2018] A. Santoro, V. Latora, G. Nicosia, and V. Nicosia, Phys. Rev. Lett. 121, 128302 (2018).
  * Dinis  _et al._ [2020] L. Dinis, J. Unterberger, and D. Lacoste, Europhysics Letters 131, 60005 (2020).
  * Mazonka and Jarzynski [1999] O. Mazonka and C. Jarzynski, Exactly solvable model illustrating far-from-equilibrium predictions (1999), arXiv:9912121 [nlin.PS] .
  * Speck and Seifert [2004] T. Speck and U. Seifert, Phys. Rev. E 70, 066112 (2004).
  * Speck [2011] T. Speck, Journal of Physics A: Mathematical and Theoretical 44, 305001 (2011).
  * Miller and Mehboudi [2020] H. J. D. Miller and M. Mehboudi, Phys. Rev. Lett. 125, 260602 (2020).
  * Note [1] The linear (Gaussian) driving dynamics of each Fourier mode are decoupled and hence evolve independently and justify this decomposition.
  * Dechant and Sasa [2021] A. Dechant and S.-i. Sasa, Phys. Rev. X 11, 041061 (2021).
  * Note [2] The theoretical analysis of various many-body systems done using the field theories has experimental relevance. For instance, liquid crystals, liquid-liquid phase separation, superconducting alloys and thin films etc.
  * Dago and Bellon [2022] S. Dago and L. Bellon, Phys. Rev. Lett. 128, 070604 (2022).
  * Dago  _et al._ [2022] S. Dago, J. Pereda, S. Ciliberto, and L. Bellon, Journal of Statistical Mechanics: Theory and Experiment 2022, 053209 (2022).
  * Frim and DeWeese [2022a] A. G. Frim and M. R. DeWeese, Phys. Rev. E 105, L052103 (2022a).
  * Frim and DeWeese [2022b] A. G. Frim and M. R. DeWeese, Phys. Rev. Lett. 128, 230601 (2022b).
  * Note [3] Our framework can be seen as one stroke of the traditional four-stroke engine. It is unclear how one can optimize the efficiency over a cycle which is expected to have more sophisticated optimal driving protocols.
  * Gardiner [2009] C. W. Gardiner,  _Stochastic Methods: A Handbook for the Natural and Social Sciences_ (Springer, 2009).

Supplemental Material
Atul Tanaji Mohite

_Department of Theoretical Physics and Center for Biophysics, Saarland
University, Saarbrücken, Germany_

###### Contents

  1. I Initial and Final Equilibrium probability distribution
    1. I.1 GL functional for Model A dynamics
    2. I.2 GL functional for Model B dynamics
  2. II Optimal transport theory and stochastic thermodynamics
    1. II.1 Optimal transport theory for fields
    2. II.2 Optimal driving protocol
    3. II.3 Minimum heat corresponding to the optimal protocol
    4. II.4 Minimum work corresponding to the optimal protocol
  3. III Analytical expression of wasserstein distance for model A dynamics
    1. III.1 Infinite frequency cutoff
    2. III.2 Finite frequency cutoff
      1. III.2.1 One dimensional space
      2. III.2.2 Two dimensional space
      3. III.2.3 Three dimensional space
  4. IV Analytical expression of Wasserstein distance for model B dynamics
    1. IV.1 Infinite frequency cutoff
    2. IV.2 Finite frequency cutoff
      1. IV.2.1 Two dimensional space
  5. V Boundary term for the heat
    1. V.1 Two dimensional space
  6. VI Pareto optimization with mean and standard deviation of the work
    1. VI.1 Single Fourier mode case
    2. VI.2 Multiple Fourier mode case
  7. VII Numerical scheme and parameters
  8. VIII Other plots
  9. IX Parameters used for plots in the maintext

##  I Initial and Final Equilibrium probability distribution

The equilibrium probability distribution for order parameter $\phi$ can be
computed using standard methods. We restrict ourselves to a small noise limit
for the fluctuations of the order parameter $\phi$. This ensures that the
equilibrium probability distribution $\mathcal{P}([\phi_{q}])$ is Gaussian.
Thus computing the mean and the variance of the probability distribution
ensures a complete probability distribution.

###  I.1 GL functional for Model A dynamics

In this section, we derive correlation functions for Model A dynamics. We use
a GL functional which encapsulates a discontinuous first-order phase
transition. The GL functional is given by:

$\mathcal{H}[\phi]=\int\bigg{[}\frac{a}{2}\phi^{2}-\frac{c}{3}\phi^{3}+\frac{b}{4}\phi^{4}+\frac{\kappa}{2}(\nabla\phi)^{2}\bigg{]}dx$
(12)

Using small noise approximation for the field
$\phi(x)=\bar{\phi_{0}}+\sqrt{T}\delta\phi(x)$, we get the following equation
for the GL functional:

$\mathcal{H}[\phi]=\int\bigg{[}\frac{a}{2}\bar{\phi_{0}}^{2}-\frac{c}{3}\bar{\phi_{0}}^{3}+\frac{b}{4}\bar{\phi_{0}}^{4}\bigg{]}dx+\frac{T}{2}\int\frac{dq}{2\pi}\big{(}a-2c\bar{\phi_{0}}+3b\bar{\phi_{0}}^{2}+\kappa
q^{2}\big{)}\delta\phi(q)^{2}$ (13)

---

---

Figure 3: (a) Free energy functional as a function of density, for
$a>\frac{2c^{2}}{9b}$, $\phi=0$ is the stable global minimum of the free
energy functional. For $a<\frac{2c^{2}}{9b}$, another minimum of the free
energy functional with $\phi=\phi_{2}=\frac{c+\sqrt{c^{2}-4ab}}{2b}$ is the
global minimum. (b) Phase diagram for the first-order phase transition,
equation of the parabola is given by $9ab=2c^{2}$

Using the expression above, and the equipartition theorem, the correlation
function can be written as follows:

$\mathit{C}(q)=\langle\delta\phi(q)^{2}\rangle=\frac{1}{a-2c\bar{\phi_{0}}+3b\bar{\phi_{0}}^{2}+\kappa
q^{2}}=\frac{1}{\alpha_{0}+\kappa q^{2}}$ (14)

Where $\langle.\rangle$ denotes the sum over all noise realisations. Hence

$\mathit{C}_{i}(q)=\frac{1}{\alpha_{i}+\kappa q^{2}};\hskip
56.9055pt\mathit{C}_{f}(q)=\frac{1}{\alpha_{f}+\kappa q^{2}}$ (15)

From [23] the initial and final values of the order parameter corresponding to
the stable minimum are given by:

$\phi_{i}=0;\hskip
56.9055pt\phi_{f}=\frac{c_{f}+\sqrt{c_{f}^{2}-4a_{f}b}}{2b}$ (16)

Using eqs. 14, 15 and 16, $\alpha_{i}$ and $\alpha_{f}$ corresponding to the
initial and final control parameters are given by the following equation:

$\alpha_{i}=a_{i}$ (17a)
$\begin{split}\alpha_{f}&=a_{f}-2c_{f}\phi_{2}+3b\phi_{2}^{2}=2b\phi_{2}^{2}-c_{f}\phi_{2}=c_{f}\phi_{2}-2a_{f}\\\
&=\frac{c_{f}^{2}-4a_{f}b+c_{f}\sqrt{c_{f}^{2}-4a_{f}b}}{2b}\end{split}$ (17b)

Hence the initial and final Gaussian equilibrium distributions are given by
the initial and final mean $\phi_{i}$ and $\phi_{f}$, eq. 16 and the initial
and the final correlation function $C_{i}(q)$ and $C_{f}(q)$, eqs. 15 and 17.

###  I.2 GL functional for Model B dynamics

In this section, we compute correlation functions for Model B dynamics. We use
a GL functional that leads to a phase transition for a conserved order
parameter $\phi$. Thus $c=0$ is imposed, to recover the standard GL functional
to recover the liquid-liquid phase separation:

$\mathcal{H}[\phi]=\int\bigg{[}\frac{a}{2}\phi^{2}+\frac{b}{4}\phi^{4}+\frac{\kappa}{2}(\nabla\phi)^{2}\bigg{]}dx$
(18)

---

---

Figure 4: (a) Free energy functional as a function of density, for
$a>-b\bar{\phi}^{2}$, $\phi=0$ is the stable global minimum of the free energy
functional. For $a<-b\bar{\phi}^{2}$, two stable minima with a finite value of
the order parameter $\phi$ emerges. (b) Phase diagram for the liquid-liquid
phase separation, equation of the parabola is given by $a=-b\bar{\phi}^{2}$

Using small noise approximation for the field
$\phi(x)=\bar{\phi_{0}}+\sqrt{T}\delta\phi(x)$, we get the following equation
for the free energy functional:

$\mathcal{H}[\phi]=\int\bigg{[}\frac{a}{2}\bar{\phi_{0}}^{2}+\frac{b}{4}\bar{\phi_{0}}^{4}\bigg{]}dx+\frac{T}{2}\int\frac{dq}{2\pi}\big{(}a+3b\bar{\phi_{0}}^{2}+\kappa
q^{2}\big{)}\delta\phi(q)^{2}$ (19)

Using the expression above the correlation function can be written as follows:

$\mathit{C}(q)=\langle\delta\phi(q)^{2}\rangle=\frac{1}{a+3b\bar{\phi_{0}}^{2}+\kappa
q^{2}}=\frac{1}{\alpha_{0}+\kappa q^{2}}$ (20)

Where $\langle.\rangle$ denotes the sum over all noise realisations. Hence

$\mathit{C}_{i}(q)=\frac{1}{\alpha_{i}+\kappa q^{2}};\hskip
56.9055pt\mathit{C}_{f}(q)=\frac{1}{\alpha_{f}+\kappa q^{2}}$ (21)

eqs. 21 and 22. For liquid-liquid phase separation, order parameter $\phi$ is
conserved hence $\phi_{1}=\phi_{2}=\bar{\phi}_{0}$. Using eqs. 20 and 19,
$\alpha_{i}$ and $\alpha_{f}$ are given by the following equation:

$\alpha_{i}=a_{i}+3b\bar{\phi_{0}}^{2};\hskip
56.9055pt\alpha_{f}=a_{f}+3b\bar{\phi_{0}}^{2}$ (22)

Hence the initial and final Gaussian equilibrium distributions are given by
the initial and final mean $\phi_{i}=\phi_{f}=\bar{\phi}_{0}$ and the initial
and the final correlation function $C_{i}(q)$ and $C_{f}(q)$,

##  II Optimal transport theory and stochastic thermodynamics

###  II.1 Optimal transport theory for fields

In section section I, we computed the initial and final equilibrium Gaussian
probability distributions. The current section is dedicated to computing the
optimal driving protocol that minimizes the entropy production or irreversible
work to drive the probability distribution from the initial phase to the final
phase in a finite driving time $\tau$. Further, we derive the optimal driving
protocol to compute the minimum energy cost required to drive the system from
the initial phase to the final phase.

If the field for the order parameter $\phi$ follows non-conserved or conserved
dynamics, then the Langevin equation for model A or model B with an additive
noise can be written as follows:

$\dot{\phi}(x,t)=-\mu(i\nabla)^{2\beta}\bigg{(}\frac{\delta\mathcal{H}}{\delta\phi}\bigg{)}+\Lambda(x,t)$
(23a) $\langle\Lambda(x,t)\Lambda(x^{\prime},t^{\prime})\rangle=2\mu
T(i\nabla)^{2\beta}\delta(x-x^{\prime})\delta(t-t^{\prime})$ (23b)

In eq. 23, $\beta=\\{0,1\\}$ corresponds to the model A or moedel B dynamics
respectively. $\mu$ is mobility and $T$ is the temperature of heat bath. By
plugging in $\beta=1$, the following Langevin equation for model B dynamics is
recovered.

$\dot{\phi}(x,t)=\nabla.\bigg{(}\mu\nabla.\frac{\delta\mathcal{H}}{\delta\phi}+\sqrt{2T\mu}\>\epsilon(x,t)\bigg{)}$
(24)

Similarly, for $\beta=0$ the Langevin equation for model A with an additive
noise is recovered:

$\dot{\phi}(x,t)=-\mu\frac{\delta\mathcal{H}}{\delta\phi}+\sqrt{2T\mu}\>\epsilon(x,t)$
(25)

Where $\epsilon$ is a standard Gaussian white noise satisfying
$\langle\epsilon(x,t)\rangle=0$ and
$\langle\epsilon(x,t)\epsilon(x^{\prime},t^{\prime})\rangle=\delta(x-x^{\prime})\delta(t-t^{\prime})$.
eq. 24 has gradient terms, which are not trivial to handle in real space. This
issue can be solved by introducing Fourier modes, leading to a simplification
of the gradient terms. We define a driving force in the Fourier space such
that for the equilibrium case $f_{q}=-\delta\mathcal{H}/\delta\phi_{q}$, in
general for driven dynamics $f_{q}=-\delta\mathcal{H}/\delta\phi_{q}$ does not
hold as $f_{q}$ is an arbitrary driving force employed to go from the initial
to the final state. In such a case, a generic expression for the driving
protocol is considered. Further, we take a Fourier transform of eqs. 23, 24
and 25. In Fourier space, eq. 23 can be written as follows:

$\dot{\phi_{q}}(t)=\mu(q,\phi)f_{q}+\Lambda_{q}$ (26a)
$\langle\Lambda_{q}(t)\Lambda_{q^{\prime}}(t^{\prime})\rangle=2\mu
T(-qq^{\prime})^{\beta}\delta(q+q^{\prime})\delta(t-t^{\prime})$ (26b)

Where mobility $\mu(q,\phi)=q^{2}\mu$ for the model B and $\mu(q,\phi)=\mu$
for model A. Using $\beta=1$ and $\beta=0$, we can recover the eq. 27 for
model B and eq. 28 for model A respectively:

$\dot{\phi}_{q}=\mu q^{2}f_{q}+iq\sqrt{2T\mu}\>\Lambda_{q}$ (27)
$\dot{\phi}_{q}=\mu f_{q}+\sqrt{2T\mu}\>\Lambda_{q}$ (28)

The Langevin equation eq. 26a can be equivalently converted into a Fokker-
Planck-Equation for probability distribution $P(\phi_{q},t)$. The FPE can be
written as follows[76]:

$\begin{split}\partial_{t}P&=\sum_{q}\frac{\partial}{\partial\phi_{q}}\bigg{[}-\mu(q,\phi)f_{q}P\bigg{]}+\sum_{q,q^{\prime}}\langle\Lambda(q,t)\Lambda(q^{\prime},t^{\prime})\rangle\bigg{[}\frac{\partial^{2}P}{\partial\phi_{q}\partial\phi
q^{\prime}}\bigg{]}\\\
&=\sum_{q}\frac{\partial}{\partial\phi_{q}}\bigg{[}-\mu(q,\phi)f_{q}P\bigg{]}+\sum_{q,q^{\prime}}T\mu_{(}q,\phi)\delta_{q+q^{\prime},0}\bigg{[}\frac{\partial^{2}P}{\partial\phi_{q}\partial\phi_{q^{\prime}}}\bigg{]}\\\
&=\sum_{q}\frac{\partial}{\partial\phi_{q}}\bigg{[}-\mu(q,\phi)f_{q}P\bigg{]}+\sum_{q}T\mu_{(}q,\phi)\bigg{[}\frac{\partial^{2}P}{\partial\phi_{q}\partial\phi_{-q}}\bigg{]}\end{split}$
(29)

In eq. 29, we take out the common factor of $\mu(q,\phi)$. We identify the
velocity $\mathcal{V}_{q}$ such that $\mathcal{V}_{q}$ is given by the
following formula:

$\begin{split}\mathcal{V}_{q}&=\bigg{[}\mu(q,\phi)f_{q}-\frac{T\mu(q,\phi)}{P}\bigg{(}\frac{\partial}{\partial\phi_{-q}}\bigg{)}P\bigg{]}=\mu(q,\phi)\bigg{[}f_{q}-\frac{T}{P}\bigg{(}\frac{\partial
P}{\partial\phi_{-q}}\bigg{)}\bigg{]}=\mu(q,\phi)\bigg{[}f_{q}-T\bigg{(}\frac{\partial\ln{P}}{\partial\phi_{-q}}\bigg{)}\bigg{]}\end{split}$
(30)

We have rearranged the factor of $\mu(q,\phi)$ in the main expression for
$\mathcal{V}_{q}$ for the sake of convenience. The Fokker-Planck equation for
the Langevin dynamics can be written as follows:

$\begin{split}\partial_{t}P&=-\sum_{q}\frac{\partial}{\partial\phi_{q}}\bigg{[}\mathcal{V}_{q}P\bigg{]}\end{split}$
(31)

The formula for the entropy production(irreversible work) is given by the
following equation [29]:

$\begin{split}\Sigma=\sum_{q}\int_{0}^{\tau}\int_{\phi_{q}}\frac{d\tau\>d\phi_{q}}{\mu(q,\phi,t)\>T}\mathcal{V}_{q}^{2}\>P\end{split}$
(32)

The total entropy production $\Sigma$ corresponds to the irreversible work for
the finite time driving from the initial state to the final state. Hence we
want to minimize the total entropy production $\Sigma$ given by equation eq.
32 under the constraint given by the Fokker-Planck eq. 30 and eq. 31. We
consider a minimization problem for entropy production with a Lagrange
multiplier for the constraint of the Fokker-Planck equation, in other words,
we compute:

$\begin{split}\Sigma&=\inf_{\mathcal{V},f,P,\rho,\psi}\int_{0}^{\tau}d\tau\int
d\phi_{q}\Bigg{[}\sum_{q}\frac{\mathcal{V}_{q}^{2}}{\mu(q,\phi,t)\>T}P-2\rho\Bigg{(}\partial_{t}P+\sum_{q}\frac{\partial}{\partial\phi_{q}}\Big{(}\mathcal{V}_{q}P\Big{)}\Bigg{)}\\\
&\hskip
14.22636pt-\sum_{q}2\psi_{q}\Bigg{(}\mathcal{V}_{q}-\mu(q,\phi,t)f_{q}+\mu(q,\phi,t)\bigg{(}T\>\frac{\partial\ln{P}}{\partial\phi_{-q}}\bigg{)}\Bigg{)}\Bigg{]}\end{split}$
(33)

The corresponding Euler-Lagrange equations read:

$\mathcal{V}_{q}:\frac{2\mathcal{V}_{q}P}{\mu(q,\phi,t)T}+2P\frac{\partial\rho}{\partial\phi_{q}}+2\psi_{q}=0$
(34a) $f_{q}:-2\mu(q,\phi,t)\psi_{q}=0$ (34b) $P:\hskip
5.69046pt\sum_{q}\frac{\mathcal{V}_{q}^{2}}{\mu(q,\phi,t)T}+2\partial_{t}\rho+\sum_{q}2\mathcal{V}_{q}\frac{\partial\rho}{\partial\phi_{q}}+\sum_{q}\frac{T\mu(q,\phi,t)}{P^{2}}\psi_{q}\frac{\partial
P}{\partial\phi_{-q}}-T\mu(q,\phi,t)\sum_{q}\frac{\partial}{\partial\phi_{-q}}\bigg{(}\frac{\psi_{q}}{P}\bigg{)}$
(34c)

Using eq. 34b and substituting it in eq. 34a and eq. 34c, we can simplify the
equations into the following form:

$\mathcal{V}_{q}=-\mu(q,\phi,t)T\frac{\partial\rho}{\partial\phi_{q}}$ (35a)
$\partial_{t}\rho=-\sum_{q}\frac{1}{2}\mu(q,\phi,t)T\bigg{(}\frac{\partial\rho}{\partial\phi_{q}}\bigg{)}^{2}=-\frac{\mathcal{V}^{2}}{2\mu(q,\phi,t)T}$
(35b)
$\partial_{t}P=\mu(q,\phi,t)T\sum_{q}\frac{\partial}{\partial\phi_{q}}\bigg{(}P\frac{\partial\rho}{\partial\phi_{q}}\bigg{)}$
(35c)

eq. 35 is analogous to [40]. The gist of the argument is to map eq. 35 onto
Monge-Kantorovich problem [35]. Monge-Kantorovich optimal transport problem
follows the minimization of the total entropy production $\Sigma^{*}$ for the
constraints given by the following equations.

$\Sigma^{*}=\int_{0}^{1}\int_{\chi_{q}}ds\>d\chi_{q}\parallel\nabla_{\chi_{q}}\tilde{\rho}(\chi_{q},s)\parallel^{2}\>P$
(36a)
$\partial_{s}\tilde{\rho}=-\frac{1}{2}\parallel\nabla_{\chi_{q}}\tilde{\rho}(\chi_{q},s)\parallel^{2}$
(36b)
$\partial_{s}\tilde{P}(\chi_{q},s)=\nabla_{\chi_{q}}\Big{(}\nabla_{\chi_{q}}\tilde{\rho}(\chi_{q},s)\tilde{P}(\chi_{q},s)\Big{)}$
(36c)

To map the set of eq. 35 onto the set of eq. 36, we scale the time and space.
The scaling is given by the following set of equations:

$s(t)=\frac{\int_{0}^{t}\mu(q,\phi,t)T(t)\>dt}{\int_{0}^{\tau}\mu(q,\phi,t)T(t)\>dt}=\frac{t}{\tau}\hskip
42.67912pt\chi_{q}=\frac{\phi_{q}}{\sqrt{\int_{0}^{\tau}\mu(q,\phi,t)T(t)\>dt}}$
(37)

For model A $\mu(q,\phi,t)=\mu$, on the other hand for model B,
$\mu(q,\phi,t)=\mu q^{2}$. Hence the relation between the entropy production
and Wasserstein distance can be used exactly for model A. For model B, the
mapping to Monge-Kantorovich [40, 35] includes not only the scaling of time
but also $\phi$ such that in q-space,

$\chi_{q}=\frac{\phi_{q}}{\sqrt{\int_{0}^{\tau}\mu(q,t)T}}=\frac{\phi_{q}}{\sqrt{\tau\mu
q^{2}T}}$ (38)

Thus we will use the mean and variance of the probability distribution for
$\chi$ for model B. The normalization due to mobility is independent of the
$q$ for the model A. Hence we can utilize the relation for the lower bound on
the entropy production with slight modifications given by the q-dependent
mobility for the model B. Another consequence of this is that for model B, the
term corresponding to the mean term does not contribute as
$\bar{\phi}_{i}=\bar{\phi}_{f}$. The lower bound on entropy production is
given by an Information theoretic measure called ’Wasserstein distance’
$\mathcal{W}_{2}(\mathcal{P}_{i}[\phi_{q}],\mathcal{P}_{f}[\phi_{q}])^{2}$

$\Sigma\geq\mathcal{W}_{2}\big{(}\mathcal{P}_{i}[\chi_{q}],\mathcal{P}_{f}[\chi_{q}]\big{)}^{2}$
(39)

Where the lower bound on the
$\mathcal{W}_{2}(\mathcal{P}_{i}[\phi_{q}],\mathcal{P}_{f}[\phi_{q}])^{2}$ is
given by the following expression:

$\begin{split}\mathcal{W}_{2}\big{(}\mathcal{P}_{i}[\chi_{q}],\mathcal{P}_{f}[\chi_{q}]\big{)}^{2}&\geq\Bigg{\\{}\Big{(}E_{f}[\chi_{q}]-E_{i}[\chi_{q}]\Big{)}^{2}+\bigg{(}\sqrt{E_{f}\Big{[}\big{(}\chi_{q}-E_{f}[\chi_{q}]\big{)}^{2}\Big{]}}-\sqrt{E_{i}\Big{[}\big{(}\chi_{q}-E_{i}[\chi_{q}]\big{)}^{2}\Big{]}}\bigg{)}^{2}\Bigg{\\}}\end{split}$
(40)

Where $E_{f}$ and $E_{i}$ in eq. 40 denotes expectation value over the final
and initial probability distributions respectively. The first term in eq. 40
is the contribution corresponding to the change of mean and the second term is
the contribution corresponding to the change of the co-variance of the
probability distribution. The inequality in eq. 40 is saturated when both
$\mathcal{P}_{i}[\chi_{q}]$ and $\mathcal{P}_{i}[\chi_{q}]$ are Gaussian,
which corresponds to the small noise limit case utilised to compute the
Gaussian equilibrium probability distributions in section I. Using eqs. 39 and
38, the entropy production for model A is given by

$\Sigma\geq\frac{\mathcal{W}_{2}\big{(}\mathcal{P}_{i}[\phi_{q}],\mathcal{P}_{f}[\phi_{q}]\big{)}^{2}}{\tau\>\mu\>T}$
(41)

Similarly, using eqs. 39 and 38, the entropy production for model B is given
by

$\Sigma\geq\frac{\mathcal{W}_{2}\big{(}\mathcal{P}_{i}[\phi_{q}],\mathcal{P}_{f}[\phi_{q}]\big{)}^{2}}{\tau\>\mu
q^{2}\>T};\hskip 56.9055pt\chi_{q}=\frac{\phi_{q}}{q}$ (42)

We have obtained the expression for the Wasserstein distance for model A and
model B dynamics respectively which gives a lower bound on entropy production
or irreversible work. We further generalize the results from [40, 35]. Using
the optimal transport map, we get expressions for the mean ($\phi^{*}(t)$) and
co-variance ($C^{*}(q,t)$) of the Fourier mode $\phi_{q}$ at any intermediate
time $t$. The exact expression is given by the following optimal transport map
equations:

$\phi^{*}(q,t)=\frac{t}{\tau}\bar{\phi_{f}}\>\delta_{q,0}+\bigg{(}1-\frac{t}{\tau}\bigg{)}\bar{\phi_{i}}\>\delta_{q,0}$
(43a)
$C^{*}(q,t)=\Bigg{(}\frac{t}{\tau}\sqrt{C_{f}(q)}+\bigg{(}1-\frac{t}{\tau}\bigg{)}\sqrt{C_{i}(q)}\Bigg{)}^{2}$
(43b)

The linearity of the optimal transport map ensures that if the initial and
final probability distributions are Gaussian then at any intermediate time
$t$, the probability distribution is Gaussian[35].

###  II.2 Optimal driving protocol

We assume linear ansatz for the driving protocol
$f_{q}(t)=K(q,t)\phi_{q}+a(q,t)$ as optimal transport map ensures gaussianity
of the probability distribution at intermediate time $t$, given the initial
and final probability distributions are Gaussian. Thus the validity of the
linear driving ansatz hold true. The Langevin equation for linear driving of
each independent Fourier mode is written as follows:

$\partial_{t}\phi_{q}=\mu(q)\Big{(}K(q,t)\mathbf{\phi}_{q}+a(q,t)\Big{)}+\Lambda_{q}$
(44)

The Fokker-Plank equation for the probability distribution $P(\phi_{q},t)$ at
an intermediate driving time $t$ can be written as follows:

$\partial_{t}P(\mathbf{\phi}_{q},t)=-\mu(q)\frac{\partial}{\partial\phi_{q}}\bigg{[}K(q,t)\mathbf{\phi}_{q}+a(t)-T\frac{\partial}{\partial\phi_{-q}}\bigg{]}P(\mathbf{\phi}_{q},t)$
(45)

The optimal force $f_{q}(t)$ satisfies the relation
$f_{q}(t)=K(q,t)\phi_{q}+a(q,t)$. Using Ito’s lemma, The equations for the
time evolution of mean and variance reads:

$d_{t}\phi_{q}^{*}(t)=\mu(q)\Big{(}K(q,t)\phi_{q}+a(q,t)\Big{)}$ (46a)
$d_{t}C^{*}(q,t)=\mu(q)\Big{(}2K(q,t)C^{*}(q,t)+2T\Big{)}$ (46b)

Set of eq. 46 should be consistent with the optimal transport map given by eq.
43 at any intermediate driving time $t$. Using eq. 43 and plugging it into eq.
46 leads to a closed-form expression for the optimal driving protocol. The
final expression for the optimal driving force
$f_{q}^{*}(t)=K(q,t)\phi_{q}+a(t)$ is given as follows:

$f_{q}^{*}(t)=\frac{d_{t}\phi^{*}_{q}(t)}{\mu(q)}-\frac{1}{C^{*}(q,t)}\bigg{(}T-\frac{1}{2\mu(q)}d_{t}C^{*}(q,t)\bigg{)}(\phi_{q}-\phi^{*}_{q}(t))$
(47)

For $q\neq 0$, from eq. 43a $\phi_{q}^{*}(t)=0$ which implies $a^{*}(q,t)=0$,
hence for $q\neq 0$ the expression for the optimal force can be simplified to
the following form:

$\begin{split}f_{q}^{*}(t)&=K^{*}(q,t)\phi_{q}=-\frac{1}{C^{*}(q,t)}\bigg{(}T-\frac{1}{2\mu(q)}d_{t}C^{*}(q,t)\bigg{)}\phi_{q}\\\
&=\Bigg{[}\bigg{(}\frac{\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}}{\mu(q)\>\tau\>\sqrt{C^{*}(q,t)}}\bigg{)}-\frac{T}{C^{*}(q,t)}\Bigg{]}\phi_{q}\end{split}$
(48)

This implies that the equation for the dynamics of the Fourier mode is given
by the following expression:

$\partial_{t}\phi_{q}(t)=\mu(q)f_{q}(t)+\eta_{q}(t)=\mu(q)K^{*}(q,t)\phi_{q}(t)+\eta_{q}(t)$
(49)

Where $\eta_{q}(t)$ is a noise term for the Langevin equation, the exact
expression for a non-conserved noise can be written by
$\eta_{q}(t)=\sqrt{2T\mu}\>\Lambda_{q}(t)$. Similarly for a conserved noise
$\eta_{q}(t)=iq\>\sqrt{2T\mu}\>\Lambda_{q}(t)$.

For the case of $q=0$,
$\phi_{0}^{*}(t)=(t\bar{\phi}_{f}/\tau+(1-t/\tau)\bar{\phi}_{i})$, hence
$\partial_{t}\phi^{*}(t)=(\bar{\phi}_{f}-\bar{\phi}_{i})/\tau$, hence
$a(0,t)=(\bar{\phi}_{f}-\bar{\phi}_{i})/(\tau\mu(0))$. Rearranging eq. 47, we
can deduce that for $q=0$, $\delta\phi_{0}(t)=\phi_{0}(t)-\phi_{0}^{*}(t)$
follows the same driven dynamics as $\phi_{q}(t)$ with a trap stiffness
$K^{*}(q,t)$, the explicit equation is written as follows:

$\partial_{t}\delta\phi_{0}(t)=\mu(0)K^{*}(0,t)\delta\phi_{0}(t)+\eta_{0}(t)$
(50a)
$\begin{split}\partial_{t}\phi_{0}(t)=\mu(0)f_{0}(t)+\eta_{0}(t)&=\partial_{t}\phi_{0}^{*}(t)+\mu(0)K^{*}(0,t)(\phi_{0}(t)-\phi_{0}^{*}(t))+\eta_{q}(t)\\\
&=\mu(0)a^{*}(0,t)+\mu(0)K^{*}(0,t)\delta\phi_{0}(t)+\eta_{0}(t)\\\
&=\mu(0)a^{*}(0,t)+\partial_{t}\delta\phi_{0}(t)\end{split}$ (50b)

Where $K^{*}(q,t)$ is given by the eq. 48 and
$a^{*}(0,t)=\partial_{t}\phi_{0}^{*}(t)/\mu(0,t)=(\bar{\phi}_{f}-\bar{\phi}_{i})/(\tau\mu(0))$.
In other words, the dynamics of the zeroth Fourier mode are offset-ed by mean
term $\phi_{0}^{*}(t)$. eq. 50a and eq. 50b are equivalent ways of writing the
dynamics of the zeroth Fourier mode. Using eq. 50a is easier for the sake of
calculations and numerical simulations.

###  II.3 Minimum heat corresponding to the optimal protocol

Having formulated the optimal driving protocol for the dynamics of the Fourier
modes, we can use the standard framework of Stochastic thermodynamics to
compute the heat production for the system. The probability distribution for
the realisation of a trajectory can be written as follows:

$\mathbb{P}\thicksim e^{-\mathcal{A}}$ (51)

Where $\mathcal{A}$ in the eq. 51 is the action defined as follows:

$\mathcal{A}=\sum_{q}\int_{0}^{\tau}dt\frac{1}{4T\mu(q)}\Big{[}\dot{\phi}_{q}(t)-\mu(q)f^{*}_{q}(t)\Big{]}^{2}$
(52)

In stochastic thermodynamics, the heat for a single trajectory is given by
$\mathcal{Q}=T\ln{(\mathbb{P}_{+}/\mathbb{P}_{-})}$ [29]. Where
$\mathbb{P}_{+}$ and $\mathbb{P}_{-}$ are the probability distribution for the
forward and time-reversed process. The average heat over all trajectory
realisations is given by:

$\begin{split}\langle\mathcal{Q}\rangle&=T\left\langle\ln{\frac{\mathbb{P_{+}}}{\mathbb{P_{-}}}}\right\rangle=-T\left\langle
A_{+}-A_{-}\right\rangle=\sum_{q}\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{q}(t)\>f^{*}_{q}(t)\right\rangle\\\
&=\sum_{q\neq
0}\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{q}(t)\>f^{*}_{q}(t)\>\right\rangle+\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{0}(t)\>f^{*}_{0}(t)\>\right\rangle\\\
&=\sum_{q\neq
0}\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{q}(t)\>\phi_{q}(t)\>\right\rangle
K^{*}(q,t)+\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{0}(t)\>f^{*}_{0}(t)\>\right\rangle\\\
\end{split}$ (53)

In eq. 53$,\langle.\rangle$ signifies the sum over all noise realisation of
the trajectory. Using eq. 50b, the expression for the heat can be resolved to
the following form:

$\begin{split}\left\langle\mathcal{Q}\right\rangle&=\sum_{q\neq
0}\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{q}(t)\>\phi_{q}(t)\>\right\rangle
K^{*}(q,t)+\int_{0}^{\tau}dt\>\Big{\langle}\big{(}\mu(0)a^{*}(0,t)+\partial_{t}\delta\phi_{0}(t)\big{)}\>\big{(}a^{*}(0,t)+K^{*}(0,t)\delta\phi_{0}(t)\big{)}\>\Big{\rangle}\\\
&=\sum_{q\neq
0}\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{q}(t)\>\phi_{q}(t)\>\right\rangle
K^{*}(q,t)\\\ &\hskip
28.45274pt+\int_{0}^{\tau}dt\>\bigg{(}\mu(0)a^{*}(t)^{2}+a^{*}(t)\langle\partial_{t}\delta\phi_{0}(t)\rangle+\mu(0)a^{*}(0,t)K^{*}(0,t)\left\langle\delta\phi_{0}(t)\right\rangle+K^{*}(0,t)\left\langle\partial_{t}\delta\phi_{0}\>\delta\phi_{0}\right\rangle\bigg{)}\\\
&=\sum_{q\neq
0}\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{q}(t)\>\phi_{q}(t)\>\right\rangle
K^{*}(q,t)\\\ &\hskip
28.45274pt+\int_{0}^{\tau}dt\>\bigg{(}\mu(0)a^{*}(t)^{2}+\mu(0)a^{*}(t)\left\langle\partial_{t}\delta\phi_{0}(t)\right\rangle+\mu(0)a^{*}(t)K^{*}(0,t)\left\langle\delta\phi_{0}(t)\right\rangle+K^{*}(0,t)\left\langle\partial_{t}\delta\phi_{0}\>\delta\phi_{0}\right\rangle\bigg{)}\\\
&=\sum_{q\neq
0}\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{q}(t)\>\phi_{q}(t)\>\right\rangle
K^{*}(q,t)+\int_{0}^{\tau}dt\>\bigg{(}\mu(0)a^{*}(t)^{2}+K^{*}(0,t)\left\langle\partial_{t}\delta\phi_{0}\>\delta\phi_{0}\right\rangle\bigg{)}\\\
\end{split}$ (54)

Using $a^{*}(0,t)=(\bar{\phi}_{f}-\bar{\phi}_{i})/\tau\mu(0)$ and rearranging
terms leads to the simplified expression for heat:

$\begin{split}\langle\mathcal{Q}\rangle&=\sum_{q\neq
0}\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{q}(t)\>\phi_{q}(t)\>\right\rangle
K^{*}(q,t)+\int_{0}^{\tau}dt\>\bigg{(}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\mu(0)\tau^{2}}+K^{*}(0,t)\left\langle\partial_{t}\delta\phi_{0}\>\delta\phi_{0}\right\rangle\bigg{)}\\\
&=\sum_{q}\int_{0}^{\tau}dt\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\mu(q)\tau^{2}}\tau+\sum_{q\neq
0}\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{q}(t)\>\phi_{q}(t)\>\right\rangle
K^{*}(q,t)+\int_{0}^{\tau}dt\>\bigg{(}K^{*}(0,t)\left\langle\partial_{t}\delta\phi_{0}\>\delta\phi_{0}\right\rangle\bigg{)}\\\
&=\sum_{q}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\mu(q)\tau}+\sum_{q\neq
0}\int_{0}^{\tau}dt\>\left\langle\dot{\phi}_{q}(t)\>\phi_{q}(t)\>\right\rangle
K^{*}(q,t)+\int_{0}^{\tau}dt\>\bigg{(}K^{*}(0,t)\left\langle\partial_{t}\delta\phi_{0}\>\delta\phi_{0}\right\rangle\bigg{)}\\\
\end{split}$ (55)

Using eqs. 49 and 50 and Ito’s lemma, the expression for
$\langle\phi_{q}(t)^{2}K(q,t)\rangle$ reads:

$\begin{split}\frac{d}{dt}\left\langle\phi(q,t)^{2}K^{*}(q,t)\right\rangle&=\partial_{t}\left\langle\phi(q,t)^{2}\>K^{*}(q,t)\right\rangle\\\
&=2\left\langle\phi(q,t)\dot{\phi}(q,t)\right\rangle
K^{*}(q,t)+\dot{K}^{*}(q,t)\>\left\langle\phi(q,t)^{2}\right\rangle\end{split}$
(56)

The first term in the LHS of eq. 56 is used to simplify eq. 55, integration
over time leads to the following equation:

$\begin{split}\left\langle\mathbf{\mathcal{Q}}\right\rangle&=\sum_{q}\Bigg{\\{}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\mu(q)\tau}+\Big{[}\frac{K^{*}(q,t)}{2}\left\langle\phi_{q}(t)^{2}\right\rangle\Big{]}_{0}^{\tau}-\int_{0}^{\tau}dt\bigg{[}\frac{\dot{K}^{*}(q,t)}{2}\left\langle\phi_{q}(t)^{2}\right\rangle\bigg{]}\Bigg{\\}}\\\
&=\sum_{q}\Bigg{\\{}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\mu(q)\tau}+\Big{[}\frac{K^{*}(q,t)}{2}C(q,t)\Big{]}_{0}^{\tau}-\int_{0}^{\tau}dt\bigg{[}\frac{\dot{K}^{*}(q,t)}{2}C(q,t)\bigg{]}\Bigg{\\}}\\\
\end{split}$ (57)

$\left\langle\phi_{q}(t)^{2}\right\rangle=C(q,t)$, eq. 43b and eq. 48 are used
to simplify the eq. 57.

$\begin{split}\langle\mathbf{\mathcal{Q}}\rangle&=\sum_{q}\Bigg{\\{}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\mu(q)\tau}+\frac{\big{(}\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}\big{)}^{2}}{2\mu(q)\tau}-\int_{0}^{\tau}dt\bigg{[}-\frac{\big{(}\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}\big{)}^{2}}{2\mu(q)\tau^{2}}+\frac{T(\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)})}{\tau
C^{\frac{1}{2}}(q,t)}\bigg{]}\Bigg{\\}}\\\ \end{split}$ (58)

After doing basic mathematical steps like integration and rearrangement of
terms gives an equation for the average heat is given by eq. 59.

$\begin{split}\langle\mathbf{\mathcal{Q}}\rangle&=\sum_{q}\Bigg{\\{}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\mu(q)\tau}+\frac{\big{(}\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}\big{)}^{2}}{\mu(q)\tau}-\int_{0}^{\tau}dt\bigg{[}\frac{T\big{(}\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}\big{)}}{\tau
C^{\frac{1}{2}}(q,t)}\bigg{]}\Bigg{\\}}\\\
&=\sum_{q}\Bigg{\\{}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\mu(q)\tau}+\frac{\big{(}\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}\big{)}^{2}}{\mu(q)\tau}+\frac{T}{2}\log\bigg{(}\frac{C_{i}(q)}{C_{f}(q)}\bigg{)}\Bigg{\\}}\\\
\end{split}$ (59)

The form of eq. 59 is the same as the wasserstein distance computed using the
framework of optimal mass transport problem in eqs. 39, 40 and 41, the only
difference is that the $T\log(\sqrt{C_{i}(q)}/\sqrt{C_{f}(q)})$ term in eq. 59
gives the boundary term correction to heat which is the same as
$\Delta\mathcal{F}$ computed for the Gaussian probability distribution. The
boundary term does not depend on mobility $\mu(q,\phi)$ or driving time
$\tau$, signifying that the boundary term contribution $\Delta\mathcal{F}$ to
heat is independent of the model dynamics which should be the case. Similarly,
the first two terms are finite time contribution of the entropy production or
irreversible work. The contribution of the entropy production (irreversible
work) is equal to the wasserstein distance computed in eqs. 39, 40 and 41. For
model A, $\mu(q)=\mu$, thus heat can be written as follows:

$\begin{split}\langle\mathbf{\mathcal{Q}}\rangle&=\sum_{q}\Bigg{\\{}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\tau\mu}+\frac{\big{(}\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}\big{)}^{2}}{\tau\mu}+\frac{T}{2}\log\bigg{(}\frac{C_{i}(q)}{C_{f}(q)}\bigg{)}\Bigg{\\}}\\\
\end{split}$ (60)

For model B, $\mu(q)=\mu q^{2}$, this annihilates the dynamics corresponding
to the Zeroth Fourier mode. Using the decomposition for the heat, the heat is
given by:

$\begin{split}\langle\mathbf{\mathcal{Q}}\rangle&=\sum_{q\neq
0}\Bigg{\\{}\frac{\big{(}\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}\big{)}^{2}}{\tau\mu
q^{2}}+\frac{T}{2}\log\bigg{(}\frac{C_{i}(q)}{C_{f}(q)}\bigg{)}\Bigg{\\}}\\\
\end{split}$ (61)

eqs. 60 and 61 are consistent with eq. 40. In this section, we started with
formulating a minimisation problem to compute the minimum entropy production
for driving the field $\phi$ from the initial to the final phase. Using the
framework of the Optimal mass transport problem, we showed that the entropy
production ( irreversible work ) is equal to the Information theoretical
distance measure called ’Wasserstein distance’. Further using the consistency
condition of the optimal transport map, we computed the optimal driving
protocol. We used the optimal driving protocol to compute the average heat and
average work. We showed that the optimal driving protocol yields the
equivalence between the entropy production (irreversible work) and
’Wasserstein distance.’

###  II.4 Minimum work corresponding to the optimal protocol

Using the first law of thermodynamics ($\langle
W\rangle=\langle\mathcal{Q}\rangle-\langle\Delta\mathcal{H}\rangle$), we can
compute the average total work. Using eq. 57 and the first law of
thermodynamics, the equation for the total average work reads as follows:

$\begin{split}\langle
W\rangle&=\sum_{q}\Bigg{\\{}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\mu(q)\tau}+\bigg{[}\frac{\Delta_{t}K^{*}(q,t)}{2}C(q,t)\bigg{]}_{0}^{\tau}-\int_{0}^{\tau}dt\bigg{[}\frac{\dot{K}^{*}(q,t)}{2}C(q,t)\bigg{]}\Bigg{\\}}\\\
\end{split}$ (62)

Where we have used
$\langle\Delta\mathcal{H}\rangle=\frac{1}{2}k_{f}(q)C_{f}(q)-\frac{1}{2}k_{i}(q)C_{i}(q)$
such that $k_{i}(q)$ and $k_{f}(q)$ are the initial and final trap stiffness
of the $q^{th}$ Fourier modes. Thus $\Delta_{0}=K_{0^{+}}^{*}(q)-k_{i}(q)$ and
$\Delta_{\tau}=K_{\tau^{-}}^{*}(q)-k_{f}(q)$ are jumps at the initial and the
final point of the protocol. For a Gaussian equilibrium distribution,
$C_{i}(q)=1/k_{i}(q);C_{f}(q)=1/k_{f}(q)$, which leads to
$\langle\Delta\mathcal{H}\rangle=0$.

$\begin{split}\langle
W\rangle&=\sum_{q}\Bigg{\\{}\frac{(\bar{\phi}_{f}-\bar{\phi}_{i})^{2}\delta_{q,0}}{\tau\mu}+\frac{\big{(}\sqrt{C_{f}(q)}-\sqrt{C_{i}(q)}\big{)}^{2}}{\tau\mu}+\frac{T}{2}\log\bigg{(}\frac{C_{i}(q)}{C_{f}(q)}\bigg{)}\Bigg{\\}}\\\
\end{split}$ (63)

The first two terms correspond to the irreversible work which scales inversely
with driving time $\tau$. The third term is the quasistatic work contribution
$\Delta\mathcal{F}$.

The expression for the stochastic work of the $q^{th}$ Fourier mode defined in
the main text is given by:

$W_{q}=\Delta_{0}K\frac{\phi_{q}(0)^{2}}{2}+\Delta_{\tau}K\frac{\phi_{q}(\tau)^{2}}{2}-\int_{0}^{\tau}\partial_{t}K(q,t)\frac{\phi_{q}^{2}(t)}{2}$
(64)

Equivalently we could have started with the definition of the stochastic work
eq. 64, using $\langle W\rangle=\sum_{q}\langle W_{q}\rangle$ and
$C(q,t)=\langle\phi_{q}^{2}(t)\rangle$ would lead to eq. 62 and further
simplification to eq. 63.

##  III Analytical expression of wasserstein distance for model A dynamics

###  III.1 Infinite frequency cutoff

Building on the previous sections, in this section we aim to compute closed-
form expressions for the lower bound on entropy production (irreversible work)
for model A and B driving dynamics in a one-dimensional, two-dimensional and
three-dimensional space. Using the expressions for the variance of the Fourier
modes computed in section I.1 and the Wasserstein distance given by eq. 40.
The lower bound on entropy production for model A is given by:

$\begin{split}\frac{\mathcal{W}_{2}\Big{(}P[\bar{\phi_{f}}],P[\bar{\phi_{i}}]\Big{)}^{2}}{\tau\mu}&\geq\frac{V\parallel\bar{\phi_{f}}-\bar{\phi_{i}}\parallel^{2}}{\tau\mu}+\frac{VT}{\tau\mu}\int_{-\infty}^{\infty}\frac{dq}{2\pi}\Bigg{[}\frac{1}{\sqrt{\alpha_{f}+\kappa
q^{2}}}-\frac{1}{\sqrt{\alpha_{i}+\kappa q^{2}}}\Bigg{]}^{2}\\\
&\geq\frac{V\parallel\bar{\phi_{f}}-\bar{\phi_{i}}\parallel^{2}}{\tau\mu}+\frac{VT}{\tau\mu}\int_{-\infty}^{\infty}\frac{dq}{2\pi^{d}}\Bigg{[}\frac{1}{\alpha_{f}+\kappa
q^{2}}+\frac{1}{\alpha_{i}+\kappa q^{2}}\\\ &\hskip
85.35826pt-\frac{2}{\sqrt{(\alpha_{f}+\kappa q^{2})(\alpha_{i}+\kappa
q^{2})}}\Bigg{]}\end{split}$ (65)

After integrating eq. 65 wasserstein distance can be reduced to the following
form:

$\begin{split}\frac{\mathcal{W}_{2}\Big{(}P[\bar{\phi_{f}}],P[\bar{\phi_{i}}]\Big{)}^{2}}{\tau\mu}&\geq\frac{V\parallel\bar{\phi_{f}}-\bar{\phi_{i}}\parallel^{2}}{\tau\mu}+\frac{VT}{2\pi\tau\mu}\Bigg{[}\frac{\pi}{\sqrt{\kappa\alpha_{f}}}+\frac{\pi}{\sqrt{\kappa\alpha_{i}}}-\frac{4}{\sqrt{\kappa\alpha_{f}}}K\Bigg{(}\sqrt{1-\frac{\alpha_{i}}{\alpha_{f}}}\Bigg{)}\Bigg{]}\end{split}$
(66)

$\alpha_{i}$ and $\alpha_{f}$ are the dummy variables defined in section I.1.
$K(k)$ is the complete elliptical integral of the first kind given by the
following equation:

$K(k)=\int_{0}^{\frac{\pi}{2}}\frac{d\theta}{\sqrt{1-k^{2}sin\theta^{2}}}$
(67)

Hence Wasserstein distance is always positive. For $\alpha_{f}=\alpha_{i}$,
the second term does not contribute to the Wasserstein distance. Hence the
minimum Wasserstein distance is zero.

###  III.2 Finite frequency cutoff

A generalized expression for the Wasserstein distance in d-dimensional space
is given by the following equation.

$\begin{split}\frac{\mathcal{W}_{2}\Big{(}P[\bar{\phi_{f}}],P[\bar{\phi_{i}}]\Big{)}^{2}}{\tau}&\geq\frac{V\parallel\bar{\phi_{f}}-\bar{\phi_{i}}\parallel^{2}}{\tau}+\frac{VTS_{d-1}}{(2\pi)^{d}\>\tau}\int_{0}^{q_{0}}q^{d-1}dq\Bigg{[}\frac{1}{\alpha_{f}+\kappa
q^{2}}+\frac{1}{\alpha_{i}+\kappa q^{2}}\\\ &\hskip
142.26378pt-\frac{2}{\sqrt{(\alpha_{f}+\kappa q^{2})(\alpha_{i}+\kappa
q^{2})}}\Bigg{]}\end{split}$ (68)

$S_{d-1}$ corresponds to integration over a d-dimensional unit sphere. The
integration over the $q$ is divergent for $d>1$, hence we need to impose a
cutoff $q_{0}$ for the range of Fourier frequency modes.

####  III.2.1 One dimensional space

If Fourier frequencies in eq. 65 are confined to the range of values
$\\{-q_{0},q_{0}\\}$ instead of $\\{-\infty,\infty\\}$, the generalized form
of the eq. 66 can be written as follows:

$\begin{split}\frac{\mathcal{W}_{2}\Big{(}P[\bar{\phi_{f}}],P[\bar{\phi_{i}}]\Big{)}^{2}}{\tau}&\geq\frac{V\parallel\bar{\phi_{f}}-\bar{\phi_{i}}\parallel^{2}}{\tau}+\frac{VT}{2\pi\tau}\Bigg{[}\frac{2}{\sqrt{k\alpha_{f}}}\tan^{-1}{\bigg{(}\sqrt{\frac{\kappa}{\alpha_{f}}}q_{0}\bigg{)}}+\frac{2}{\sqrt{\kappa\alpha_{i}}}\tan^{-1}\bigg{(}{\sqrt{\frac{\kappa}{\alpha_{i}}}q_{0}}\bigg{)}\\\
&\hskip
85.35826pt-\frac{4}{\sqrt{\kappa\alpha_{f}}}F\Bigg{(}\tan^{-1}\bigg{(}q_{0}\sqrt{\frac{\kappa}{\alpha_{i}}}\bigg{)},\sqrt{1-\frac{\alpha_{i}}{\alpha_{f}}}\Bigg{)}\Bigg{]}\end{split}$
(69)

Where $F(\phi,k)$ is defined as the incomplete elliptical integral of the
first kind:

$F(\phi,k)=\int_{0}^{\phi}\frac{d\theta}{\sqrt{1-k^{2}sin\theta^{2}}}$ (70)

eq. 70 can be equivalently written as:

$F(\phi,k)=\int_{0}^{tan(\phi)}\frac{dq}{\sqrt{(1+q^{2})(1+k^{\prime
2}q^{2})}}$ (71)

Where $k^{\prime}$ in eq. 71 is given by the complementary elliptical modules
$k^{\prime 2}=1-k^{2}$. eq. 71 is can be used to obtain eqs. 66 and 69 in the
limit $q_{0}=\infty$.

####  III.2.2 Two dimensional space

For a two-dimensional space, the closed-form analytical expression for the
Wasserstein distance is given by:

$\begin{split}\frac{\mathcal{W}_{2}\Big{(}P[\bar{\phi_{f}}],P[\bar{\phi_{i}}]\Big{)}^{2}}{\tau}&\geq\frac{V\parallel\bar{\phi_{f}}-\bar{\phi_{i}}\parallel^{2}}{\tau}+\frac{VTS_{1}}{8\pi^{2}\>\tau\kappa}\Bigg{[}\log\bigg{(}1+\frac{\kappa
q_{0}^{2}}{\alpha_{f}}\bigg{)}+\log\bigg{(}1+\frac{\kappa
q_{0}^{2}}{\alpha_{i}}\bigg{)}\\\ &\hskip
28.45274pt-4\log\bigg{(}\frac{\sqrt{\alpha_{f}+\kappa
q_{0}^{2}}+\sqrt{\alpha_{i}+\kappa
q_{0}^{2}}}{\sqrt{\alpha_{f}}+\sqrt{\alpha_{i}}}\bigg{)}\Bigg{]}\end{split}$
(72)

By definition, $S_{1}$ is equal to $2\pi$.

####  III.2.3 Three dimensional space

Using eq. 68, the closed form analytical expression for Wasserstein distance
in a three-dimensional space is given by:

$\begin{split}\frac{\mathcal{W}_{2}\Big{(}P[\bar{\phi_{f}}],P[\bar{\phi_{i}}]\Big{)}^{2}}{\tau}&\geq\frac{V\parallel\bar{\phi_{f}}-\bar{\phi_{i}}\parallel^{2}}{\tau}+\frac{VTS_{2}}{8\pi^{3}\>\tau}\Bigg{[}\frac{2q_{0}}{\kappa}-\frac{\sqrt{\alpha_{i}}}{\kappa^{\frac{3}{2}}}\tan^{-1}\bigg{(}q_{0}\sqrt{\frac{\kappa}{\alpha_{i}}}\bigg{)}-\frac{\sqrt{\alpha_{f}}}{\kappa^{\frac{3}{2}}}\tan^{-1}\bigg{(}q_{0}\sqrt{\frac{\kappa}{\alpha_{f}}}\bigg{)}\\\
&+\frac{2\sqrt{\alpha_{i}}}{\kappa^{\frac{3}{2}}}E\Bigg{(}\tan^{-1}\bigg{(}q_{0}\sqrt{\frac{k}{\alpha_{i}}}\bigg{)},\sqrt{1-\frac{\alpha_{i}}{\alpha_{f}}}\Bigg{)}-\frac{q_{0}\sqrt{1+q_{0}^{2}\Big{(}1-\frac{\alpha_{i}}{\alpha_{f}}\Big{)}}}{\kappa^{\frac{3}{2}}(1+q_{0}^{2})}\Bigg{]}\end{split}$
(73)

Where $S_{2}=4\pi$.

##  IV Analytical expression of Wasserstein distance for model B dynamics

###  IV.1 Infinite frequency cutoff

Due to infrared divergence, the closed-form expression for Wasserstein
distance can not be computed. The expression for the Wasserstein distance in
eq. 62 is given by the following expression:

$\begin{split}\frac{\mathcal{W}_{2}\Big{(}P[\bar{\phi_{f}}],P[\bar{\phi_{i}}]\Big{)}^{2}}{\tau}&\geq\frac{VT}{\tau}\int_{-\infty}^{\infty}\frac{q^{d-3}dq}{(2\pi)^{d}}\Bigg{[}\frac{1}{\sqrt{\alpha_{f}+\kappa
q^{2}}}-\frac{1}{\sqrt{\alpha_{i}+\kappa q^{2}}}\Bigg{]}^{2}\end{split}$ (74)

###  IV.2 Finite frequency cutoff

Analogously we can generalize the Wasserstein distance for model B in a
d-dimensional space. The generalisation of Wasserstein distance for a
d-dimensional space gives the following equation.

$\begin{split}\frac{\mathcal{W}_{2}\Big{(}P[\bar{\phi_{f}}],P[\bar{\phi_{i}}]\Big{)}^{2}}{\tau}&\geq\frac{VTS_{d-1}}{(2\pi)^{d}\>\tau}\int_{\frac{q_{0}}{N}}^{q_{0}}q^{d-3}dq\Bigg{[}\frac{1}{\alpha_{f}+\kappa
q^{2}}+\frac{1}{\alpha_{i}+\kappa q^{2}}-\frac{2}{\sqrt{(\alpha_{f}+\kappa
q^{2})(\alpha_{i}+\kappa q^{2})}}\Bigg{]}\end{split}$ (75)

Where $S_{d-1}$ is given by integration over a d-dimensional unit sphere.
Integration over $q$ is divergent, hence we need to impose a cutoff
$[q_{i},q_{f}]$ such that $[q_{i},q_{f}]=[\pi/N,\pi]$ for the range of Fourier
frequency modes.

####  IV.2.1 Two dimensional space

For a two-dimensional space, the analytical expression for the Wasserstein
distance is given by the following equation:

$\begin{split}\frac{\mathcal{W}_{2}\Big{(}P[\bar{\phi_{f}}],P[\bar{\phi_{i}}]\Big{)}^{2}}{\tau}&\geq\frac{VTS_{1}}{8\pi^{2}\>\tau}\Bigg{[}\frac{1}{\alpha_{f}}\log\bigg{(}\frac{\alpha_{f}N^{2}+\kappa\pi^{2}}{\alpha_{f}+\kappa\pi^{2}}\bigg{)}+\frac{1}{\alpha_{i}}\log\bigg{(}\frac{\alpha_{i}N^{2}+\kappa\pi^{2}}{\alpha_{i}+\kappa\pi^{2}}\bigg{)}\vspace{1cm}\\\
&\hskip
28.45274pt-\frac{2}{\sqrt{\alpha_{i}\alpha_{f}}}\log\bigg{(}\frac{\sqrt{\alpha_{i}(\alpha_{f}+\kappa
q_{0}^{2})}-\sqrt{\alpha_{f}(\alpha_{i}+\kappa
q_{0}^{2})}}{\sqrt{\alpha_{i}(\alpha_{f}+\kappa
q_{0}^{2})}+\sqrt{\alpha_{f}(\alpha_{i}+\kappa
q_{0}^{2})}}\bigg{)}\Bigg{|}_{\frac{\pi}{N}}^{\pi}\Bigg{]}\\\
&=\frac{VTS_{1}}{8\pi^{2}\>\tau}\Bigg{[}\frac{1}{\alpha_{f}}\log\bigg{(}\frac{\alpha_{f}N^{2}+\kappa\pi^{2}}{\alpha_{f}+\kappa\pi^{2}}\bigg{)}+\frac{1}{\alpha_{i}}\log\bigg{(}\frac{\alpha_{i}N^{2}+\kappa\pi^{2}}{\alpha_{i}+\kappa\pi^{2}}\bigg{)}\vspace{1cm}\\\
&\hskip
14.22636pt+\frac{4}{\sqrt{\alpha_{i}\alpha_{f}}}\tanh^{-1}\bigg{(}\frac{\alpha_{i}(\alpha_{f}+\kappa\pi^{2})}{\alpha_{f}(\alpha_{i}+\kappa\pi^{2})}\bigg{)}-\frac{4}{\sqrt{\alpha_{i}\alpha_{f}}}\tanh^{-1}\bigg{(}\frac{\alpha_{i}(\alpha_{f}N^{2}+\kappa\pi^{2})}{\alpha_{f}(\alpha_{i}N^{2}+\kappa\pi^{2})}\bigg{)}\Bigg{]}\\\
&=\frac{VTS_{1}}{8\pi^{2}\>\tau}\Bigg{[}\frac{1}{\alpha_{f}}\log\bigg{(}\frac{\alpha_{f}N^{2}+\kappa\pi^{2}}{\alpha_{f}+\kappa\pi^{2}}\bigg{)}+\frac{1}{\alpha_{i}}\log\bigg{(}\frac{\alpha_{i}N^{2}+\kappa\pi^{2}}{\alpha_{i}+\kappa\pi^{2}}\bigg{)}\vspace{1cm}\\\
&\hskip
28.45274pt+\frac{2}{\sqrt{\alpha_{i}\alpha_{f}}}\log\Bigg{(}\frac{\sqrt{1+\frac{\kappa\pi^{2}}{\alpha_{i}}}-\sqrt{1+\frac{\kappa\pi^{2}}{\alpha_{f}}}}{\sqrt{1+\frac{\kappa\pi^{2}}{\alpha_{i}}}+\sqrt{1+\frac{\kappa\pi^{2}}{\alpha_{f}}}}\Bigg{)}-\frac{2}{\sqrt{\alpha_{i}\alpha_{f}}}\log\Bigg{(}\frac{\sqrt{1+\frac{\kappa\pi^{2}}{\alpha_{i}N^{2}}}-\sqrt{1+\frac{\kappa\pi^{2}}{\alpha_{f}N^{2}}}}{\sqrt{1+\frac{\kappa\pi^{2}}{\alpha_{i}N^{2}}}+\sqrt{1+\frac{\kappa\pi^{2}}{\alpha_{f}N^{2}}}}\Bigg{)}\Bigg{]}\end{split}$
(76)

##  V Boundary term for the heat

The boundary term for heat dissipation is given by the following equation:

$\begin{split}T\Delta
S&=\sum_{\mathbf{q}}T\log\bigg{(}\frac{\sqrt{C_{f}(q)}}{\sqrt{C_{i}(q)}}\bigg{)}=\sum_{q}\frac{T}{2}\log\bigg{(}\frac{C_{f}(q)}{C_{i}(q)}\bigg{)}\\\
&=\sum_{q}\frac{T}{2}\log\bigg{(}\frac{\alpha_{i}+\kappa
q^{2}}{\alpha_{f}+\kappa
q^{2}}\bigg{)}=\frac{VT}{2}\int_{0}^{\pi}\frac{q^{d-1}dq}{(2\pi)^{d}}\log\bigg{(}\frac{\alpha_{i}+\kappa
q^{2}}{\alpha_{f}+\kappa q^{2}}\bigg{)}\end{split}$ (77)

###  V.1 Two dimensional space

For a two-dimensional space, the expression can be generalized to the
following form:

$\begin{split}T\Delta
S&=\frac{VTS_{d-1}}{2}\int_{0}^{\pi}\frac{q^{d-1}\>dq}{(2\pi)^{d}}\log\bigg{(}\frac{\alpha_{i}+\kappa
q^{2}}{\alpha_{f}+\kappa q^{2}}\bigg{)}\\\
&=\frac{VT\>2\pi}{8\pi^{2}}\int_{0}^{\pi}dq\>q\log\bigg{(}\frac{\alpha_{i}+\kappa
q^{2}}{\alpha_{f}+\kappa q^{2}}\bigg{)}\\\
&=\frac{VT}{4\pi}\int_{0}^{\pi}dq\>q\log\bigg{(}\frac{\alpha_{i}+\kappa
q^{2}}{\alpha_{f}+\kappa q^{2}}\bigg{)}\\\ &=\frac{VT}{8\pi
k}\Bigg{[}\alpha_{i}\log\bigg{(}1+\frac{\kappa\pi^{2}}{\alpha_{i}}\bigg{)}-\alpha_{f}\log\bigg{(}1+\frac{\kappa\pi^{2}}{\alpha_{f}}\bigg{)}+\kappa\pi^{2}\log\bigg{(}\frac{\alpha_{i}+\kappa\pi^{2}}{\alpha_{f}+\kappa\pi^{2}}\bigg{)}\Bigg{]}\end{split}$
(78)

##  VI Pareto optimization with mean and standard deviation of the work

###  VI.1 Single Fourier mode case

The Langevin equation for the $q^{th}$ Fourier mode is given by the following
equation:

$\begin{split}\partial_{t}\phi_{q}(t)&=\mu(q,t)f_{q}(t)+\sqrt{2T\mu(q,t)}\>\Lambda_{q}(t)\\\
&=\mu(q,t)K(q,t)\phi_{q}(t)+\sqrt{2T\mu(q,t)}\>\Lambda_{q}(t)\end{split}$ (79)

The framework of optimal protocol discussed so far concerns optimising the
average work for driving. One could instead try to optimize a linear
combination of the mean work $\langle W_{q}\rangle$ and the variance of the
work, $\sqrt{\langle\big{(}W_{q}-\langle
W_{q}\rangle\big{)}^{2}\rangle}=\sigma_{q}$ as previously done in [52]. The
linear cost function for the optimisation problem is defined as follows:

$\mathcal{J}^{\lambda}_{q}=\lambda\langle W_{q}\rangle+(1-\lambda)\sigma_{q}$
(80)

Where $\lambda$ in eq. 80 characterises weight given to the mean or
fluctuations of the work. $\lambda=1$ corresponds to minimising the average
work which yields the exact analytical results discussed in previous sections.
On the other hand, $\lambda=1$ corresponds to minimising the work fluctuations
only. The optimisation problem with eq. 80 is reduced to computing the value
of $\langle W_{q}\rangle$ and $\sigma_{q}$ that minimises cost function
$\mathcal{J}^{\lambda}_{q}$. This is addressed by numerically optimizing the
optimum protocol that minimizes the cost function $\mathcal{J}^{\lambda}_{q}$
in eq. 80. For any given protocol $K(q,t)$, computing the values of $\langle
W_{q}\rangle$ and $\sigma_{q}$ is equivalent to solving the following set of
ODE’s:

$\begin{split}&\partial_{t}\langle\omega_{q}(t)\rangle=-\frac{\mu(q)\dot{K}(q,t)}{2}\langle\phi_{q}(t)^{2}\rangle\\\
&\partial_{t}\langle\phi_{q}(t)^{2}\rangle=2\mu(q)K(q,t)\langle\phi_{q}(t)^{2}\rangle+2T\mu(q)\\\
&\partial_{t}\langle\phi_{q}(t)^{2}\omega_{q}(t)\rangle=2T\mu(q)\langle\omega_{q}(t)\rangle-\frac{1}{2}\mu(q)\dot{K}(q,t)\langle\phi_{q}(t)^{4}\rangle+2\mu(q)K(q,t)\langle\phi_{q}(t)^{2}\omega_{q}(t)\rangle\\\
&\partial_{t}\langle\omega_{q}(t)^{2}\rangle=-\mu(q)\dot{K}(q,t)\langle\phi_{q}(t)^{2}\omega_{q}(t)\rangle\\\
&\partial_{t}\langle\phi_{q}(t)^{4}\rangle=4\mu(q)K(q,t)\langle\phi_{q}(t)^{4}\rangle+12T\mu(q)\langle\phi_{q}(t)^{2}\rangle\\\
&\partial_{t}\langle\phi_{q}(t)^{2}\phi_{q}(0)^{2}\rangle=2\mu(q)K(q,t)\langle\phi_{q}(t)^{2}\phi_{q}(0)^{2}\rangle+2T\mu(q)\langle\phi_{q}(0)^{2}\rangle\\\
&\partial_{t}\langle\phi_{q}(0)^{2}\omega_{q}(t)\rangle=-\frac{1}{2}\mu(q)\dot{K}(q,t)\langle\phi_{q}(t)^{2}\phi_{q}(0)^{2}\rangle\end{split}$
(81)

The set of eq. 81 is a generalisation of the previously proposed equations in
[52]. For Model A dynamics $\mu(q)=\mu$. For Model B dynamics, $\mu(q)=\mu
q^{2}$. Where the work done for a single noise realization of one protocol is
written as:

$W_{q}=\omega_{q}(\tau)+\frac{\Delta_{\tau}K(q,\tau)}{2}\phi_{q}(\tau)^{2}-\frac{\Delta_{0}K(q,0)}{2}\phi_{q}(0)^{2};\hskip
28.45274pt\omega_{q}(\tau)=-\int_{0}^{\tau}\frac{\dot{K}(q,t)}{2}\phi_{q}(t)^{2}\>dt$
(82)

$\Delta_{\tau}=K_{f}(q)-K_{\tau^{-}}(q)$ and
$\Delta_{0}=K_{0^{+}}(q)-K_{i}(q)$ are defined as the jumps at the final and
initial point. $q^{th}$ Fourier mode satisfies eq. 79 for driving from a
initial $K_{i}(q)=-(\alpha_{i}+\kappa q^{2})$ to a final trap stiffness
$K_{f}(q)=-(\alpha_{f}+\kappa q^{2})$. Taking the average over noise
realization of the work, we write the equations for the first and second
moments of the work:

$\langle
W_{q}\rangle=\langle\omega_{q}(\tau)\rangle+\frac{\Delta_{\tau}K(q,\tau)}{2}\langle\phi_{q}(\tau)^{2}\rangle-\frac{\Delta_{0}K(q,0)}{2}\langle\phi_{q}(0)^{2}\rangle$
(83a) $\begin{split}\langle
W_{q}^{2}\rangle&=\langle\omega_{q}(\tau)^{2}\rangle+\Delta_{\tau}K(q,\tau)\langle\phi_{q}(\tau)^{2}\omega_{q}(\tau)\rangle-\Delta_{0}K(q,0)\langle\phi_{q}(0)^{2}\omega_{q}(\tau)\rangle\\\
&+\frac{1}{4}\Big{[}\Delta_{0}K(q,0)^{2}\langle\phi_{q}(0)^{4}\rangle+\Delta_{\tau}K(q,\tau)^{2}\langle\phi_{q}(\tau)^{4}\rangle-2\Delta_{0}K(q,0)\>\Delta_{\tau}K(q,\tau)\langle\phi_{q}(\tau)^{2}\>\phi_{q}(0)^{2}\rangle\Big{]}\end{split}$
(83b)

The set of eqs. 81, 82 and 83 are simulated using Monte-Carlo simulation to
find an optimum protocol that minimises the cost function eq. 80.

###  VI.2 Multiple Fourier mode case

In the previous section, we discussed Pareto optimisation of the mean and
variance for each independent Fourier mode. The total work done $W$ is a
stochastic variable given by the sum of each independent Fourier mode.

$W=\sum_{q}W_{q}$ (84)

The mean and variance of the stochastic variable $W$ satisfy the following
equations:

$\langle W\rangle=\sum_{q}\langle W_{q}\rangle$ (85a)
$\sigma=\sqrt{\langle\big{(}W-\langle
W\rangle\big{)}^{2}\rangle}=\sqrt{\sum_{q}\sigma_{q}^{2}}$ (85b)

The cost function eq. 80 is modified to the following equation:

$\mathcal{J}^{\lambda}=\lambda\langle
W\rangle+(1-\lambda)\sigma=\lambda\sum_{q}\langle
W_{q}\rangle+(1-\lambda)\sqrt{\sum_{q}\sigma_{q}^{2}}$ (86)

Each Fourier mode evolves independently as previously given by eqs. 81, 82 and
83. The sole difference for the multi-mode optimisation is the cost function
eq. 86 has a dependence on all Fourier modes and cannot be decoupled to an
independent optimisation problem for each Fourier mode. The dynamics of each
Fourier mode are still decoupled. Hence the mean and variance of the work
distribution is computed using the set of eqs. 81, 82 and 83. We can decouple
the total cost function as a sum of independent cost functions by redefining
$\mathcal{J}^{\lambda}_{q}$ as a linear combination of average work and
variance of work, refer to the main text. Importantly, the structure of the
optimisation problem for a linear combination of average work and standard
deviation of the work is preserved. Thus preserving the single Fourier mode
phase transition to multiple Fourier mode case as well.

##  VII Numerical scheme and parameters

We used the Monte-Carlo simulation and discretisation to solve Pareto
optimisation as previously done in [52]. The ODEs are integrated by
discretizing the protocol into 100 points with linear interpolations. As a
sanity check, our scheme was able to reproduce the results obtained in [52].
At each Monte-Carlo step, the scaled protocol was stochastically updated by
the following equation
$\mu(q)K_{n}(q,t)\rightarrow\mu(q)\big{(}K_{n}(q,t)+l_{r}\>max(K_{i},K_{f})\>\epsilon(t)\big{)}$.
The protocols minimising the cost function were accepted for the next Monte-
Carlo step. $K_{n}(q,t)$ denotes the $n^{th}$ point of the discreatised
protocol $K(q,t)$. $l_{r}$ denotes the learning rate for stochastic gradient
descent. $\epsilon(t)$ generates a random number from a Gaussian distribution
with zero mean and unit standard deviation which randomly changes the $n^{th}$
point of the protocol $K(q,t)$, $\epsilon(t)\in\mathcal{N}(0,1)$. The
$l_{r}=0.005$ was used to simulate $3*10^{5}$ Monte-Carlo steps. $10^{5}$
hysteresis steps were used to characterize the observed first-order phase
transition.

##  VIII Other plots

Figure 5: (a) For small values of $\mu$, there exists a similar kind of change
of slope for cost function $\mathcal{J}^{0}_{\lambda}$. Unlike higher values
of $\mu$ a discontinuity in the $d\mathcal{J}^{0}_{\lambda}/d\lambda$ is
missing hence this is not characterised as a phase transition. It can be seen
in the plot the change of slope is rather continuous and hysteresis is missing
too. (b) The protocols are plotted for different values of $\lambda$ at the
higher values of $\mu$. Two different classes of protocols corresponding to
mean $\lambda=1$ (topmost red line) and variance dominating $\lambda=0$
(bottom-most blue line) are observed. For the intermediate values of
$\lambda$, the protocols are denoted by interpolation between the red and blue
colours. At $\lambda_{c}$, a sudden jump is observed for the optimal protocol
which corresponds to a first-order phase transition quoted in the main text.
The variance optimizing class of protocols ( $\lambda=0$) are characterized by
the uniform interpolation from the initial to the final state. We observe that
they exhibit the property $K^{*}(t)=min\\{K_{i},K_{f}\\}$. Through ad-hoc
analysis of eqs. 81 and 82, it can be observed that the variance of the work
scales as $K^{3}(t)$ to the leading order, in comparison to the mean work
scales as $K^{2}(t)$. A more systematic theoretical analysis needs to be
deployed to understand this interplay. We are unaware of closed-form
expressions for the mean and variance of the work for the systems operating
far away from the linear response regime. (c) Numerically obtained phase
diagram $\mu-\lambda$ for Model A is plotted for different Fourier modes.
Increasing the Fourier frequency $q$, increases the critical value of mobility
$\mu_{c}$ to observe a phase transition. (d) Numerically obtained phase
diagram $\mu-\lambda$ for Model B is plotted for different Fourier modes.
Model B with mobility $\mu$ is equivalent to model A with $\mu^{{}^{\prime}}$
such that $\mu q^{2}=\mu^{{}^{\prime}}$. In other words, fig(d) is equivalent
to plotting fig (c) with a scaled $\mu$ axis so that scaling is different for
different Fourier modes.

##  IX Parameters used for plots in the maintext

Parameters for plots in fig.(2) of the main text:
(a) $a_{i}=1.0$, $c_{i}=0.0$, $a_{f}=-1.0$, $c_{f}=2.12$, $\phi_{i}=0.0$,
$\phi_{f}=2.52$, $b=1.0$, $\kappa=0.5$, $N=64$ lattice sites with $\Delta
x=1$, which implies $\alpha_{i}=1.0$, $\alpha_{f}=7.34$.
(b), (c), (d), (e) $\mu=4$, $\alpha_{i}=5$, $\alpha_{f}=0.2$, $q=0$, $\tau=1$
(f) $\mu=7$, $\alpha_{i}=5$, $\alpha_{f}=0.2$, $q=l\pi/N$, $\tau=1$, $N=16$