# Supplement:
Thermodynamic Geometry of Nonequilibrium Fluctuations in Cyclically Driven
Transport

Zi Wang  Center for Phononics and Thermal Energy Science, China-EU Joint Lab
on Nanophononics, Shanghai Key Laboratory of Special Artificial Microstructure
Materials and Technology, School of Physics Science and Engineering, Tongji
University, Shanghai 200092, China  Jie Ren  Corresponding Email:
<EMAIL_ADDRESS>Center for Phononics and Thermal Energy Science, China-EU
Joint Lab on Nanophononics, Shanghai Key Laboratory of Special Artificial
Microstructure Materials and Technology, School of Physics Science and
Engineering, Tongji University, Shanghai 200092, China

###### Abstract

In Sec. I, we provide the details of the full counting statistics (FCS) in the
Markovian dynamics. In Sec. II, we present the derivation of our geometric
separation of FCS in a periodically driven thermal device and show its
manifestation as a metric term. In Sec. III, the formulae for average currents
are derived based on our geometric theory on FCS. In Sec. IV, we illustrate
our geometric non-adiabatic control principle. In Sec. V, we derive our
geometric thermodynamic uncertainty relation (TUR) and its implication on the
fluctuation of entropy production. Finally, the details for the two models
exemplified in the main text are provided in Sec. VI.

##  I General Markovian Dynamics

In this section, we discuss the basic theory for the generating function of a
fluctuating current.

We consider a periodically driven discrete Markovian system coupled to several
reservoirs. Its probability distribution $p_{i}$ is labeled by $i$ ($1\leq
i\leq N$), with $N$ being the overall number of states. The rate along the
transition $j\to i$ induced by the $\nu$-th reservoir is $k_{ij}^{\nu}$
($i\neq j$), whose time-dependence is kept implicit for clarity. As demanded
by the thermodynamic consistency, we assume the local detailed balance
condition $k_{ij}^{\nu}/k_{ji}^{\nu}=e^{\beta_{\nu}(E_{j}-E_{i})}$, with
$\beta_{\nu}:=1/T_{\nu}$ being the inverse temperature of the $\nu$-th
reservoir (Boltzmann constant is set to $1$) and $E_{i}$ the energy level of
the system. The master equation of $p_{i}$ is
$\partial_{t}\ket{p(t)}=\hat{L}\ket{p(t)}$, with
$L_{ij}:=\sum_{\nu}(k_{ij}^{\nu}-\delta_{ij}\sum_{j\neq i}k_{ji}^{\nu})$
preserving probability by $\sum_{i}L_{ij}=0$.

The transition $k_{ij}^{\nu}$ is associated with an accumulation of exchange
current $\Delta Q=d_{ij}^{\nu}$ ($d_{ij}^{\nu}=-d_{ji}^{\nu}$). The
accumulated current $Q$ can be particle number, energy, entropy, work, and so
on, which in the special case of $d_{ij}^{\nu}=\ln(k_{ij}^{\nu}/k_{ji}^{\nu})$
corresponds to the reservoir entropy production $\Sigma$. To consider the
fluctuation of this current, we take advantage of the full counting statistics
formalism. We twist the operator $\hat{L}$ to $\hat{L}_{\chi}$, with
$L_{\chi,ij}=\sum_{\nu}k_{ij}^{\nu}e^{\chi d_{ij}^{\nu}}$ for $i\neq j$ and
$L_{\chi,ii}=L_{ii}$. Here, $\chi$ is the counting field of $Q$. Taking the
derivative of the cumulant generating function (CGF)

${\mathcal{G}}=\ln{\mathcal{Z}}=\ln\braket{1}{p_{\chi}(t)}$ (S1)

with respect to $\chi$ yields the $n$-th cumulant of the accumulated current
$Q$ during $[0,t]$

$\braket{Q^{n}}_{c}=\frac{\partial^{n}{\mathcal{G}}}{\partial\chi^{n}}|_{\chi=0}.$
(S2)

Here, $\ket{p_{\chi}(t)}$ is generated by the equation of motion
$\partial_{t}\ket{p_{\chi}(t)}=\hat{L}_{\chi}\ket{p_{\chi}(t)}$ and $\bra{1}$
is a vector with all elements being $1$. We note that the characteristic
function ${\mathcal{Z}}$ satisfies ${\mathcal{Z}}|_{\chi=0}=1$ due to the
probability normalization condition $\braket{1}{p_{\chi}}|_{\chi=0}=1$. As a
result, the CGF has the property ${\mathcal{G}}|_{\chi=0}=0$.

This is equivalent to the path integral formalism. By assigning the weight to
a trajectory $\omega=x(t)$ ($0\leq t\leq\tau$)

$P[\omega]=p(x_{0})e^{-\int_{0}^{\tau}dt\sum_{i\neq
j;\nu}[\delta_{x(t),i}k_{ji}^{\nu}(t)-\dot{m}_{ji}^{\nu}\ln
k_{ji}^{\nu}(t)]},$ (S3)

with $p(x_{0})$ being the initial distribution at $t=0$ and $m_{ij}^{\nu}$
counting the accumulated number of transition $k_{ij}^{\nu}$, the CGF is also
expressed as

${\mathcal{G}}=\ln\int d\omega p(x_{0})e^{-\int_{0}^{\tau}dt\sum_{i\neq
j;\nu}[\delta_{x(t),i}k_{ji}^{\nu}(t)-\dot{m}_{ji}^{\nu}(\chi
d_{ji}^{\nu}(t)+\ln k_{ji}^{\nu}(t))]}.$ (S4)

The overdamped Brownian dynamics can be taken as a specific limit case of this
discrete transition dynamics [1].

##  II Dynamic and Geometric Components of the CGF

Here, we derive the general geometric components of the current fluctuation
and its metric component.

We start from the equation of motion
$\partial_{t}\ket{p_{\chi}(t)}=\hat{L}_{\chi}({\bm{\Lambda}}(t))\ket{p_{\chi}(t)}$,
where ${\bm{\Lambda}}(t)$ is the periodically driven parameter vector with
period $\tau_{p}$. Suppose that $\hat{L}_{\chi}$ can be diagonalized
$\hat{L}_{\chi}=\sum_{n}E_{n}\ket{r_{n}}\bra{l_{n}}$, with $E_{n}$,
$\ket{r_{n}}$ and $\bra{l_{n}}$ being respectively the eigenvalue, right- and
left-eigenvectors satisfying $\hat{L}_{\chi}\ket{r_{n}}=E_{n}\ket{r_{n}}$ and
$\bra{l_{n}}\hat{L}_{\chi}=E_{n}\bra{l_{n}}$ and the orthonormal condition
$\braket{l_{m}}{r_{n}}=\delta_{m,n}$. The eigenvectors $\ket{r_{n}}$ and
$\bra{l_{n}}$ and the eigenvalues $E_{n}$ are all functions of the counting
parameter $\chi$. $n=0$ corresponds to the steady state. We also note that by
taking $\chi=0$, the twisted operator $\hat{L}_{\chi}$ reduces to $\hat{L}$,
$\ket{r_{0}}$ to the steady state distribution $\ket{\pi}$ and $\bra{l_{0}}$
to the one-vector $\bra{1}$. After a large number of periods, the system
enters its cyclic state. According to the Floquet theorem, the system state is
of the form

$\ket{p_{\chi}(t)}=e^{{\mathcal{G}}(t)}\ket{\phi(t)},$ (S5)

where $\ket{\phi(t+\tau_{p})}=\ket{\phi(t)}$ is a cyclic state and
$\ket{p_{\chi}(t)}$ only accumulates an CGF phase ${\mathcal{G}}(\tau_{p})$
during one driving period. Inserting this ansatz into the equation of motion,
the dynamics of $\ket{\phi(t)}$ is governed by

$(\hat{L}_{\chi}-\partial_{t})\ket{\phi(t)}=(\partial_{t}{\mathcal{G}}(t))\ket{\phi(t)}.$
(S6)

The operator $\partial_{t}$ can be seen as a perturbation to the steady state.
In the absence of $\partial_{t}$, Eq. (S6) is solved by the state
$\ket{\phi(t)}=\ket{r_{0}}$ and the equation for CGF

$\partial_{t}{\mathcal{G}}_{\rm dyn}(t)=E_{0}(t),$ (S7)

which corresponds to the dynamic-phase-like contribution

${\mathcal{G}}_{\rm dyn}(\tau_{\rm p})=\int_{0}^{\tau_{\rm p}}dtE_{0}(t).$
(S8)

Denoting the deviation of ${\mathcal{G}}(t)$ from ${\mathcal{G}}_{\rm dyn}(t)$
as ${\mathcal{G}}_{\rm geo}(t)={\mathcal{G}}(t)-{\mathcal{G}}_{\rm dyn}(t)$,
which is later shown to be a general geometric phase like contribution
applicable to arbitrary driving frequency, we can write Eq. (S6) as

$[\hat{L}_{\chi}-E_{0}(t)]\ket{\phi(t)}=[\partial_{t}+\partial_{t}{\mathcal{G}}_{\rm
geo}(t)]\ket{\phi(t)}.$ (S9)

Taking the inner product of both sides with $\bra{l_{0}}$, the left hand side
is zero. Thus, the right hand side is orthogonal to $\ket{r_{0}}$ and we
derive the equation for ${\mathcal{G}}_{\rm geo}$ as
$(\partial_{t}{\mathcal{G}}_{\rm
geo})\braket{l_{0}}{\phi}=-\braket{l_{0}}{\partial_{t}\phi}$. This leaves a
gauge degree of freedom for $\ket{\phi}$. If we make the gauge transform
$\ket{\phi}\to e^{\alpha(t)}\ket{\phi}$, $\partial_{t}{\mathcal{G}}_{\rm geo}$
is changed to $\partial_{t}{\mathcal{G}}_{\rm geo}+\partial_{t}\alpha(t)$,
ensuring the gauge-invariance of ${\mathcal{G}}_{\rm geo}(\tau_{\rm p})$ since
$\int_{0}^{\tau_{\rm p}}dt\partial_{t}\alpha=0$.

We can fix the gauge of $\ket{\phi(t)}$ by requiring its form

$\ket{\phi}=\ket{r_{0}}+\ket{\phi_{\bot}},$ (S10)

where $\braket{l_{0}}{\phi_{\bot}}=0$. Substituting this expression into Eq.
(S9), we obtain

