\section{Derivation of the Thermodynamic Metric}

We consider a system with coordinates $\boldsymbol{x}\in \mathbb{R}} \def\NN{\mathbb{N}} \def\ZZ{\mathbb{Z}^d$ subject to an ergodic Markovian dynamics.
A particle trajectory in which the $i$th particle is at $\boldsymbol{x}^{(i)}_0$ at $t=0$ is denoted $\boldsymbol{X}^{(i)}(t; \boldsymbol{x}_0^{(i)})$.
Within this framework, the density
\begin{equation}
  \rho(\boldsymbol{x}, t) = \avg{\delta(\boldsymbol{X}(t;\boldsymbol{x}_0)-\boldsymbol{x})}
\end{equation}
evolves according to a linear partial differential equation
\begin{equation}
  \partial_t \rho(\boldsymbol{x},t) = \mathcal{L}^{\dagger} \rho(\boldsymbol{x},t).
\end{equation}
For example, we consider the microscopic dynamics specified by the stochastic differential equation
\begin{equation}
  d\boldsymbol{X}(t) = -\nabla U(\boldsymbol{X}(t), \boldsymbol{\lambda}(t)) dt + \sqrt{2\beta^{-1}} d\boldsymbol{W}_t,
\end{equation}
where $\boldsymbol{\lambda}(t) \in \mathbb{R}} \def\NN{\mathbb{N}} \def\ZZ{\mathbb{Z}^k$ is a vector of control parameters conjugate to the generalized forces $\boldsymbol{\Theta}(\boldsymbol{x},t) = -\partial_{\boldsymbol{\lambda}} U(\boldsymbol{x}, \boldsymbol{\lambda}(t))$, $\beta$ is the inverse temperature, and $U$ is the potential energy of system, and $\boldsymbol{W}$ is a standard Wiener process.
In this case, the distribution satisfies the evolution equation
\begin{equation}
  \partial_t \rho(\boldsymbol{x},t) = -\nabla \cdot (\boldsymbol{v}(\boldsymbol{x},t) \rho(\boldsymbol{x},t))
  \label{eq:afp}
\end{equation}
where the local velocity $\boldsymbol{v}$ is explicitly
\begin{equation}
  \boldsymbol{v}(\boldsymbol{x}, t) = -\nabla U(\boldsymbol{x}, \boldsymbol{\lambda}(t)) - \beta^{-1} \nabla \log \rho(\boldsymbol{x},t)
\end{equation}
meaning that
\begin{equation}
  \mathcal{L}^\dagger = -\nabla U \cdot \nabla + \beta^{-1} \Delta.
\end{equation}

In this work, we focus on the problem of transforming an initial distribution, which could be the equilibrium distribution with the initial control parameters, e.g.,
\begin{equation}
  \rho_A(\boldsymbol{x}) = Z^{-1}(\boldsymbol{\lambda}(0)) e^{-\beta U(\boldsymbol{x},\boldsymbol{\lambda}(0))},
\end{equation}
into a target distribution $\rho_B(\boldsymbol{x})$ which is the distribution that the system relaxes to when the control parameters are fixed at $\boldsymbol{\lambda}(t_{\rm f})$.
The protocol has a finite duration $0<t_{\rm f}< \infty$.
More general distributions $\rho_A$ and $\rho_B$ could also be considered.
We employ the standard stochastic thermodynamics definitions of work,
\begin{equation}
  W[\boldsymbol{X}(t)] = -\int_0^{t_{\rm f}} \boldsymbol{\Theta}(\boldsymbol{X}(t,\boldsymbol{x}_0)) \cdot \dot{\boldsymbol{\lambda}}(t) dt,
  \label{eq:work}
\end{equation}
and heat,
\begin{equation}
  \mathcal{Q}[\boldsymbol{X}(t)] = -\int_0^{t_{\rm f}} \nabla U(\boldsymbol{X}(t,\boldsymbol{x}_0)) \circ d\boldsymbol{X}(t).
  \label{eq:heat}
\end{equation}
In this formulation, the work and the heat yield the total change in energy when the system is driven by the protocol, meaning that the first law of thermodynamics is simply,
\begin{equation}
    \Delta U = W - \mathcal{Q}.
\end{equation}
Using these expressions and \eqref{eq:afp}, the total entropy change is given by~\cite{aurell_refined_2012},
\begin{equation}
  \Delta \Sigma_{\rm tot} = \int_0^{t_{\rm f}} \beta(t) \int_{\Omega} \boldsymbol{v}(\boldsymbol{x},t)^T \boldsymbol{v}(\boldsymbol{x},t) \rho(\boldsymbol{x},t) d\boldsymbol{x}\ dt
  \label{eq:astot}
\end{equation}
where the excess entropy dissipated to the environment is given by the average dissipated heat
\begin{equation}
    \Delta \Sigma = \beta(\Delta U - \avg{W}) = \beta \avg{\mathcal{Q}}.
\end{equation}
The entropy change of the system is a state function and has no explicit dependence on the protocol used.

To find a minimum dissipation protocol $\boldsymbol{\lambda}_*$, we must solve the following minimization problem:
\begin{equation}
  \boldsymbol{\lambda}_* = \argmin_{\boldsymbol{\lambda}:[0,t_{\rm f}]\to \mathbb{R}} \def\NN{\mathbb{N}} \def\ZZ{\mathbb{Z}^k} \Delta \Sigma_{\rm tot}[\boldsymbol{\lambda}] \textrm{ subj. to } \rho(\cdot,0) = \rho_A,\quad \rho(\cdot,t_{\rm f}) = \rho_B.
  \label{eq:bb}
\end{equation}
As discussed in the main text, this is the Benamou-Brenier formulation of the optimal transport problem~\cite{benamou_monge_2002}.
It should be noted that the rigorous formulation of the relation between the minimizer \eqref{eq:bb} and the optimal transportation plan requires that the time-dependent local velocity fields and density $(\boldsymbol{v}, \rho)$ are allowed to vary over all possible functions that satisfy~\eqref{eq:fp} and the boundary conditions.
While we place no restrictions on the control parameters $\boldsymbol{\lambda}$, in practice we represent $\boldsymbol{\lambda}$ with a neural network, so we do restrict the function class to an extent, though sufficiently wide neural networks should allow us to find expressive enough protocols to ensure that this restriction is unimportant.

We do not know the distribution $\rho(\boldsymbol{x},t)$ at all intermediate times without solving~\eqref{eq:afp} because, as we change the protocol, the system is driven away from equilibrium.
However, for sufficiently weak driving, we can solve the equation for the dynamics perturbatively.
We first assume that the rate of driving is slow; explicitly, $\dot{\boldsymbol{\lambda}} \equiv \epsilon \dot{\boldsymbol{\lambda}}_\epsilon$ so we can write the adjoint linear operator driving the dynamics as
\begin{equation}
  \mathcal{L}^{\dagger} = \mathcal{L}_0^{\dagger} + \epsilon \mathcal{L}_1^{\dagger}
\end{equation}
with
\begin{equation}
  \mathcal{L}_1^\dagger = \dot{\boldsymbol{\lambda}}_\epsilon \cdot \partial_{\boldsymbol{\lambda}}.
\end{equation}
Expanding the solution to first order in $\epsilon$,
\begin{equation}
  \rho(\boldsymbol{x},t) = \rho_0(\boldsymbol{x},t) + \epsilon \rho_1(\boldsymbol{x},t) + \mathcal{O}(\epsilon^2),
\end{equation}
we see that
\begin{equation}
  \rho_0(\boldsymbol{x},t) = e^{-\beta U(\boldsymbol{x},\boldsymbol{\lambda}(t))}/Z(\boldsymbol{\lambda}(t))
\end{equation}
and the correction term satisfies
\begin{equation}
  \partial_t \rho_1(\boldsymbol{x},t) = \mathcal{L}_0^\dagger \rho_1(\boldsymbol{x},t) + \mathcal{L}_1^\dagger \rho_0(\boldsymbol{x},t)
\end{equation}
at order $\epsilon$.
Using the method of characteristics and the initial condition $\rho_1(\boldsymbol{x},0) = 0$ and assuming that the protocol is advanced from time $t$, we obtain
\begin{equation}
  \rho_1(\boldsymbol{x},t) = \int_0^\tau e^{\mathcal{L}_0^\dagger(\tau-s)} \dot{\boldsymbol{\lambda}}_\epsilon(t) \cdot \partial_{\boldsymbol{\lambda}(t)} \rho_0(\boldsymbol{x}^{\boldsymbol{\lambda}}(s),t) ds.
\end{equation}
Here $\boldsymbol{x}^{\boldsymbol{\lambda}}$ denotes the fast-timescale relaxation of the system to slow driving and $\boldsymbol{x}^{\boldsymbol{\lambda}}(0) = \boldsymbol{x}_0$; the protocol is effectively fixed over the relaxation time of $\boldsymbol{x}^{\boldsymbol{\lambda}}$ due to the assumed timescale separation.
In the integral, $\tau$ denotes the relaxation of the fast process, which typically can be taken to infinity assuming that the integral converges sufficiently quickly.
Writing the equilibrium distribution as
\begin{equation}
  \rho_0(\boldsymbol{x},t) = e^{-\beta U(\boldsymbol{x},\boldsymbol{\lambda}(t)) + \beta F(\boldsymbol{\lambda}(t))},
\end{equation}
where $F$ is the equilibrium free energy with the control parameters fixed,
\begin{equation}
  F(\boldsymbol{\lambda}) = -\beta^{-1} \log \int_{\Omega} e^{-\beta U(\boldsymbol{x},\boldsymbol{\lambda})} d\boldsymbol{x}.
\end{equation}
The derivative with respect to $\boldsymbol{\lambda}$ yields the deviation of the generalized forces from their average value at a particular point $\boldsymbol{x}$, that is, 
\begin{equation}
  \partial_{\boldsymbol{\lambda}(t)} \rho_0(\boldsymbol{x},t) = \beta\left( \boldsymbol{\Theta}(\boldsymbol{x}(s),t) - \avg{\boldsymbol{\Theta}(\boldsymbol{x}(0),t)}_{\boldsymbol{\lambda}(t)} \right) \rho_0(\boldsymbol{x},t).
\end{equation}
Because the protocol advances slowly relative to the timescale of relaxation of the system, we can extend $\tau\to \infty$ so that
\begin{equation}
  \epsilon \rho_1(\boldsymbol{x},t) = \epsilon \beta \rho_0(\boldsymbol{x}, t) \int_0^\infty \delta \boldsymbol{\Theta}(\boldsymbol{x}^{\boldsymbol{\lambda}}(s),t) \cdot \dot \boldsymbol{\lambda}(t) ds.
\end{equation}
The excess dissipation associated with the protocol can be computed by computing the average work associated with the perturbation.
The contribution from the total change in energy does not depend on the protocol and contributes the leading term, so we are left with, at order $\epsilon$, 
\begin{equation}
\begin{aligned}
      -\avg{W} &= \int_\Omega \int_0^{t_{\rm f}} \dot{\boldsymbol{\lambda}}(t) \cdot \partial_{\boldsymbol{\lambda}} U(\boldsymbol{x}, \boldsymbol{\lambda}(t)) \rho(\boldsymbol{x}, t)\ d\boldsymbol{x}\ dt \\
      & = \beta \int_{\Omega} \int_0^{t_{\rm f}} \dot{\boldsymbol{\lambda}}^T(t)\left( \int_0^\infty \delta \boldsymbol{\Theta}(\boldsymbol{x}^{\boldsymbol{\lambda}}(s),t) \delta \boldsymbol{\Theta}^T(\boldsymbol{x}^{\boldsymbol{\lambda}}_0,t) ds \right) \rho_0(\boldsymbol{x}^{\boldsymbol{\lambda}}_0,t)  \dot{\boldsymbol{\lambda}}(t)\ d\boldsymbol{x}^{\boldsymbol{\lambda}}_0\ dt.
\end{aligned}
\end{equation}
This expression optimizes \eqref{eq:astot} and recovers the thermodynamic metric of~\citet{sivak_thermodynamic_2012}.
The dissipation along a protocol $\boldsymbol{\lambda}:[0,t_{\rm f}]$ is proportional to the path functional
\begin{equation}
    L[\boldsymbol{\lambda}]:= \int_0^{t_{\rm f}} \dot{\boldsymbol{\lambda}}^T(t) \zeta(\boldsymbol{\lambda}(t)) \dot{\boldsymbol{\lambda}}(t) dt,
\end{equation}
where 
\begin{equation}
    \zeta(\boldsymbol{\lambda}(t)) = \beta \int_0^\infty \avg{\delta \boldsymbol{\Theta}(\boldsymbol{x}(s),t) \delta \boldsymbol{\Theta}^T(\boldsymbol{x}(0),t)}_{\boldsymbol{\lambda}} ds;
\end{equation}
the subscript $\boldsymbol{\lambda}$ denotes that the average for the initial condition is taken from the equilibrium distribution with the control parameters fixed at $\boldsymbol{\lambda}(t)$.
Because $\zeta$ is a positive semi-definite form, it is a Riemannian metric and all associated path lengths are non-negative.

\section{Computational Details for  Protocol Optimization}
We consider two different approaches for learning protocols that minimize \eqref{eq:stot}: we either specify target intermediate distributions and learn a protocol $\boldsymbol{\lambda}_*$ that drives the system through these distributions or specify a protocol $\boldsymbol{\lambda}$ and determine $\phi_*$, the speed at which the protocol is traversed.

We use automatic differentiation to learn a time-dependent $\boldsymbol{\lambda}_*$ that drives the system through target intermediate distributions. Here, $\boldsymbol{\lambda}_*:\mathbb{R}} \def\NN{\mathbb{N}} \def\ZZ{\mathbb{Z}^1 \to \mathbb{R}} \def\NN{\mathbb{N}} \def\ZZ{\mathbb{Z}^2$ is a single hidden layer neural network. The input to $\boldsymbol{\lambda}_*$ is the current time of the protocol and the output is the temperature and the mechanical parameter being modulated. Our general algorithm is shown below in Algorithm~\ref{alg:autodiff}.

\begin{algorithm}[H]
    \caption{Automatic Differentiation to Optimize $\boldsymbol{\lambda}_*$} \label{alg:autodiff}
    \begin{algorithmic}
    \State {\bfseries Initialize protocol $\boldsymbol{\lambda}_*$}
    \For {$e = 1, \dots, n_{\rm{epochs}}$}
    \State Initialize state ${\boldsymbol{x}_0}$
        \For {$t < \tau$} 
            \State{$\boldsymbol{X}^{\boldsymbol{\lambda}}(t+ \Delta t) = \boldsymbol{X}^{\boldsymbol{\lambda}}(t) + \dot{\boldsymbol{X}}^{\boldsymbol{\lambda}}(t) \Delta t$} \Comment{Use desired integration scheme}
            \State{$t \leftarrow t + \Delta t$}
            \If{end of interval}
                \State{Compute $\mathfrak{L}(t)$ and update $\boldsymbol{\lambda}_*$}
                \State{Clear gradient information}
            \EndIf
        \EndFor
    \EndFor
    \end{algorithmic}
\end{algorithm}

    
The automatic differentiation approach introduced above can be extended to any integration scheme that stores the exact gradients along the integration step. Alternatively, an adjoint method can be used, where gradients are computed by integrating backwards in time \cite{chen_neural_2018}, thus incurring a lower memory cost. However, for the systems we investigate here, it is tractable to store the gradients and optimize $\boldsymbol{\lambda}_*$ by directly backpropagating through the dynamics. 

When the total number of time steps is large for a longer protocol duration $\tau$, it is impractical to differentiate through the entire trajectory. During backpropagation, gradient information from each time step is multiplied, via the chain rule, leading to exploding or vanishing gradients \cite{metz_gradients_2021}. To mitigate this, we split our trajectory into $N$ short-time intervals and instead carry out intermediate optimization steps during a cycle. In practice, we see that an interval length, $M = \frac{\tau}{N dt}$, of around 5-10 steps is ideal, with larger $M$ for larger $\tau$ being more expedient. See \url{https://github.com/rotskoff-group/wasserstein-interpolation} for code and the parameters used.




\section{Nanoscale Brownian Engine}
We investigated a system consisting of a one-dimensional overdamped Brownian particle in a heat bath of temperature $T$ confined by a harmonic potential with time-dependent strength $k$,
\begin{equation}
    U_{\rm harmonic}(X(t), \boldsymbol{\lambda}(t)) = \frac{1}{2} k_t X_i^2,
    \label{eq:harmonic}
\end{equation}
with the equation of motion evolving according to an overdamped Langevin equation
\begin{equation}
    dX(t) = -\nabla U_{\rm harmonic}(X(t), \boldsymbol{\lambda}(t)) dt + \sqrt{2\beta^{-1}} dW_t,
    \label{eq:brownian_eom}
\end{equation}
where $W$ is standard Wiener process.

The Brownian engine is driven according to the Stirling cycle, where each stage is of equal duration $\tau/4$. We can fully characterize the cycle given $T_h$, $T_c$, $k_h$ and $k_l$, the maximum and minimum allowed temperature and trap strength respectively. 
Additionally, because the equilibrium distributions are Gaussian with $\mu = 0$ and $\sigma^2 = T/k$, we can exactly specify intermediate distributions along the cycle, because the displacement interpolation that minimizes the Wasserstein distance is Gaussian at each intermediate time in this case~\cite{peyre_computational_2019}. 
We consider the four Gaussian distributions that happen at the beginning of each stage, which have variance $\sigma^2_l$, $\sigma^2_h$, $\sigma^2_h$, and $\sigma^2_l$ respectively, where $\sigma^2_l$ = $T_h/k_h$ and $\sigma^2_h$ = $T_c/k_l$. Along each stage, we compute Wasserstein geodesics $\rho_*(t)$ that interpolate the corresponding endpoints. This geodesic defines a set of intermediate distributions, and we learn a protocol $\boldsymbol{\lambda}_*$ that can successfully drive the engine so that the observed distribution $\hat{\rho}(t)$ is the same as the target distribution $\rho_*(t)$. Because the endpoint distributions are Gaussian, the distributions along the geodesic $\rho_*(t)$ are likewise Gaussian with $\mu_t = 0$ and variance
\begin{equation}
    \sigma_t^2 = \biggl(\frac{t_f - t}{t_f} \sigma_i + \frac{t}{t_f} \sigma_f\biggr)^2,
\end{equation}
where $\sigma^2_i$ is the variance of $\rho_*(0)$ and $\sigma^2_f$ is the variance of $\rho_*(t_f)$, with $0 \leq t \leq t_f$.


We compare the trained $\boldsymbol{\lambda}_*$ to $\boldsymbol{\lambda}_{\rm LR}$, a protocol determined using a linear response approximation \cite{brandner_thermodynamics_2015} that likewise drives the Brownian engine according to the Stirling cycle. To ensure that the protocols were comparable, the maximum and minimum allowed temperatures and harmonic potential strengths were the same for $\boldsymbol{\lambda}_{\rm LR}$ and $\boldsymbol{\lambda}_*$. We used an Euler-Maruyama integration scheme to simulate the dynamics of the Brownian engine and implemented this in Python using the \texttt{PyTorch} and \texttt{torchsde} packages. Because the distributions of this system are 
 Gaussian, we optimized a Mean Squared Error (MSE) loss 
 \begin{equation}
     \mathfrak{L} (t) = (\sigma^2_{\hat{\rho}(t)} - \sigma^2_{\rho_*(t)})^2 + (\mu_{\hat{\rho}(t)} - \mu_{\rho_*(t)})^2.
 \end{equation}
For more complex distributions, it will be expedient to use loss functions, such as the Kullback-Leibler (KL) divergence to estimate the instantaneous derivation of the empirical distribution from the target displacement interpolant. 
Here, $\hat{\rho}$, the empirical distribution of positions simulated under $\boldsymbol{\lambda}_*$, was determined by simulating 1000 independent Brownian engines in parallel.




We compute the heat produced along a cycle as 
\begin{equation}
    Q =  \int -\partial_x U_{\rm harmonic}(X_t, \boldsymbol{\lambda}(t)) \circ dX_{t},
    \label{eq:brownian_heat}
\end{equation}
where $\circ$ is the time-symmetric Stratonovich product.
Similarly, the work done on the system along a trajectory is 
\begin{equation}
    W =  \int \partial_k U_{\rm harmonic}(X_t, \boldsymbol{\lambda}(t))\ dk.
    \label{eq:brownian_work}
\end{equation}
To account for the time varying temperature, we note that we compute the dissipation along the trajectory as  
\begin{equation}
    \beta Q =  \int \frac{-\partial_x U_{\rm harmonic}(X_t, \boldsymbol{\lambda}(t)) \circ dX_{t}}{k_{\rm B} T_t}.
    \label{eq:brownian_dissip}
\end{equation}
With these definitions, we define the efficiency $\eta$ of the engine as
\begin{equation}
    \eta =  \frac{\avg{W}}{\avg{Q_h}},
    \label{eq:efficiencyq}
\end{equation}
where 
\begin{equation}
    Q_h =  \int \gamma_t \bigl(-\partial_x U_{\rm harmonic}(X_t, \boldsymbol{\lambda}(t))\bigr) \circ dX_{t},
    \label{eq:brownian_heat_hot}
\end{equation}
and 
\begin{equation}
    \gamma_t = \frac{T_h(T_l - T_t)}{T_t(T_l - T_h)}
    \label{eq:brownian_gamma}
\end{equation}
as defined in \cite{brandner_thermodynamics_2015}. See Fig~\ref{fig:models} a) and Fig~\ref{fig:w2vL} inset for the dissipation and efficiency of the Brownian engine driven by these $\boldsymbol{\lambda}_*$ and $\boldsymbol{\lambda}_{\rm LR}$. Finally, we compute the squared Wasserstein distance $\mathcal{W}_2^2$ along the isothermal steps of the cycle
\begin{equation}
    \mathcal{W}_2^2 = \sum_0^{\tau/(4 dt)}\bigl((\sigma_{(i + 1)dt} - \sigma_{(i)dt})^2 + ((\mu_{(i + 1)dt} - \mu_{(i)dt})^2\bigr)
\end{equation}
and plot the relationship between the dissipation, $\beta Q$, and $\mathcal{W}_2^2/\tau$ for the isothermal compression step of the Stirling cycle in Fig.~\ref{fig:w2vL}.




\section{Qubit Engine}
We model a system of a single superconducting qubit evolving according to Lindbladian dynamics as introduced in \cite{brandner_thermodynamic_2020}. The Hamiltonian of this system is 
\begin{equation}
    H_v = -\frac{\hbar\Omega}{2}\bigl(\epsilon\sigma_x + \sqrt{V^2 - \epsilon^2}\sigma_z\bigr).
\end{equation}
Here, $\epsilon$ is the coherence parameter, and we fix $\epsilon=0.6$ for the experiments we run. The energy scale is denoted by $\hbar\Omega$, and $\sigma_x$ and $\sigma_z$ correspond to the standard Pauli matrices. We investigate protocols where the level-splitting parameter $V$ of the qubit and the temperature $T$ of the environment can be modulated.
We use the definitions outlined in \cite{brandner_thermodynamic_2020} to simulate and characterize the thermodynamics of the system and present it here again for convenience.


The generalized Bloch equations are 
\begin{equation}
    \begin{pmatrix} \dot{r}_x (t) \\ \dot{r}_y (t) \\ \dot{r}_z (t) \end{pmatrix} = 
    \begin{pmatrix} -k^+_{\boldsymbol{\lambda}_t} & -\Omega V_t & -\theta'_{V_t}\dot{V}_t 
                    \\ \Omega V_t & -k^+_{\boldsymbol{\lambda}_t} & 0
                    \\ \theta'_{V_t}\dot{V}_t & 0 & -2k^+_{\boldsymbol{\lambda}_t}
    \end{pmatrix} 
    \begin{pmatrix} r_x (t)\\ r_y (t)\\ r_z (t)\end{pmatrix}
    + \begin{pmatrix} 0 \\ 0 \\ k^-_{\boldsymbol{\lambda}_t} \end{pmatrix},
\end{equation}  
where 
\begin{equation}
    k^{\pm}_{\boldsymbol{\lambda}_t} \equiv \Gamma\Omega V \frac{1 \pm \exp{\hbar\Omega V/T}}{\exp[\hbar\Omega V/T] - 1}
\end{equation}
and
\begin{equation}
    \theta'_V \equiv \frac{\epsilon}{V\sqrt{V^2 - \epsilon^2}}.
\end{equation}
We compute the heat produced by the engine as
\begin{equation}
    Q = \frac{1}{2} \int_0^{\tau} \bigl(4 r_t  \arctanh (2r_t) + \ln[1/4 - r_t^2]\bigr)\ \dot{T}_t \ dt
\end{equation}
and the work done by the system as 
\begin{equation}
    W = -\hbar\Omega\int_0^{\tau} \bigl(-V_t\theta'_{V_t}r_t^x + r_t^z\bigr)\dot{V}_t \ dt,
\end{equation}
where the length of the Bloch vector 
$r_t \equiv \sqrt{(r_t^x)^2 + (r_t^y)^2 + (r_t^z)^2}$.
With these definitions, we compute the dissipation as
\begin{equation}
    \beta Q = \frac{1}{2} \int_0^{\tau} \frac{1}{T_t} \bigl(4 r_t  \arctanh (2r_t) + \ln[1/4 - r_t^2]\bigr)\ \dot{T}_t \ dt.
\end{equation}
Finally, as in \cite{brandner_thermodynamic_2020} we define the efficiency, 
\begin{equation}
    \eta =  \frac{W}{Q}.
    \label{eq:efficiency}
\end{equation}


From the Bloch vector, we can compute the density matrix $\rho_t$ as 
\begin{equation}
    \rho_t \equiv \mathbb{I}/2 + r_t^x \Pi^x_{V_t} + r_t^y \Pi^y_{V_t}  + r_t^z \Pi^z_{V_t},
\end{equation}
where 
\begin{align*} 
\Pi^x_{V_t} &= \lvert E_V^+ \rangle\langle E_V^- \rvert +  \lvert E_V^- \rangle\langle E_V^+ \rvert, \\ 
\Pi^y_{V_t} &= i \lvert E_V^-\rangle\langle E_V^+ \rvert - i\lvert E_V^+ \rangle\langle E_V^- \rvert, \\ 
\Pi^z_{V_t} &= \lvert E_V^+ \rangle\langle E_V^+ \rvert -  \lvert E_V^- \rangle\langle E_V^- \rvert.
\end{align*}
Lastly, we define 
\begin{align*}
    E_V^+ &=  \begin{pmatrix} \phantom{-} \sin[\theta_V/2] \\ -\cos[\theta_V/2] \end{pmatrix}, \\
    E_V^- &=  \begin{pmatrix} \cos[\theta_V/2] \\ \sin[\theta_V/2] \end{pmatrix}, \\
\end{align*}
with $\sin[\theta_V/2] = \epsilon/\sqrt{2V^2 + 2V\sqrt{V^2 - \epsilon^2}}$.


We first consider a base protocol $\boldsymbol{\lambda}(\bar{t}) = [T_{\bar{t}}, V_{\bar{t}}]$, where $\bar{t} = t/\tau$, and $\tau$ is the duration of the protocol. As in \cite{brandner_thermodynamic_2020}, we define $\boldsymbol{\lambda}(\bar{t})$ to be sinusoidal,
\begin{equation}
    \boldsymbol{\lambda}(\bar{t}) = \hbar\Omega[(1 + \sin^2(\pi\Omega\bar{t}), 1 + \sin^2(\pi\Omega \bar{t} + \pi/4)].
\end{equation}
To compute the optimal speed function $\phi_*$ that minimizes \eqref{eq:stot}, we first simulate the qubit engine under $\boldsymbol{\lambda}$ in the slow-driving regime with $\tau_{\rm slow} = 50$. Using the periodic steady-state $\rho_t$, we compute an ``instantaneous length'' of the path
\begin{equation}
    dL_{t_i} = \sqrt{\bigl(\rho_{t_i+dt}(\lvert 1 \rangle) - \rho_{t_i}(\lvert 1 \rangle)\bigr)^2 + \bigl(\rho_{t_i+dt}(\lvert 0 \rangle)  -\rho_{t_i}(\lvert 0 \rangle)\bigr)^2},
\end{equation}
where $0\leq i\leq\tau_{\rm slow}/dt$.
We then compute the discretized time-derivative of the speed function
\begin{equation}
    \dot{\phi}_{*_i} = \biggl(\frac{dt}{\tau_{\rm slow}}\biggr)\frac{L}{dL_{t_i}},
    \label{eq:speedot}
\end{equation}
where $L = \sum_{i} dL_{t_i}$. 
This leads to an effective time $\bar{t}_i = \sum_{k = 0}^{i}\dot{\phi}_{*_k}$, where $\dot{\phi}_* (\bar{t}_i) = \dot{\phi}_{*_i}$  and $0 \leq \bar{t}_i \leq 1$. Finally, the continuous function $\dot{\phi}_*(\bar{t})$ can be recovered via linear interpolation and $\phi_*(\bar{t})$ by integration. 
Intuitively, when the instantaneous path length $dL_{t_i}$ is larger than the average path length $\bigl(dt/\tau_{\rm slow}\bigr)\mathcal{L}$ for a time step, the protocol moves slower than a constant speed protocol in order to minimize the dissipative cost; when the instantaneous path length is smaller than the average path length the protocol moves faster than a constant speed protocol.


\begin{figure}[ht]
    \centering
    \includegraphics[width=0.45\linewidth]{qubit_efficiency.pdf}
    \includegraphics[width=0.45\linewidth]{qubit_densities.pdf}
    \caption{The efficiency (a) and steady-state densities at $\tau = 5$ (b) of the qubit engine for various speed functions for a prespecified sinusoidal protocol. The efficiency of the engine is higher under a speed function determined via a linear response approach (green triangle) compared to the speed function determined by minimizing \eqref{eq:stot} (dark blue circle). However, there is no guarantee of the intermediate states that the system visits when learning speed functions for a prespecified protocol.}
    \label{fig:qubit_dissip_eff_speed}
\end{figure}


We computed $\phi_{\rm LR} (t)$ using the linear response approximation presented in \cite{brandner_thermodynamic_2020}. In practice, the procedure for computing $\phi_{\rm LR} (t)$ is the same as computing $\phi_{*} (t)$, with a different definition of the length metric $\mathcal{L}$ (cf. Eq. S45 in \cite{brandner_thermodynamic_2020}). Fig.~\ref{fig:models} b) shows a comparison of the dissipation profiles of the qubit engine driven by $\boldsymbol{\lambda}(\phi_{*} (t))$, $\boldsymbol{\lambda}(\phi_{\rm LR} (t))$, and $\boldsymbol{\lambda}(t)$. We also plot the efficiency $\eta$ of the qubit engine driven under these protocols and observe that the efficiency remains the highest under $\boldsymbol{\lambda}(\phi_{\rm LR} (t))$ (Fig.~\ref{fig:qubit_dissip_eff_speed}). Optimizing \eqref{eq:stot} serves to minimize the dissipation, and there is no guarantee of a higher efficiency. Finally, we plot the steady-state density of the qubit engine driven by the three protocols, as shown in Fig.~\ref{fig:qubit_dissip_eff_speed}.


As evident in Fig~\ref{fig:qubit_dissip_eff_speed} b), there is no guarantee of the intermediate states that a system will visit when determining a speed function for a prespecified protocol $\boldsymbol{\lambda}$. Next, we instead consider learning a $\boldsymbol{\lambda}_*(t)$ using Algorithm~\ref{alg:autodiff} that cyclically drives the system through a set of prespecified intermediate distribution. We specified four intermediate distributions as the steady-state distributions the qubit engine driven under $\boldsymbol{\lambda}(t)$ at $t = 0, \tau/4, \tau/2,$ and $3\tau/4$ in the quasistatic limit. We then computed four intermediate Wasserstein geodesics that interpolated these endpoints, where
\begin{equation}
    \rho (\lvert 1 \rangle, t) = \biggl(\frac{t_f - t}{t_f} \rho (\lvert 1 \rangle, t_f) + \frac{t}{t_f} \rho (\lvert 1 \rangle, 0)\biggr),
\end{equation}
$0\leq t\leq t_f$ and $\rho (\lvert 0 \rangle, t) = 1 -  \rho (\lvert 1 \rangle, t)$. 
Finally, we carried out the automatic differentiation optimization presented in Algorithm \ref{alg:autodiff} using a Mean Squared Error (MSE) objective
\begin{equation}
    \mathfrak{L}(t)= \frac{1}{2} \bigl(\rho_* (\lvert 0 \rangle, t) - \hat{\rho} (\lvert 0 \rangle,  t)\bigr)^2 + \bigl(\rho_* (\lvert 1 \rangle, t) - \hat{\rho} (\lvert 1 \rangle, t)\bigr)^2.
\end{equation}

For each of the four intervals of the protocol, $[0, \tau/4]$, $[\tau/4, \tau/2]$, $[\tau/2, 3\tau/4]$, and $[3\tau/4, \tau]$, we constrained the maximum and minimum $T$ and $V$ for $\boldsymbol{\lambda}_*$ to be the same as the maximum and minimum $T$ and $V$ for the base sinusoidal protocol $\boldsymbol{\lambda}$. Finally, we define $\boldsymbol{\lambda}_{\rm linear}$, a protocol which linearly interpolates $\boldsymbol{\lambda}(0)$, $\boldsymbol{\lambda}(\tau/4)$, $\boldsymbol{\lambda}(\tau/2)$, $\boldsymbol{\lambda}(3\tau/4)$, and $\boldsymbol{\lambda}(\tau)$. We compare the dissipation and efficiencies of the qubit engine driven by $\boldsymbol{\lambda}_*$ and $\boldsymbol{\lambda}_{\rm linear}$ (see Fig~\ref{fig:qubit_int}). We observe that under $\boldsymbol{\lambda}_*$ the qubit engine is less dissipative, but also slightly less efficient. Importantly, the former observation is consistent with~\eqref{eq:stot}; minimizing ~\eqref{eq:stot} does not necessarily guarantee higher efficiency. 

We compare these two protocols because we observe that the protocols generally have similar intermediate states for $\tau > 1$ (see Fig.~\ref{fig:qubit_int}). For the fastest driving, $\tau = 1$, both protocols deviate significantly from the target distribution as it is infeasible to learn a protocol $\boldsymbol{\lambda}_*$ that can achieve the target distributions on this timescale. However, we see that under $\boldsymbol{\lambda}_{\rm linear}$, the steady-state distribution of the qubit engine deviates slightly from the target Wasserstein geodesic especially when $t > \tau/2$. This ultimately results in higher dissipation when driven by $\boldsymbol{\lambda}_{\rm linear}$ compared to $\boldsymbol{\lambda}_{*}$.



\begin{figure}[ht]
    \centering
    \includegraphics[width=0.45\linewidth]{qubit_dissipation_intermediate_states.pdf}
    \includegraphics[width=0.45\linewidth]{qubit_efficiency_intermediate.pdf}
    \includegraphics[width=0.45\linewidth]{qubit_densities_intermediate.pdf}
    \caption{The dissipation (a), efficiency (b) and steady-state densities at $\tau=5$ (c) for protocols optimized to traverse through target intermediate distributions compared to a linear protocol. The dissipation is lower under the optimal protocol that minimizes \eqref{eq:stot} (dark blue circle) compared to a linear protocol (red triangle). For fastest-driving ($\tau = 1$), both protocols deviate significantly from intermediate target distributions, resulting in a high dissipative cost. For $\tau > 1$ optimal protocol achieves target intermediate distributions resulting in a lower dissipative cost. Still, the efficiency is higher for the  linear protocol compared to the optimal protocol. Finally, the optimal protocol closely realizes the target steady-state distribution, while the linear protocol slightly deviates from the target intermediate distributions for $t > \tau/2$ resulting in a higher dissipative cost.}
    \label{fig:qubit_int}
\end{figure}