$(\hat{L}_{\chi}-E_{0})\ket{\phi_{\bot}}=\ket{\partial_{t}r_{0}}+\ket{\partial_{t}\phi_{\bot}}+(\partial_{t}{\mathcal{G}}_{\rm
geo})(\ket{r_{0}}+\ket{\phi_{\bot}}),$ (S11)

the inner product of which and $\bra{l_{0}}$ provides an equation for
${\mathcal{G}}_{\rm geo}$:

$\begin{split}\partial_{t}{\mathcal{G}}_{\rm
geo}&=-\braket{l_{0}}{\partial_{t}r_{0}}-\braket{l_{0}}{\partial_{t}\phi_{\bot}}:=\partial_{t}{\mathcal{G}}_{\rm
curv}+\partial_{t}{\mathcal{G}}_{\rm metr},\end{split}$ (S12)

where we identify the curvature term and the metric term. They describe
respectively the adiabatic and nonadiabatic effect of the driving. Through a
whole period, the curvature term is

${\mathcal{G}}_{\rm curv}(\tau_{\rm
p})=\oint_{\partial\Omega}d\Lambda_{\mu}A_{\mu}:=-\oint_{\partial\Omega}d\Lambda_{\mu}\braket{l_{0}}{\partial_{\mu}r_{0}},$
(S13)

where we define the geometric connection
$A_{\mu}=-\braket{l_{0}}{\partial_{\mu}r_{0}}$. Using the Stokes formula, this
is equivalent to

${\mathcal{G}}_{\rm curv}(\tau_{\rm p})=\int_{\Omega}dS_{\mu\nu}F_{\mu\nu},$
(S14)

in terms of the geometric curvature
$F_{\mu\nu}:=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}=-\braket{\partial_{\mu}l_{0}}{\partial_{\nu}r_{0}}+\braket{\partial_{\nu}l_{0}}{\partial_{\mu}r_{0}}$.

Now we focus on the metric term ${\mathcal{G}}_{\rm metr}$. We first note

$\begin{split}{\mathcal{G}}_{\rm metr}(\tau_{\rm p})&=-\int_{0}^{\tau_{\rm
p}}dt\braket{l_{0}}{\partial_{t}\phi_{\bot}}\\\ &=\int_{0}^{\tau_{\rm
p}}dt\braket{\partial_{t}l_{0}}{\phi_{\bot}}.\end{split}$ (S15)

Since Eq. (S11) can be formally solved as

$\begin{split}\ket{\phi_{\bot}}&=(\hat{L}_{\chi}-E_{0})^{+}\ket{\partial_{t}\phi}+(\partial_{t}{\mathcal{G}}_{\rm
geo})(\hat{L}_{\chi}-E_{0})^{+}\ket{\phi_{\bot}}\\\
&=(\hat{L}_{\chi}-E_{0})^{+}\ket{\partial_{t}\phi}-(\hat{L}_{\chi}-E_{0})^{+}\ket{\phi}\braket{l_{0}}{\partial_{t}\phi},\end{split}$
(S16)

with the pseudo-inverse being $(\hat{L}_{\chi}-E_{0})^{+}:=\sum_{n\neq
0}\frac{1}{E_{n}-E_{0}}\ket{r_{n}}\bra{l_{n}}$, we derive the metric
expression

$\begin{split}{\mathcal{G}}_{\rm metr}(\tau_{\rm p})&=\int_{0}^{\tau_{\rm
p}}dtg_{tt}:=\int_{0}^{\tau_{\rm
p}}dt\braket{\partial_{t}l_{0}}{\hat{G}}{\partial_{t}\phi}.\end{split}$ (S17)

Here, we define the metric tensor as

$\hat{G}=(\hat{L}_{\chi}-E_{0})^{+}-(\hat{L}_{\chi}-E_{0})^{+}\ket{\phi}\bra{l_{0}}.$
(S18)

This metric structure also endows a pseudo-Riemannian manifold in the driven
parameter space. In terms of this, ${\mathcal{G}}_{\rm metr}$ can be written
as

$\displaystyle{\mathcal{G}}_{\rm metr}$ $\displaystyle=$
$\displaystyle\int_{0}^{\tau_{\rm
p}}g_{\mu\nu}\dot{\Lambda}_{\mu}\dot{\Lambda}_{\nu}dt,$ (S19)
$\displaystyle\text{with}\;\;g_{\mu\nu}$ $\displaystyle:=$
$\displaystyle\frac{1}{2}\left[\braket{\partial_{\mu}l_{0}}{\hat{G}}{\partial_{\nu}\phi}+\braket{\partial_{\nu}l_{0}}{\hat{G}}{\partial_{\mu}\phi}\right].$

In contrast to the Riemannian manifold, the metric $[g_{\mu\nu}]$ is not
necessarily positive-definite. This sacrifice makes our formalism applicable
to the whole fluctuation properties of an arbitrary current.

When the near-adiabatic regime is concerned, the Floquet state $\ket{\phi}$
can be approximated, to the leading order, as $\ket{r_{0}}$. This provides the
nonadiabatic metric structure

$\begin{split}\mathfrak{g}_{\mu\nu}&=\frac{1}{2}\left[\braket{\partial_{\mu}l_{0}}{\hat{G}}{\partial_{\nu}r_{0}}+\braket{\partial_{\nu}l_{0}}{\hat{G}}{\partial_{\mu}r_{0}}\right]\\\
&=\sum_{n\neq
0}\frac{\braket{\partial_{\mu}l_{0}}{r_{n}}\braket{l_{n}}{\partial_{\nu}r_{0}}+\braket{\partial_{\nu}l_{0}}{r_{n}}\braket{l_{n}}{\partial_{\mu}r_{0}}}{2(E_{n}-E_{0})},\end{split}$
(S20)

which describes the leading finite-time effect in the near-adiabatic regime.
Away from this regime, the exact $\ket{\phi}$ can be obtained by using the
Dyson-series-like self-consistent equation Eq. (S16) with iteration.

Therefore, as a summary the CGF has the decomposition:

${\mathcal{G}}={\mathcal{G}}_{\rm dyn}+{\mathcal{G}}_{\rm
geo}={\mathcal{G}}_{\rm dyn}+{\mathcal{G}}_{\rm curv}+{\mathcal{G}}_{\rm
metr}.$ (S21)

We note that the above formalism is applicable to arbitrary time antisymmetric
current in general nonequilibrium conditions. This provides the basis for
analyzing the fluctuation properties of a finite-time cyclically driven
thermal device.

##  III Average geometric currents

Here, we derive the average current based on our FCS theory.

In this section, $\ket{p}$ and $\hat{L}$ without the subscript $\chi$ means
$\ket{p_{\chi}}|_{\chi=0}$ and $\hat{L}_{\chi}|_{\chi=0}$.

First, we consider the total accumulated current and show that the mean value
$\braket{Q}=\partial_{\chi}{\mathcal{G}}(\tau_{p})|_{\chi=0}$ is also
expressed as

$\braket{Q}=\int_{0}^{\tau_{p}}dt\braket{1}{\hat{J}(t)}{p(t)},$ (S22)

where the current density operator is given by
$\hat{J}=\partial_{\chi}\hat{L}_{\chi}|_{\chi=0}$. Eq. (S22) is important in
studying both the average current in stochastic systems and the current in
deterministic systems. Defining the propagator for the twisted dynamics
$\partial_{t}\ket{p_{\chi}(t)}=\hat{L}_{\chi}\ket{p_{\chi}(t)}$ as

$\hat{U}_{\chi}(t_{2},t_{1})={\mathcal{T}}e^{\int_{t_{1}}^{t_{2}}dt\hat{L}_{\chi}(t)},$
(S23)

in terms of the time-ordering operator ${\mathcal{T}}$, we note that the
Floquet state Eq. (S5) is the eigenstate of $\hat{U}_{\chi}(\tau_{p},0)$

$\hat{U}_{\chi}(\tau_{p},0)\ket{p_{\chi}(0)}=e^{{\mathcal{G}}(\tau_{p})}\ket{p_{\chi}(0)}.$
(S24)

We differentiate both side with respect to $\chi$, set $\chi=0$ and take the
inner product with $\bra{1}$. The left hand side is

$\begin{split}&\braket{1}{[\partial_{\chi}\hat{U}_{\chi}(\tau_{p},0)]}{p_{\chi}(0)}|_{\chi=0}+\braket{1}{\hat{U}_{\chi}(\tau_{p},0)}{\partial_{\chi}p_{\chi}(0)}|_{\chi=0}\\\
&=\int_{0}^{\tau_{p}}dt\braket{1}{\hat{J}}{p(t)}+\braket{1}{\partial_{\chi}p_{\chi}(0)}|_{\chi=0},\end{split}$
(S25)

using
$\partial_{\chi}\hat{U}_{\chi}(\tau_{p},0)=\int_{0}^{\tau_{p}}dt\hat{U}_{\chi}(\tau_{p},t)[\partial_{\chi}\hat{L}_{\chi}(t)]\hat{U}_{\chi}(t,0)$
and $\bra{1}\hat{U}_{\chi}|_{\chi=0}=\bra{1}$. The right hand side is

$\begin{split}&[\partial_{\chi}{\mathcal{G}}(\tau_{p})]e^{{\mathcal{G}}(\tau_{p})}\braket{1}{p_{\chi}(0)}|_{\chi=0}+e^{{\mathcal{G}}(\tau_{p})}\braket{1}{\partial_{\chi}p_{\chi}(0)}|_{\chi=0}\\\
&=\braket{Q}+\braket{1}{\partial_{\chi}p_{\chi}(0)}|_{\chi=0},\end{split}$
(S26)

using ${\mathcal{G}}(\tau_{p})|_{\chi=0}=0$, the definition
$\braket{Q}:=\partial_{\chi}{\mathcal{G}}(\tau_{p})|_{\chi=0}$ and the
normalization condition
$\braket{1}{p_{\chi}(0)}|_{\chi=0}=\braket{1}{p(0)}=1$. Combining these
results, we arrive at the identity for average total accumulated current Eq.
(S22).

Second, we consider the dynamic current component. To prove

$\braket{Q_{\rm dyn}}:=\partial_{\chi}{\mathcal{G}}_{\rm
dyn}(\tau_{p})|_{\chi=0}=\int_{0}^{\tau_{p}}dt\braket{1}{\hat{J}}{\pi(t)},$
(S27)

with $\ket{\pi}=\ket{r_{0}}|_{\chi=0}$ being the instantaneous steady state,
we start from the eigenvalue equation

$\hat{L}_{\chi}\ket{r_{0}}=E_{0}\ket{r_{0}}.$ (S28)

Calculating the derivative at $\chi=0$ and taking inner product with
$\bra{1}$, the left hand side is
$\braket{1}{\hat{J}}{\pi}+\braket{1}{\hat{L}_{\chi}}{\partial_{\chi}r_{0}}|_{\chi=0}=\braket{1}{\hat{J}}{\pi}$
using the property $\bra{1}\hat{L}=0$, while the right hand side is
$\partial_{\chi}E_{0}|_{\chi=0}$ since $\ket{\pi}$ is normalized and
$E_{0}|_{\chi=0}=0$. Therefore, we have the result
$\partial_{\chi}E_{0}|_{\chi=0}=\braket{1}{\hat{J}(t)}{\pi(t)}$ and
consequently Eq. (S27).

Thirdly, we consider the general geometric current. For the cyclic state
$\ket{\phi(t)}$, we note that the corresponding average current density
$-\partial_{\chi}\braket{l_{0}}{\partial_{t}\phi}|_{\chi=0}$ is identical to
$-\braket{\partial_{\chi}l_{0}}{\partial_{t}\phi}|_{\chi=0}$ up to a complete
differential term
$-\braket{l_{0}}{\partial_{\chi}\partial_{t}\phi}|_{\chi=0}=-\partial_{t}\braket{1}{\partial_{\chi}\phi}|_{\chi=0}$,
for any cyclic $\ket{\phi}$. This complete differential term makes no
contribution when integrating over a whole period. By further showing the
identity

$\bra{\partial_{\chi}l_{0}}\hat{L}|_{\chi=0}+\bra{1}\hat{J}=\bra{1}(\partial_{\chi}E_{0})|_{\chi=0},$
(S29)

and therefore

$\begin{split}\bra{\partial_{\chi}l_{0}}|_{\chi=0}&=\braket{\partial_{\chi}l_{0}}{r_{0}}|_{\chi=0}\bra{1}-\bra{1}\hat{J}\hat{L}^{+},\end{split}$
(S30)

where $\hat{L}^{+}$ is the pseudo-inverse of $\hat{L}$, we prove that our
formulation in terms of CGF provides the expression for average accumulated
geometric current

$\begin{split}\braket{Q_{\rm geo}}&:=\partial_{\chi}{\mathcal{G}}_{\rm
geo}|_{\chi=0}=-\int_{0}^{\tau_{\rm
p}}dt\partial_{\chi}\braket{l_{0}}{\partial_{t}\phi}|_{\chi=0}\\\
&=-\int_{0}^{\tau_{\rm
p}}dt\braket{\partial_{\chi}l_{0}}{\partial_{t}\phi}|_{\chi=0}\\\
&=-\int_{0}^{\tau_{\rm
p}}dt\left[\braket{\partial_{\chi}l_{0}}{r_{0}}\braket{1}{\partial_{t}\phi}-\braket{1}{\hat{J}\hat{L}^{+}}{\partial_{t}\phi}\right]|_{\chi=0}\\\
&=\int_{0}^{\tau_{\rm
p}}dt\braket{1}{\hat{J}\hat{L}^{+}}{\partial_{t}p(t)}.\end{split}$ (S31)

The fourth equality is due to the constant normalization condition
$\partial_{t}\braket{1}{\phi}|_{\chi=0}=0$. This expression of $\braket{Q_{\rm
geo}}$ simplifies to the adiabatic geometric current

$\begin{split}\braket{Q_{\rm
curv}}&=\oint_{\partial\Omega}d\Lambda_{\mu}\braket{1}{\hat{J}\hat{L}^{+}}{\partial_{\mu}\pi(t)}\\\
&=\int_{\Omega}d\Lambda_{\mu}d\Lambda_{\nu}[\braket{1}{\partial_{\mu}(\hat{J}\hat{L}^{+})}{\partial_{\nu}\pi}-\braket{1}{\partial_{\nu}(\hat{J}\hat{L}^{+})}{\partial_{\mu}\pi}],\end{split}$
(S32)

and the nonadiabatic geometric current

$\begin{split}\braket{Q_{\rm metr}}&=\int_{0}^{\tau_{\rm
p}}dt\braket{1}{\hat{J}\hat{L}^{+}}{\partial_{t}\phi_{\bot}}|_{\chi=0}\\\
&=\int_{0}^{\tau_{\rm
p}}dt\braket{1}{\hat{J}\hat{L}^{+}\partial_{t}[(\hat{L}_{\chi}-E_{0})^{+}}{\partial_{t}\phi}]_{\chi=0}\\\
&-\int_{0}^{\tau_{\rm
p}}dt\braket{1}{\hat{J}\hat{L}^{+}\partial_{t}[(\hat{L}_{\chi}-E_{0})^{+}}{\phi}\braket{l_{0}}{\partial_{t}\phi}]_{\chi=0}\\\
&=\int_{0}^{\tau_{\rm
p}}dt\braket{1}{\hat{J}\hat{L}^{+}\partial_{t}[\hat{L}^{+}}{\partial_{t}p}]\\\
&=\int_{0}^{\tau_{\rm
p}}dt\dot{\Lambda}_{\mu}\dot{\Lambda}_{\nu}\frac{1}{2}\left[\bra{1}\hat{J}\hat{L}^{+}\partial_{\mu}(\hat{L}^{+}\ket{\partial_{\nu}p})+(\mu\leftrightarrow\nu)\right]\\\
&:=\int_{0}^{\tau_{\rm
p}}dt\dot{\Lambda}_{\mu}\dot{\Lambda}_{\nu}g_{\mu\nu}^{Q}.\end{split}$ (S33)

Here, $g_{\mu\nu}^{Q}$ is the metric related to the average current $Q$.

Correspondingly, the near-adiabatic metric for average $Q$ is obtained by
replacing $\ket{p}$ in Eq. (S33) by $\ket{\pi}$:

$\mathfrak{g}_{\mu\nu}^{Q}=\frac{1}{2}\left[\bra{1}\hat{J}\hat{L}^{+}\partial_{\mu}(\hat{L}^{+}\ket{\partial_{\nu}\pi})+(\mu\leftrightarrow\nu)\right].$
(S34)

Therefore, we have shown the decomposition for the average currents

$\braket{Q}=\braket{Q_{\rm dyn}}+\braket{Q_{\rm geo}}=\braket{Q_{\rm
dyn}}+\braket{Q_{\rm curv}}+\braket{Q_{\rm metr}}.$ (S35)

We note here that the above derived average current components are part of our
general geometric formulation. Our CGF theory is indispensable in analyzing
the fluctuations of finite-time periodically driven systems.

##  IV Non-adiabatic control over current pump

Here, we discuss the non-adiabatic geometric control over current pump effect.

The metric $\mathfrak{g}_{\mu\nu}^{Q}$ for an arbitrary current $Q$ is not
guaranteed to be positive definite. It can have both positive and negative
eigenvalues. If we pass through the positive (negative) eigenvector
directions, the non-adiabatic geometric current is along its positive
(negative) direction. Between the negative and positive eigenvectors, there
must exist a direction where
$\mathfrak{g}_{\mu\nu}\dot{\Lambda}_{\mu}\dot{\Lambda}_{\nu}=0$ is satisfied.
If we design our driving protocol to be along these zero-vector directions, we
can eliminate the non-adiabatic pump effect. This constructs a non-adiabatic
geometric control over the current pump effect.

Here, we demonstrate this control principle with a two-level-system (TLS)
model in Fig. S1. It is composed of a TLS of energy difference $\omega$ and is
coupled to two reservoirs of the same temperature $T$. We drive the parameter
$\bm{\Lambda}:=(\omega,T)^{\rm T}$. By counting the heat current $Q$ from the
left reservoir into the TLS, we obtain the twisted master equation operator

$\hat{L}_{\chi}=\begin{pmatrix}-2n&(1+n)(e^{-\chi\omega}+1)\\\
n(e^{\chi\omega}+1)&-2(1+n)\end{pmatrix},$ (S36)

where $\chi$ is the counting parameter and $n=1/(e^{\omega/T}-1)$ is the Bose-
Einstein distribution. Using Eq. (S20), we obtain the metric in parameter
space as

$\mathfrak{g}^{Q}=\partial_{\chi}\mathfrak{g}|_{\chi=0}=\frac{1}{4T^{2}}{\rm
csch}^{3}(\frac{\omega}{T})\sinh^{4}(\frac{\omega}{2T})\begin{pmatrix}-2T&\omega\\\
\omega&0\end{pmatrix}.$ (S37)

It indeed has both positive and negative eigenvalues and we illustrate its
corresponding zero-vector direction field in Fig. S1. This non-adiabatic
control principle can also be used to regulate other types of generic currents
and their cumulants.

Figure S1:  The non-adiabatic geometric control over current pump. In the
parameter space composed of system energy level difference $\omega$ and the
bath temperature $T$, the arrows illustrates the direction of zero non-
adiabatic directions, along which
$\mathfrak{g}_{\mu\nu}^{Q}\dot{\Lambda}_{\mu}\dot{\Lambda}_{\nu}=0$ is
satisfied. By designing contours along this zero direction field, as shown in
the red curve, we can eliminate the non-adiabatic effect of current pump
effect.

##  V Derivation of Geometric TURs

The geometric TUR for cyclically driven thermal device is derived in this
section.

We set the counting parameter $\chi=0$ in this proof. Starting from
$\partial_{t}\ket{p(t)}=\hat{L}(t)\ket{p(t)}$, we obtain a self-consistent
equation for $\ket{\delta(t)}=\ket{p(t)}-\ket{\pi(t)}$:

$\ket{\delta(t)}=\hat{L}^{+}\ket{\partial_{t}\pi(t)}+\hat{L}^{+}\ket{\partial_{t}\delta},$
(S38)

where $\ket{\pi}$ is the instantaneous steady state satisfying
$\hat{L}(t)\ket{\pi(t)}=0$ and $\hat{L}^{+}$ is the pseudo-inverse of
$\hat{L}$. $\hat{L}^{+}$ is defined by
$\hat{L}^{+}\hat{L}=\hat{L}\hat{L}^{+}=\hat{I}-\ket{\pi}\bra{1}$, with
$\hat{I}$ being the identity operator and $\bra{1}$ the one-vector. By
recurrently substituting, it is obvious to obtain the adiabatic and
nonadiabatic corrections to the distribution $\ket{p}=\ket{\pi}+\ket{p_{\rm
curv}}+\ket{p_{\rm metr}}$ respectively, where $\ket{p_{\rm
curv}}=\hat{L}^{+}\ket{\partial_{t}\pi(t)}$ and $\ket{p_{\rm
metr}}=\hat{L}^{+}\partial_{t}(\hat{L}^{+}\ket{\partial_{t}\pi(t)})$.

The current operator $\hat{J}$ is defined by
$J_{ij}=\sum_{\mu}k_{ij}^{\mu}d_{ij}^{\mu}$ in terms of the transition rates
$k_{ij}^{\mu}$ and the antisymmetric tensor $d_{ij}^{\mu}=-d_{ji}^{\mu}$. The
component $d_{ij}^{\mu}$ is the stochastic current accompanying the transition
$j\to i$ induced by the reservoir $\mu$. In these notations, the adiabatic and
nonadiabatic geometric currents are

$\begin{split}\braket{Q_{\rm
curv}}&=\int_{\Omega}dS_{\mu\nu}F_{\mu\nu}^{Q},\\\ \braket{Q_{\rm
metr}}&=\int_{0}^{\tau_{p}}dt\dot{\Lambda}_{\mu}\dot{\Lambda}_{\nu}\mathfrak{g}_{\mu\nu}^{Q},\end{split}$
(S39)

where
$F_{\mu\nu}^{Q}=\braket{1}{\partial_{\mu}(\hat{J}\hat{L}^{+})}{\partial_{\nu}\pi}-\braket{1}{\partial_{\nu}(\hat{J}\hat{L}^{+})}{\partial_{\mu}\pi}$,
$\mathfrak{g}_{\mu\nu}^{Q}=\frac{1}{2}\left[\bra{1}\hat{J}\hat{L}^{+}\partial_{\mu}(\hat{L}^{+}\ket{\partial_{\nu}\pi})+(\mu\leftrightarrow\nu)\right]$.
We note here that these expression for the average current can also be derived
from Eq. (S14) and Eq. (S20) by taking derivative to $\chi$ and then setting
$\chi\to 0$, as shown in Sec. III.

Now, let us start our proof of Geometric TURs from the fluctuation-response
inequality [2], which is equivalent to the classical Cramer-Rao bound in the
linear response regime. Similar method has been adopted by [3, 4] in deriving
their TURs. Our method here is similar to those used by Ref. [2] in proving
their generalized TURs. We consider a virtual perturbation on the transition
rates, with $\theta$ being its perturbation strength. Considering the
Kullback-Leibler (KL) divergence of the perturbed and unperturbed trajectory
weight to the order of $\theta^{2}$, one can write

$\begin{split}D_{\rm KL}&=\int d\omega
P_{\theta}(\omega)\ln(\frac{P_{\theta}(\omega)}{P_{0}(\omega)})\\\
&=\frac{\theta^{2}}{2}\int d\omega P_{0}(\omega)[\partial_{\theta}\ln
P_{\theta}(\omega)|_{\theta\to 0}]^{2}+{\mathcal{O}}(\theta^{3})\\\
&=-\frac{\theta^{2}}{2}\int d\omega P_{0}(\omega)[\partial_{\theta}^{2}\ln
P_{\theta}(\omega)|_{\theta\to 0}]+{\mathcal{O}}(\theta^{3})\\\
&=-\frac{\theta^{2}}{2}\braket{\partial_{\theta}^{2}\ln
P_{\theta}}{{}_{\theta\to 0}}_{0}+{\mathcal{O}}(\theta^{3}).\end{split}$ (S40)

For the discrete Markovian systems, the weight for the trajectory
$\omega={x(t)}$ is given by [3]

$P_{\theta}(\omega)=p(x_{0})e^{-\int_{0}^{\tau_{p}}dt\sum_{i\neq
j;\mu}[\delta_{x(t),i}k_{ji}^{\mu}(\theta,t)-\dot{m}_{ji}^{\mu}\ln
k_{ji}^{\mu}(\theta,t)]},$ (S41)

where $p(x_{0})$ is the initial condition and $m_{ij}^{\mu}$ is the
accumulated number of transitions. The basic average quantities are given by
$\braket{\delta_{x(t),i}}_{\theta}=p_{i}(\theta,t)$ and
$\braket{\dot{m}_{ij}^{\mu}}_{\theta}=p_{j}(\theta,t)k_{ij}^{\mu}(\theta,t)$.
We choose the form of the virtual perturbation to be
$k_{ij}^{\mu}(\theta,t)=k_{ij}^{\mu}e^{\theta\alpha_{ij}^{\mu}(t)}$. By
showing $\partial_{\theta}^{2}\ln P_{\theta}|_{\theta\to 0}=-\sum_{i\neq
j;\mu}\int_{0}^{\tau_{p}}dt\delta_{x(t),i}k_{ji}^{\mu}(\alpha_{ji}^{\mu})^{2}$,
we express the KL divergence as

$\begin{split}D_{\rm
KL}&\approx-\frac{\theta^{2}}{2}\braket{\partial_{\theta}^{2}\ln
P_{\theta}}{{}_{\theta\to 0}}_{0}\\\ &=\frac{\theta^{2}}{2}\sum_{i\neq
j;\mu}\int_{0}^{\tau_{p}}dt\braket{\delta_{x(t),i}}_{0}k_{ji}^{\mu}(t)[\alpha_{ji}^{\mu}(t)]^{2}\\\
&=\frac{\theta^{2}}{2}\sum_{i\neq
j;\mu}\int_{0}^{\tau_{p}}dtp_{i}(t)k_{ji}^{\mu}(t)[\alpha_{ji}^{\mu}(t)]^{2}\\\
&=\frac{\theta^{2}}{2}\sum_{i>j;\mu}\int_{0}^{\tau_{p}}dt\frac{[j_{ji}^{\mu}(t)]^{2}}{t_{ji}^{\mu}(t)}.\end{split}$
(S42)

Here we select $\alpha_{ij}^{\mu}=j_{ij}^{\mu}/t_{ij}^{\mu}$ in the last line,
where we define the detailed current
$j_{ij}^{\mu}=p_{j}k_{ij}^{\mu}-p_{i}k_{ji}^{\mu}$ and the activity
$t_{ij}^{\mu}=p_{j}k_{ij}^{\mu}+p_{i}k_{ji}^{\mu}$. By using the log-sum
inequality $\frac{2(b-a)^{2}}{a+b}\leq(b-a)(\ln b-\ln a)$ valid for an
arbitrary pair of $a$ and $b$, $D_{\rm KL}$ provides a lower bound on the
entropy production $\braket{\Sigma}_{0}=\int_{0}^{\tau_{p}}dt\sum_{i\neq
j;\mu}p_{j}k_{ij}^{\mu}\ln[(p_{j}k_{ij}^{\mu})/(p_{i}k_{ji}^{\mu})]$, as given
by Ref. [2],

$D_{\rm KL}\leq\frac{\theta^{2}}{4}\braket{\Sigma}_{0}.$ (S43)

Also, since $(b-a)^{2}/(a+b)\leq a+b$ for positive $a$ and $b$, the KL
divergence is also bounded by

$D_{\rm KL}\leq\frac{\theta^{2}}{2}\braket{A}_{0},$ (S44)

with the average activity being
$\braket{A}_{0}=\int_{0}^{\tau_{p}}dt\sum_{i\neq
j;\mu}(p_{j}k_{ij}^{\mu}+p_{i}k_{ji}^{\mu})$.

Using the form of $\alpha_{ij}^{\mu}$, it is simple to show that [2], to the
first order of $\theta$,

$\begin{split}\hat{L}(\theta,t)\ket{p(0,t)}&=(1+\theta)\hat{L}(0,t)\ket{p(0,t)}+\mathcal{O}(\theta^{2})\\\
&=(1+\theta)\partial_{t}\ket{p(0,t)}+\mathcal{O}(\theta^{2}).\end{split}$
(S45)

By applying $\hat{L}^{+}(\theta,t)$ on both sides, we have

$\ket{p(0,t)}=\ket{\pi(\theta,t)}+(1+\theta)\hat{L}^{+}(\theta,t)\ket{\partial_{t}p(0,t)}.$
(S46)

By applying the adiabatic perturbation method, we derive
$\ket{p(0,t)}=\ket{\pi(\theta,t)}+(1+\theta)\hat{L}^{+}(\theta,t)\ket{\pi(\theta,t)}+(1+\theta)^{2}\hat{L}^{+}(\theta,t)\partial_{t}(\hat{L}^{+}(\theta,t)\ket{\partial_{t}\pi(\theta,t)})+{\mathcal{O}}(\theta^{2})=\ket{p(\theta,t)}+\theta\ket{p_{\rm
curv}(\theta,t)}+2\theta\ket{p_{\rm
metr}(\theta,t)}+{\mathcal{O}}(\theta^{2})$. Therefore, we have

$\begin{split}\partial_{\theta}\ket{p(\theta,t)}|_{\theta\to
0}&=\lim_{\theta\to 0}\frac{\ket{p(\theta,t)}-\ket{p(0,t)}}{\theta}\\\
&=-\ket{p_{\rm curv}(0,t)}-2\ket{p_{\rm metr}(0,t)}.\end{split}$ (S47)

Also, we can obtain
$\braket{1}{\partial_{\theta}\hat{J}(\theta)}{{}_{\theta\to
0}|p(0,t)}=\braket{1}{\hat{J}(0,t)}{p(0,t)}$ by using the form of the
perturbed transition rates $k_{ij}^{\mu}(\theta,t)$.

In this way, we show that the linear response of $\braket{Q}$ due to the
perturbation is given by

$\begin{split}\lim_{\theta\to
0}\frac{\braket{Q}_{\theta}-\braket{Q}_{0}}{\theta}&=\int_{0}^{\tau_{p}}dt\partial_{\theta}\braket{1}{\hat{J}(\theta,t)}{p(\theta,t)}|_{\theta\to
0}\\\ &=\braket{Q_{\rm dyn}}_{0}-\braket{Q_{\rm metr}}_{0}.\end{split}$ (S48)

According to the linear fluctuation-response inequality
$(\braket{Q}_{\theta}-\braket{Q}_{0})^{2}\leq 2D_{\rm KL}\braket{Q^{2}}_{c}$,
first propose by Ref. [2], we reach the final result

$\braket{\Sigma}\braket{Q^{2}}_{c}\geq 2(\braket{Q_{\rm dyn}}-\braket{Q_{\rm
metr}})^{2},$ (S49)

where we omit the subscript $0$ since all averages are with regard to the
unperturbed dynamics. The average entropy production is thus bounded by
$\Sigma_{g}$

$\braket{\Sigma}\geq\frac{2(\braket{Q_{\rm dyn}}-\braket{Q_{\rm
metr}})^{2}}{\braket{Q^{2}}_{c}}:=\Sigma_{\rm g}.$ (S50)

This is the so called Geometric TURs, one of our main results.

Similarly, the average activity is bounded by

$\braket{A}\geq\frac{(\braket{Q_{\rm dyn}}-\braket{Q_{\rm
metr}})^{2}}{\braket{Q^{2}}_{c}}:=A_{\rm g}.$ (S51)

We further discuss some consequences of Eq. (S49). If the instantaneous steady
state $\ket{\pi(t)}$ is actually the equilibrium state, the dynamic part
$\braket{\Sigma_{\rm dyn}}$ and $\braket{Q_{\rm dyn}}$ and the adiabatic
curvature part $\braket{\Sigma_{\rm curv}}$ vanish. The TUR simplifies to

$\braket{\Sigma_{\rm metr}}\braket{Q^{2}}_{c}\geq 2\braket{Q_{\rm metr}}^{2}.$
(S52)

If we consider the current $Q$ to be $\Sigma$, this inequality further
simplifies to

$\braket{\Sigma^{2}}_{c}\geq 2\braket{\Sigma_{\rm
metr}}\geq\frac{2{\mathcal{L}}^{2}}{\tau_{p}},$ (S53)

bounding the fluctuation of entropy production with the thermodynamic length
given by the metric of $\braket{\Sigma_{\rm metr}}$, i.e., for
$\braket{\Sigma_{\rm
metr}}=\int_{0}^{\tau_{p}}dtg_{\mu\nu}^{\Sigma}\dot{\Lambda}_{\mu}\dot{\Lambda}_{\nu}$,
${\mathcal{L}}:=\oint_{\partial\Omega}dt\sqrt{g_{\mu\nu}^{\Sigma}\dot{\Lambda}_{\mu}\dot{\Lambda}_{\nu}}$.
Here, the second inequality can actually saturated by optimizing the
parametrization of protocol. We note ${\mathcal{L}}$ is independent of the
detailed protocol but only dependent on the path in the parameter space.
According to this inequality, the underlying thermodynamic length is bounded
from above by an arbitrary parametrization of the path in the parameter space

${\mathcal{L}}\leq\sqrt{\tau_{p}\braket{\Sigma^{2}}_{c}/2}.$ (S54)

This form a basis for future inference of the thermodynamic length.

##  VI Models

We treat the theoretical details of our two models with continuous and
discrete degrees of freedom in this section.

###  VI.1 Chiral current in the nonequilibrium tricycle

We consider a system composed of three quantum dots with tunable energy
levels. Electrons can tunnel between one of the dots and an electron thermal
reservoir. We also indirectly couple these dots with three photonic/phononic
baths. The bosonic baths compensate the energy difference between dots and
renders the system a nonequilibrium quantum tricycle. To be specific, the
system Hamiltonian
$\hat{H}=\hat{H}_{S}+\hat{H}_{R}+\hat{H}_{SR}+\hat{H}_{B}+\hat{H}_{SB}$ is
composed of the three quantum dot levels
$\hat{H}_{S}=\sum_{n=1}^{3}\epsilon_{n}\hat{c}_{n}^{\dagger}\hat{c}_{n}$, the
electron reservoirs
$\hat{H}_{R}=\sum_{k}\epsilon_{k}\hat{d}_{k}^{\dagger}\hat{d}_{k}$, the
tunneling term
$\hat{H}_{SR}=\sum_{k}t_{k}(\hat{d}_{k}^{\dagger}\hat{c}_{3}+\hat{c}_{3}^{\dagger}\hat{d}_{k})$,
the Bosonic thermal bath
$\hat{H}_{B}=\sum_{\nu=1;k}^{\nu=3}\epsilon_{\nu,k}\hat{a}_{\nu,k}^{\dagger}\hat{a}_{\nu,k}$,
and the system-bath coupling term
$\hat{H}_{SB}=\sum_{\nu=1;k}^{\nu=3}r_{\nu,k}(\hat{a}_{\nu,k}+\hat{a}_{\nu,k}^{\dagger})(\hat{c}_{\nu}^{\dagger}\hat{c}_{\nu+1}+\hat{c}_{\nu+1}^{\dagger}\hat{c}_{\nu})$.
Here, $\nu=4$ denotes the same site as $\nu=1$. We restrict ourselves to the
weak coupling regime. Under the standard Born-Markov approximation, we obtain
the transition rates between different states.

By concentrating on the strong Coulomb blockade regime, the master equation is
$\partial_{t}\ket{p_{\bm{\chi}}}=\hat{L}_{\bm{\chi}}\ket{p_{\bm{\chi}}}$, with
$\ket{p}=(p_{0},p_{1},p_{2},p_{3})^{\rm T}$ and

$\hat{L}_{\bm{\chi}}=\begin{pmatrix}-k_{+}&0&0&k_{-}e^{\chi_{2}(\epsilon_{3}-\mu)/T_{R}}\\\
0&-(k_{21}+k_{31})&k_{12}e^{\chi_{2}(\epsilon_{2}-\epsilon_{1})/T_{1}}&k_{13}e^{-\chi_{1}+\chi_{2}(\epsilon_{3}-\epsilon_{1})/T_{3}}\\\
0&k_{21}e^{-\chi_{2}(\epsilon_{2}-\epsilon_{1})/T_{1}}&-(k_{12}+k_{32})&k_{23}e^{\chi_{2}(\epsilon_{3}-\epsilon_{2})/T_{2}}\\\
k_{+}e^{-\chi_{2}(\epsilon_{3}-\mu)/T_{R}}&k_{31}e^{\chi_{1}-\chi_{2}(\epsilon_{3}-\epsilon_{1})/T_{3}}&k_{32}e^{-\chi_{2}(\epsilon_{3}-\epsilon_{2})/T_{2}}&-(k_{-}+k_{13}+k_{23})\end{pmatrix},$
(S55)

where $k_{ij}$ represents the transition from $j$ to $i$. To be concrete, we
suppose $\epsilon_{3}>\epsilon_{2}>\epsilon_{1}$. Here, the rates
$\Gamma_{R}=2\pi\sum_{k}t_{k}^{2}\delta(\epsilon-\epsilon_{k})$ and
$\Gamma_{\nu}=2\pi\sum_{k}r_{\nu,k}^{2}\delta(\epsilon-\epsilon_{\nu,k})$ are
taken as constant under the flat-band limit.
$f(\epsilon)=1/[e^{(\epsilon-\mu)/T_{R}}+1]$ and
$n_{\nu}(\epsilon)=1/(e^{\epsilon/T_{\nu}}-1)$ are respectively the Fermi and
Bose distribution.

$\begin{split}k_{+}&=\Gamma_{R}f(\epsilon_{3}),\\\
k_{-}&=\Gamma_{R}[1-f(\epsilon_{3})],\end{split}$ (S56)

$\begin{split}k_{21}&=\Gamma_{1}n_{1}(\epsilon_{2}-\epsilon_{1}),\\\
k_{12}&=\Gamma_{1}[1+n_{1}(\epsilon_{2}-\epsilon_{1})],\end{split}$ (S57)

$\begin{split}k_{32}&=\Gamma_{2}n_{2}(\epsilon_{3}-\epsilon_{2}),\\\
k_{23}&=\Gamma_{2}[1+n_{2}(\epsilon_{3}-\epsilon_{2})],\end{split}$ (S58)

$\begin{split}k_{31}&=\Gamma_{3}n_{3}(\epsilon_{3}-\epsilon_{1}),\\\
k_{13}&=\Gamma_{3}[1+n_{3}(\epsilon_{3}-\epsilon_{1})].\end{split}$ (S59)

We use $\chi_{1}$ to count the chiral current and $\chi_{2}$ to count the
entropy production.

Similarly, with regard to the dynamic activity of the whole system, which is
the total number of bidirectional transitions between each pair of states, we
can also count it with $\chi_{3}$, with the corresponding

$\hat{L}_{\bm{\chi}}=\begin{pmatrix}-k_{+}&0&0&k_{-}e^{\chi_{3}}\\\
0&-(k_{21}+k_{31})&k_{12}e^{\chi_{3}}&k_{13}e^{\chi_{3}}\\\
0&k_{21}e^{\chi_{3}}&-(k_{12}+k_{32})&k_{23}e^{\chi_{3}}\\\
k_{+}e^{\chi_{3}}&k_{31}e^{\chi_{3}}&k_{32}e^{\chi_{3}}&-(k_{-}+k_{13}+k_{23})\end{pmatrix},$
(S60)

We drive the system according to the protocol
$\epsilon_{3}(t)=1+\delta\sin(2\pi t/\tau_{p})$,
$\epsilon_{2}(t)=0.5+\delta\sin(2\pi t/\tau_{p}+\phi)$,
$\epsilon_{1}(t)=\delta\sin(2\pi t/\tau_{p}+2\phi)$ ($0\leq\delta<0.25$), and
$\Gamma_{\nu}=\Gamma_{R}=\Gamma$. The parameters $\delta$, $\phi$, $\tau_{p}$,
$\Gamma_{B}$, $T_{\nu}$ ($\nu=1,2,3$), $T_{R}$ and $\mu$ determine our model.
The characteristic time scale of this system is $\tau_{c}\approx 1/\Gamma$.
The pumped chiral current is shown in Fig. S2. $T_{3}=1.00$ corresponds to the
case with no dynamic components.

In the Fig. 2(c) of the main text, we select $T_{1}=T_{2}=T_{R}=1,\mu=5$, the
period $10^{1/2}\leq\tau_{p}\leq 100$, the driving amplitude $\delta=0.2$, the
temperature $0.97\leq T_{3}\leq 1.03$ and the phase difference $0\leq\phi\leq
2\pi$.

In the Fig. 2(d) of the main text, we select the temperatures all being
$T_{1}=T_{2}=T_{3}=T_{R}=1$, $\mu=5$, the driving amplitude
$\delta\in\\{0.05,0.10,0.15,0.20\\}$ and the phase difference $0\leq\phi\leq
2\pi$.

Figure S2: The nonequilibrium chiral current versus the driving period.
$T_{3}=1.00$ has no dynamic component and other two systems have nonzero
dynamic chiral current. Here, the parameters are given by $\phi=2\pi/3$,
$\delta=0.2$, $T_{1}=T_{2}=T_{R}=1.00$, and $\mu=5.0$.

###  VI.2 Brownian heat pump

Figure S3:  The Brownian heat pump model by cyclically driving the stiffness
of the oscillator $k$ and the coupling to the left reservoir $\gamma_{1}$. (a)
The Brownian model and its effective system. (b) The geometrically pumped heat
density $Q/\tau$ versus the inverse period $1/\tau_{\rm p}$ (the solid line).
Here, we choose the protocol $k(t)/k_{0}=1+\sin(2\pi t/\tau_{\rm
p})/2,\gamma_{1}(t)/\gamma_{2}=1+\sin(2\pi t/\tau_{\rm p}+\pi/2)/2$ and
$T_{1}=T_{2}$, with $k_{0}$ being a stiffness constant. The dot-dash line is
for the adiabatic component $\braket{Q_{\rm curv}}/\tau_{\rm p}$ and the dash
line denotes the optimal driving period. (c) The thermodynamic uncertainty
relations ($\Sigma_{\alpha}=\Sigma_{\rm b}$ or $\Sigma_{\rm g}$). Our
geometric bound $\braket{\Sigma}\geq\Sigma_{\rm g}$ is satisfied (blue
circles), while the steady state bound $\braket{\Sigma}\geq\Sigma_{\rm b}$ can
be broken (orange triangles); we only show the breaking situations of the
latter bound for clarity. (d) The geometric bound of the fluctuation of the
entropy production (Eq. (S53)). Blue circles are for constant speed driving
protocols and green triangles (on the dashed line) are for the optimal
protocols parametrized in terms of the length ${\mathcal{L}}$.

Consider a harmonic oscillator composed of two particles coupling two
reservoirs, which is illustrated in Fig. S3(a). The Langevin equation of
motion is given by

$\begin{split}\gamma_{1}\dot{x}_{1}&=k(x_{2}-x_{1})+\xi_{1},\\\
\gamma_{2}\dot{x}_{2}&=k(x_{1}-x_{2})+\xi_{2},\end{split}$ (S61)

where $k$ is the stiffness of the oscillator and $\xi_{i}$ is the Gaussian
white noise satisfying
$\braket{\xi_{i}(t_{1})\xi_{j}(t_{2})}=2\gamma_{i}T_{i}\delta_{ij}\delta(t_{1}-t_{2})$.
Due to the translational invariance of the above equations, it is equivalent
to the dynamics of a single overdamped oscillator $\gamma\dot{y}=-ky+\xi(t)$,
with $y=x_{1}-x_{2}$ being the effective degree of freedom,
$\gamma=\gamma_{1}\gamma_{2}/(\gamma_{1}+\gamma_{2})$, and
$\xi=(\gamma_{2}\xi_{1}-\gamma_{1}\xi_{2})/(\gamma_{1}+\gamma_{2})$. We note
the effective noise $\xi$ has zero mean and the variance
$\braket{\xi(t_{1})\xi(t_{2})}=2\gamma T_{\rm e}\delta(t_{1}-t_{2})$, with the
effective temperature given by $T_{\rm
e}=(\gamma_{2}T_{1}+\gamma_{1}T_{2})/(\gamma_{1}+\gamma_{2})$.

Here, the stochastic heat flowing into the system from the first reservoir is
$\dot{Q}_{1}=(-\gamma_{1}\dot{x}_{1}+\xi_{1})\dot{x}_{1}=-\frac{k^{2}}{\gamma_{1}}y^{2}+\frac{k}{\gamma_{1}}\xi_{1}y$,
and the stochastic entropy production rate
$\dot{\Sigma}=-\dot{Q}_{1}/T_{1}-\dot{Q}_{2}/T_{2}$ is given by

$\dot{\Sigma}=(\frac{1}{\gamma_{1}T_{1}}+\frac{1}{\gamma_{2}T_{2}})k^{2}y^{2}-(\frac{\xi_{1}}{\gamma_{1}T_{1}}-\frac{\xi_{2}}{\gamma_{2}T_{2}})ky.$
(S62)

Defining the joint distribution function $\ket{p(y,Q,\Sigma,t)}$, its dynamics
is governed by the stochastic Liouville equation
$\partial_{t}\ket{p(y,Q,\Sigma,t)}=\hat{\Omega}\ket{p(y,Q,\Sigma,t)}$
following the continuity relation
$\hat{\Omega}=-\partial_{y}\dot{y}-\partial_{Q_{1}}\dot{Q}_{1}-\partial_{\Sigma}\dot{\Sigma}$.
By taking an ensemble average over the reservoir noise, we derive the
generalized Fokker-Planck operator as $\hat{L}=\lim_{\Delta t\to
0}\frac{1}{\Delta t}[\int_{t}^{t+\Delta
t}dt_{1}\braket{\hat{\Omega}(t_{1})}+\int_{t}^{t+\Delta
t}dt_{1}\int_{t}^{t_{1}}dt_{2}\braket{\hat{\Omega}(t_{1})\hat{\Omega}(t_{2})}]$.
The first term is the ballistic term
$\hat{L}_{1}=\frac{k}{\gamma}\partial_{y}y+\frac{k^{2}}{\gamma_{1}}y^{2}\partial_{Q_{1}}-(\frac{1}{\gamma_{1}T_{1}}+\frac{1}{\gamma_{2}T_{2}})k^{2}y^{2}\partial_{\Sigma}$,
while the second term is the diffusion term $\hat{L}_{2}=\frac{T_{\rm
e}}{\gamma}\partial_{y}^{2}+\frac{kT_{1}}{\gamma_{1}}(\partial_{y}y+y\partial_{y})\partial_{Q_{1}}-\frac{k}{\gamma}(\partial_{y}y+y\partial_{y})\partial_{\Sigma}-\frac{2k^{2}}{\gamma_{1}}y^{2}\partial_{Q_{1}}\partial_{\Sigma}+\frac{k^{2}T_{1}}{\gamma_{1}}y^{2}\partial_{Q_{1}}^{2}+(\frac{1}{\gamma_{1}T_{1}}+\frac{1}{\gamma_{2}T_{2}})k^{2}y^{2}\partial_{\Sigma}^{2}$.

We use the counting field ${\bm{\chi}}=(\chi_{1},\chi_{2})^{\rm T}$ with two
components generating the statistics of $Q_{1}$ and $\Sigma$ respectively. By
making a Fourier-Laplace transformation $\ket{p(y,{\bm{\chi}},t)}=\int
dQd\Sigma\ket{p(y,Q,\Sigma,t)}e^{\chi_{1}Q+\chi_{2}\Sigma}$, we derive the
twisted Fokker-Planck operator as a generator for the dynamics of
$\ket{p(y,{\bm{\chi}},t)}$

$\begin{split}\hat{L}_{\bm{\chi}}=\frac{T_{\rm
e}}{\gamma}\partial_{y}^{2}+a_{1}y\partial_{y}+a_{2}y^{2}+a_{3},\end{split}$
(S63)

with the coefficients
$a_{1}=\frac{k}{\gamma}[1+2(\chi_{2}-\frac{\gamma}{\gamma_{1}}T_{1}\chi_{1})]$,
$a_{2}=k^{2}[\frac{\chi_{1}}{\gamma_{1}}(T_{1}\chi_{1}-1)+(\frac{1}{\gamma_{1}T_{1}}+\frac{1}{\gamma_{2}T_{2}})\chi_{2}(\chi_{2}+1)-\frac{2}{\gamma_{1}}\chi_{1}\chi_{2}]$
and
$a_{3}=\frac{k}{\gamma}(1+\chi_{2}-\frac{\gamma}{\gamma_{1}}T_{1}\chi_{1})$. A
similarity transformation $\hat{U}=e^{\beta ky^{2}/2}$ would bring
$\hat{L}_{\bm{\chi}}$ into a Hermitian operator
$\tilde{L}_{\bm{\chi}}=\tilde{L}_{\bm{\chi}}^{\dagger}=\hat{U}\hat{L}_{\bm{\chi}}\hat{U}^{-1}=\frac{T_{\rm
e}}{\gamma}\partial_{y}^{2}+(a_{2}-\frac{a_{1}^{2}\gamma}{4T_{\rm
e}})y^{2}+(a_{3}-\frac{a_{1}}{2})$, where $\beta=\frac{a_{1}\gamma}{2T_{\rm
e}k}$.

Analogous to the solution of quantum oscillators, in terms of the lowering
$\hat{b}=\sqrt{\epsilon}\partial_{y}+\frac{1}{2\sqrt{\epsilon}}y$ and raising
operator
$\hat{b}^{\dagger}=-\sqrt{\epsilon}\partial_{y}+\frac{1}{2\sqrt{\epsilon}}y$
($[\hat{b},\hat{b}^{\dagger}]=1$), where the factor $\epsilon=\frac{T_{\rm
e}}{\sqrt{a_{1}^{2}\gamma^{2}-4a_{2}T_{\rm e}\gamma}}$, we can write
$\tilde{L}_{\bm{\chi}}=-\frac{T_{\rm
e}}{\gamma\epsilon}\hat{b}^{\dagger}\hat{b}+(a_{3}-\frac{a_{1}}{2}-\frac{T_{\rm
e}}{2\gamma\epsilon})$. Therefore, we derive the eigenvalues of
$\hat{L}_{\bm{\chi}}$ to be $E_{n}=-n\frac{T_{\rm
e}}{\gamma\epsilon}+(a_{3}-\frac{a_{1}}{2}-\frac{T_{\rm e}}{2\gamma\epsilon})$
($E_{0}$ corresponding to the steady state) and the eigenvectors to be

$\begin{split}\ket{r_{n}}&=\hat{U}^{-1}\ket{\psi_{n}}=\frac{1}{\sqrt{n!}}\hat{U}^{-1}(\hat{b}^{\dagger})^{n}\ket{\psi_{0}},\\\
\bra{l_{n}}&=\bra{\psi_{n}}\hat{U}=\frac{1}{\sqrt{n!}}\bra{\psi_{0}}\hat{b}^{n}\hat{U},\end{split}$
(S64)

where $\ket{\psi_{0}}=(\bra{\psi_{0}})^{\dagger}$ is the ground eigenstate of
$\tilde{L}_{\bm{\chi}}$ satisfying $\hat{b}\ket{\psi_{0}}=0$. Here,
$\\{\ket{\psi_{n}}\\}$ forms an orthonormal basis
$\braket{\psi_{m}}{\psi_{n}}=\delta_{mn}$.

The ground state satisfies $\hat{b}\ket{\psi_{0}}=0$
($(\sqrt{\epsilon}\partial_{y}+\frac{1}{2\sqrt{\epsilon}}y)\ket{\psi_{0}}=0$)
and is thus given by a Gaussian distribution [5]

$\ket{\psi_{0}}=(\frac{1}{2\pi\epsilon})^{1/4}e^{-\frac{y^{2}}{4\epsilon}},$
(S65)

with an expression similar to the ground state wavefunction of the quantum
harmonic oscillator. To calculate the geometric CGF, we note the relations
$y=\sqrt{\epsilon}(\hat{b}+\hat{b}^{\dagger})$,
$\ket{\partial_{\mu}\psi_{0}}=\frac{(y^{2}-\epsilon)\partial_{\mu}\epsilon}{4\epsilon^{2}}\ket{\psi_{0}}$,
$\partial_{\mu}\hat{U}=(\partial_{\mu}\beta k)\frac{y^{2}}{2}\hat{U}$ and
$\partial_{\mu}\hat{U}^{-1}=-(\partial_{\mu}\beta
k)\frac{y^{2}}{2}\hat{U}^{-1}$. We show

$\begin{split}\ket{\partial_{\mu}r_{0}}&=-\frac{\epsilon\partial_{\mu}(\beta
k)}{2}\ket{r_{0}}+\sqrt{2}[\frac{\partial_{\mu}\epsilon}{4\epsilon}-\frac{\epsilon\partial_{\mu}(\beta
k)}{2}]\ket{r_{2}},\\\
\bra{\partial_{\mu}l_{0}}&=\frac{\epsilon\partial_{\mu}(\beta
k)}{2}\bra{l_{0}}+\sqrt{2}[\frac{\partial_{\mu}\epsilon}{4\epsilon}+\frac{\epsilon\partial_{\mu}(\beta
k)}{2}]\bra{l_{2}}.\end{split}$ (S66)

The relation
$\braket{l_{0}}{\partial_{\mu}r_{0}}+\braket{\partial_{\mu}l_{0}}{r_{0}}=0$ is
satisfied due to the normalization condition $\braket{l_{0}}{r_{0}}=1$. For
arbitrary driven parameters, the geometric connection is
$A_{\mu}=-\braket{l_{0}}{\partial_{\mu}r_{0}}=\epsilon\partial_{\mu}(\beta
k)/2$. The corresponding geometric curvature is
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}=[(\partial_{\mu}\epsilon)(\partial_{\nu}\beta
k)-(\partial_{\nu}\epsilon)(\partial_{\mu}\beta k)]/2$. The metric tensor is
given by

$g_{\mu\nu}=\frac{\gamma}{16T_{\rm
e}\epsilon}[4\epsilon^{4}(\partial_{\mu}\beta k)(\partial_{\nu}\beta
k)-(\partial_{\mu}\epsilon)(\partial_{\nu}\epsilon)].$ (S67)

To illustrate our theory, we consider a heat engine powered by driving $k$ and
$\gamma_{1}$, i.e. ${\bm{\Lambda}}=(k,\gamma_{1})^{\rm T}$. We first consider
the driving protocol $k=k_{0}[1+a\sin(2\pi t/\tau_{p})]$,
$\gamma_{1}=\gamma_{0}[1+a\sin(2\pi t/\tau_{p}+\phi)]$, and
$\gamma_{2}=\gamma_{0}$. The parameters $\tau_{p}$, $T_{1}$, $T_{2}$ and
$\phi$ can be taken as configurations of the driving protocol. We further
define a time scale of this system as $\tau_{c}=\gamma_{0}/k_{0}$, which is of
the same order of the characteristic time scale $\gamma/k$.

The instantaneous dynamic components for $\braket{Q}$, $\braket{Q^{2}}_{c}$,
$\braket{\Sigma}$ and $\braket{\Sigma^{2}}_{c}$ are given by

$\begin{split}E^{Q}&:=\partial_{\chi_{1}}E_{0}|_{{\bm{\chi}}=0}=\frac{1}{\tau_{c}}\frac{\tilde{k}(T_{1}-T_{2})}{\tilde{\gamma}_{1}+1},\end{split}$
(S68)
$\begin{split}E^{Q^{2}}&:=\partial_{\chi_{1}}^{2}E_{0}|_{{\bm{\chi}}=0}=\frac{1}{\tau_{c}}\frac{2\tilde{k}(T_{2}\tilde{\gamma}_{1}+T_{1})(T_{1}\tilde{\gamma}_{1}+T_{2})}{(\tilde{\gamma}_{1}+1)^{3}},\end{split}$
(S69)
$\begin{split}E^{\Sigma}&:=\partial_{\chi_{2}}E_{0}|_{{\bm{\chi}}=0}=\frac{1}{\tau_{c}}\frac{\tilde{k}(T_{1}-T_{2})^{2}}{T_{1}T_{2}(\tilde{\gamma}_{1}+1)},\end{split}$
(S70)
$\begin{split}E^{\Sigma^{2}}&:=\partial_{\chi_{2}}^{2}E_{0}|_{{\bm{\chi}}=0}\\\
&=\frac{1}{\tau_{c}}\frac{2\tilde{k}(T_{1}-T_{2})^{2}[T_{1}^{2}\tilde{\gamma}_{1}+T_{2}^{2}\tilde{\gamma}_{1}+T_{1}T_{2}(\tilde{\gamma}_{1}^{2}+1)]}{T_{1}^{2}T_{2}^{2}(\tilde{\gamma}_{1}+1)^{3}},\end{split}$
(S71)

where we define the dimensionless factors $\tilde{k}:=k/k_{0}$ and
$\tilde{\gamma}_{1}:=\gamma_{1}/\gamma_{0}$.

The geometric connection vectors for $\braket{Q_{\rm curv}}$, $\braket{Q_{\rm
curv}^{2}}_{c}$ and $\braket{\Sigma_{\rm curv}}$ are given by

$\begin{split}A^{Q}_{\mu}d\Lambda_{\mu}&=-\frac{(T_{2}\tilde{\gamma}_{1}+T_{1})}{2\tilde{k}(\tilde{\gamma}_{1}+1)^{2}}d\tilde{k}+\frac{T_{1}\tilde{\gamma}_{1}+T_{2}}{2(\tilde{\gamma}_{1}+1)^{3}}d\tilde{\gamma}_{1},\end{split}$
(S72)

$\begin{split}A^{Q^{2}}_{\mu}d\Lambda_{\mu}&=\frac{\tilde{\gamma}_{1}(T_{2}\tilde{\gamma}_{1}+T_{1})[T_{1}(\tilde{\gamma}_{1}-2)+3T_{2}]}{\tilde{k}(\tilde{\gamma}_{1}+1)^{4}}d\tilde{k}\\\
&+\frac{3\tilde{\gamma}_{1}(T_{1}-T_{2})(T_{1}\tilde{\gamma}_{1}+T_{2})}{(\tilde{\gamma}_{1}+1)^{5}}d\tilde{\gamma_{1}},\end{split}$
(S73)

$\begin{split}A^{\Sigma}_{\mu}d\Lambda_{\mu}&=\frac{(T_{2}\tilde{\gamma}_{1}+T_{1})(T_{1}\tilde{\gamma}_{1}+T_{2})}{2\tilde{k}T_{1}T_{2}(\tilde{\gamma}_{1}+1)^{2}}d\tilde{k}\\\
&+\frac{(T_{1}-T_{2})(T_{1}\tilde{\gamma}_{1}+T_{2})}{2T_{1}T_{2}(\tilde{\gamma}_{1}+1)^{3}}d\tilde{\gamma}_{1},\end{split}$
(S74)

$\begin{split}A^{\Sigma^{2}}_{\mu}d\Lambda_{\mu}&=\frac{3(T_{1}-T_{2})^{2}(T_{1}\tilde{\gamma}_{1}+T_{2})}{T_{1}^{2}T_{2}^{2}(\tilde{\gamma}_{1}+1)^{4}}[\frac{\tilde{\gamma}_{1}(T_{2}\tilde{\gamma}_{1}+T_{1})}{\tilde{k}}d\tilde{k}\\\
&+\frac{(T_{1}-T_{2})\tilde{\gamma}_{1}}{\tilde{\gamma}_{1}+1}d\tilde{\gamma}_{1}].\end{split}$
(S75)

Obviously, from Eq. (S70), Eq. (S71), Eq. (S74) and Eq. (S75), we can easily
show that both the dynamic and adiabatic geometric components of the entropy
production (and its variance) vanish in the isothermal case where
$T_{1}(t)=T_{2}(t)$. Particularly, although
$A^{\Sigma}_{\mu\nu}d\Lambda_{\mu}\neq 0$ in this situation, it simplifies to
$\frac{1}{2}d\ln\tilde{k}$ as a total derivative of $\tilde{k}$ and therefore
the accumulation during one cycle
$\int_{\partial\Omega}A^{\Sigma}_{\mu\nu}d\Lambda_{\mu}$ is zero.

The geometric metrics are given by

$\begin{split}\mathfrak{g}^{Q}_{\mu\nu}d\Lambda_{\mu}d\Lambda_{\nu}&=\frac{\tau_{c}\tilde{\gamma}_{1}}{4\tilde{k}^{3}(\tilde{\gamma}_{1}+1)^{3}}[-(T_{2}\tilde{\gamma}_{1}+T_{1})(d\tilde{k})^{2}+\tilde{k}T_{2}d\tilde{k}d\tilde{\gamma}_{1}+\frac{\tilde{k}^{2}(T_{1}-T_{2})}{(\tilde{\gamma}_{1}+1)^{2}}(d\tilde{\gamma}_{1})^{2}],\end{split}$
(S76)

$\begin{split}\mathfrak{g}^{Q^{2}}_{\mu\nu}d\Lambda_{\mu}d\Lambda_{\nu}&=\tau_{c}[\frac{\tilde{\gamma}_{1}(T_{2}\tilde{\gamma}_{1}+T_{1})[6T_{2}\tilde{\gamma}_{1}+T_{1}(\tilde{\gamma}_{1}^{2}-4\tilde{\gamma}_{1}+1)]}{2\tilde{k}^{3}(\tilde{\gamma}_{1}+1)^{5}}(d\tilde{k})^{2}-\frac{T_{2}\tilde{\gamma}_{1}[6T_{2}\tilde{\gamma}_{1}+T_{1}(\tilde{\gamma}_{1}^{2}-4\tilde{\gamma}_{1}+1)]}{2\tilde{k}^{2}(\tilde{\gamma}_{1}+1)^{5}}d\tilde{k}d\tilde{\gamma}_{1}\\\
&-\frac{\tilde{\gamma}_{1}[T_{1}T_{2}(-3\tilde{\gamma}_{1}^{2}+14\tilde{\gamma}_{1}-3)+T_{1}^{2}(\tilde{\gamma}_{1}^{2}-8\tilde{\gamma}_{1}+1)+T_{2}^{2}(\tilde{\gamma}_{1}^{2}-8\tilde{\gamma}_{1}+1)]}{2\tilde{k}(\tilde{\gamma}_{1}+1)^{7}}(d\tilde{\gamma}_{1})^{2}],\end{split}$
(S77)

$\mathfrak{g}^{\Sigma}_{\mu\nu}d\Lambda_{\mu}d\Lambda_{\nu}=\tau_{c}[\frac{\tilde{\gamma}_{1}(T_{2}\tilde{\gamma}_{1}+T_{1})(T_{1}\tilde{\gamma}_{1}+T_{2})}{4\tilde{k}^{3}T_{1}T_{2}(\tilde{\gamma}_{1}+1)^{3}}(d\tilde{k})^{2}+\frac{(T_{1}^{2}-T_{2}^{2})\tilde{\gamma}_{1}}{4\tilde{k}^{2}T_{1}T_{2}(\tilde{\gamma}_{1}+1)^{3}}d\tilde{k}d\tilde{\gamma}_{1}+\frac{(T_{1}-T_{2})^{2}\tilde{\gamma}_{1}}{4\tilde{k}T_{1}T_{2}(\tilde{\gamma}_{1}+1)^{5}}(d\tilde{\gamma}_{1})^{2}],$
(S78)

$\begin{split}\mathfrak{g}^{\Sigma^{2}}_{\mu\nu}d\Lambda_{\mu}d\Lambda_{\nu}&=\tau_{c}[\frac{\tilde{\gamma}_{1}(T_{2}\tilde{\gamma}_{1}+T_{1})[6T_{1}^{3}\tilde{\gamma}_{1}^{2}+6T_{2}^{3}\tilde{\gamma}_{1}+T_{1}T_{2}^{2}(7\tilde{\gamma}_{1}^{2}-10\tilde{\gamma}_{1}+1)+T_{1}^{2}T_{2}\tilde{\gamma}_{1}(\tilde{\gamma}_{1}^{2}-10\tilde{\gamma}_{1}+7)]}{2\tilde{k}^{3}T_{1}^{2}T_{2}^{2}(\tilde{\gamma}_{1}+1)^{5}}(d\tilde{k})^{2}\\\
&+\frac{\tilde{\gamma}_{1}(T_{1}^{2}-T_{2}^{2})[6T_{1}^{2}\tilde{\gamma}_{1}+6T_{2}^{2}\tilde{\gamma}_{1}+T_{1}T_{2}(\tilde{\gamma}_{1}^{2}-10\tilde{\gamma}_{1}+1)]}{2\tilde{k}^{2}T_{1}^{2}T_{2}^{2}(\tilde{\gamma}_{1}+1)^{5}}d\tilde{k}d\tilde{\gamma}_{1}\\\
&+\frac{\tilde{\gamma}_{1}(T_{1}-T_{2})^{2}[T_{1}T_{2}(3\tilde{\gamma}_{1}^{2}-14\tilde{\gamma}_{1}+3)-T_{1}^{2}(\tilde{\gamma}_{1}^{2}-8\tilde{\gamma}_{1}+1)-T_{2}^{2}(\tilde{\gamma}_{1}^{2}-8\tilde{\gamma}_{1}+1)]}{2\tilde{k}T_{1}^{2}T_{2}^{2}(\tilde{\gamma}_{1}+1)^{7}}(d\tilde{\gamma}_{1})^{2}].\end{split}$
(S79)

In the non-biased situation. $T_{1}=T_{2}$ and the metric structure for
$\braket{\Sigma}$ and $\braket{\Sigma^{2}}_{c}$ simplifies to be

$\begin{split}\mathfrak{g}^{\Sigma}_{\mu\nu}d\Lambda_{\mu}d\Lambda_{\nu}&=\tau_{c}\frac{\tilde{\gamma}_{1}}{4\tilde{k}^{3}(\tilde{\gamma}_{1}+1)}(d\tilde{k})^{2},\\\
\mathfrak{g}^{\Sigma^{2}}_{\mu\nu}d\Lambda_{\mu}d\Lambda_{\nu}&=2(\mathfrak{g}_{\Sigma})_{\mu\nu}d\Lambda_{\mu}d\Lambda_{\nu}.\end{split}$
(S80)

In this Brownian model, as an illustration, the metric structure of the
average and variance of entropy production is of the form that our bound Eq.
(S53) can actually be saturated. In numerical simulation,s we randomly choose
$10^{1/2}\leq\tau_{p}\leq 100,T_{1}=1,0.1\leq T_{2}\leq 1.9$, the driving
amplitude $a=0.5$ and the phase difference $0\leq\phi\leq 2\pi$. The results
validate the Geometric TURs in Fig S3(c).

To consider the geometric bound on the variance of entropy production in Fig.
S3(d), we fix $T_{1}=T_{2}=1$ and select the driving amplitude $a\in[0.4,0.9]$
and the phase difference $0\leq\phi\leq 2\pi$.

As shown in Fig. S3(b), the nonadiabatic heat pump effect is bounded from
above by $\braket{Q}_{\rm curv}/\tau_{\rm p}$. We can easily show that the
average heat flux $\braket{Q}/\tau_{\rm p}$ reaches its maximum
$-\braket{Q_{\rm curv}}^{2}/(4\tau_{\rm p}\braket{Q_{\rm metr}})$ at the
optimal period $-2\tau_{\rm p}\braket{Q_{\rm metr}}/\braket{Q_{\rm curv}}$. In
Fig. S3(c), independent of the driving protocols, the entropy production
$\braket{\Sigma}$ is bounded from below by $\Sigma_{\rm g}$ (blue circles) as
derived in Eq. (S50), but breaks the corresponding steady state bound
$\braket{\Sigma}\geq\Sigma_{\rm b}:=2\braket{Q}^{2}/\braket{Q^{2}}_{c}$
(orange triangles). Furthermore, in this system, the metric expression for the
average of entropy production $\mathfrak{g}_{\mu\nu}^{\Sigma}$ and entropy
variance $\mathfrak{g}_{\mu\nu}^{\Sigma^{2}}$ for $T_{1}=T_{2}$, i.e. Eq.
(S79), satisfies

$\mathfrak{g}_{\mu\nu}^{\Sigma^{2}}=2\mathfrak{g}_{\mu\nu}^{\Sigma}=\frac{\gamma_{1}\gamma_{2}}{2k^{3}(\gamma_{1}+\gamma_{2})}\begin{pmatrix}1&0\\\
0&0\end{pmatrix},$ (S81)

implying that our geometric bound Eq. (S53) is actually saturable by
reparametrizing the protocol in terms of the thermodynamic length, i.e. the
time spend around a parameter point $dt=(\tau_{\rm
p}/{\mathcal{L}})\sqrt{\mathfrak{g}_{\mu\nu}^{\Sigma}d\Lambda_{\nu}d\Lambda_{\nu}}$
[6, 7], as shown by the green triangles in Fig. S3(d). Here,
$\mathfrak{g}_{\mu\nu}^{\Sigma}:=\partial_{\chi_{2}}\mathfrak{g}_{\mu\nu}|_{{\bm{\chi}}=0}$
and
$\mathfrak{g}_{\mu\nu}^{\Sigma^{2}}:=\partial_{\chi_{2}}^{2}\mathfrak{g}_{\mu\nu}|_{{\bm{\chi}}=0}$.
The blue circles in Fig. S3(d) shows the validity of Eq. (S53) in the constant
speed protocols.

## References

  * Gingrich  _et al._ [2017] T. R. Gingrich, G. M. Rotskoff, and J. M. Horowitz, Inferring dissipation from current fluctuations, J. Phys. A 50, 184004 (2017).
  * Dechant and Sasa [2020] A. Dechant and S.-i. Sasa, Fluctuation–response inequality out of equilibrium, Proc. Natl. Acad. Sci. U.S.A. 117, 6430 (2020).
  * Liu  _et al._ [2020] K. Liu, Z. Gong, and M. Ueda, Thermodynamic uncertainty relation for arbitrary initial states, Phys. Rev. Lett. 125, 140602 (2020).
  * Koyuk and Seifert [2020] T. Koyuk and U. Seifert, Thermodynamic uncertainty relation for time-dependent driving, Phys. Rev. Lett. 125, 260604 (2020).
  * Ren  _et al._ [2012] J. Ren, S. Liu, and B. Li, Geometric heat flux for classical thermal transport in interacting open systems, Phys. Rev. Lett. 108, 210603 (2012).
  * Crooks [2007] G. E. Crooks, Measuring thermodynamic length, Phys. Rev. Lett. 99, 100602 (2007).
  * Brandner and Saito [2020] K. Brandner and K. Saito, Thermodynamic geometry of microscopic heat engines, Phys. Rev. Lett. 124, 040602 (2020).