Multidimensional entropic bound: Estimator of entropy production for Langevin dynamics with an arbitrary time-dependent protocol

Sangyun Lee School of Physics, Korea Institute for Advanced Study, Seoul, 02455, Korea Dong-Kyum Kim Data Science Group, Institute for Basic Science, Daejeon, 34126, Korea Jong-Min Park School of Physics, Korea Institute for Advanced Study, Seoul, 02455, Korea Won Kyu Kim School of Computational Sciences, Korea Institute for Advanced Study, Seoul, 02455, Korea Hyunggyu Park School of Physics, Korea Institute for Advanced Study, Seoul, 02455, Korea Quantum Universe Center, Korea Institute for Advanced Study, Seoul 02455, Korea Jae Sung Lee [email protected] School of Physics, Korea Institute for Advanced Study, Seoul, 02455, Korea

Abstract

Entropy production (EP) is a key quantity in thermodynamics, and yet measuring EP has remained a challenging task. Here we introduce an EP estimator, called multidimensional entropic bound (MEB), utilizing an ensemble of trajectories. The MEB can accurately estimate the EP of overdamped Langevin systems with an arbitrary time-dependent protocol. Moreover, it provides a unified platform to accurately estimate the EP of underdamped Langevin systems under certain conditions. In addition, the MEB is computationally efficient because optimization is unnecessary. We apply our developed estimator to three physical systems driven by time-dependent protocols pertaining to experiments using optical tweezers: a dragged Brownian particle, a pulling process of a harmonic chain, and an unfolding process of an RNA hairpin. Numerical simulations confirm the validity and efficiency of our method.

I introduction

Entropy production (EP), referring to the quantification of the irreversibility of a thermodynamic process, is one of the most fundamental thermodynamic quantities. The EP was originally identified in the Clausius form in equilibrium thermodynamics. More recently, crucial progress in the field of thermodynamics has been the extension of the EP to general nonequilibrium phenomena at the level of a single stochastic trajectory. This extension triggered a renaissance of thermodynamics, namely the establishment of stochastic thermodynamics. Based on the novel EP formulation, EP theories have been developed and extensively studied over the last two decades. An early one is the fluctuation theorem Jarzynski (1997); Crooks (1999); Seifert (2019, 2005); Collin et al. (2005), which can be understood as a generalization of the thermodynamic second law. Later developments include a group of thermodynamic trade-off relations such as the thermodynamic uncertainty relation (TUR) Barato and Seifert (2015); Gingrich et al. (2016); Horowitz and Gingrich (2017); Dechant and Sasa (2018a); Hasegawa and Vu (2019); Dechant (2019); Hasegawa and Van Vu (2019); Koyuk and Seifert (2020), the power-efficiency trade-off relation Benenti et al. (2011); Brandner et al. (2013); Shiraishi et al. (2016); Dechant and Sasa (2018b); Lee et al. (2020), and the speed limit Shiraishi et al. (2018); Nicholson et al. (2020); Vo et al. (2020); Falasco and Esposito (2020); Yoshimura and Ito (2021); Lee et al. (2022).

Subsequently, experimentally feasible methods for measuring the EP have been actively suggested and discussed Roldán and Parrondo (2010, 2012); Battle et al. (2016); Avinery et al. (2019); Martiniani et al. (2019); Li et al. (2019); Manikandan et al. (2020); Martínez et al. (2019); Kim et al. (2020); Otsubo et al. (2020, 2022); Skinner and Dunkel (2021); Kim et al. (2022); van der Meer et al. (2022); Seif et al. (2021). In fact, measuring EP is not a trivial task. It is almost impossible to measure the EP by using its definition, the logarithmic ratio of forward and time-reversal path probabilities Seifert (2012), since all path probabilities cannot be measured, especially for a continuous system. Instead, there exist two direct EP measurement methods using the “equality” for the total EP, $\Delta S^{\rm tot}$ . The first method uses the equality $\Delta S^{\rm tot}=\Delta S^{\rm sys}+Q/T$ , where $\Delta S^{\rm sys}$ is the Shannon entropy change of a system and $Q$ is dissipated heat into a reservoir at temperature $T$ Sekimoto (2010); Seifert (2012). In experiments, it is difficult to measure the amount of heat flow accurately with a calorimeter. One may calculate $Q$ from trajectory data instead, but this requires full knowledge of the external and internal forces acting on the system Sekimoto (2010). Therefore, this method is not practically useful for complicated cases such as a biological system. The second direct method uses the equality for the average EP in terms of the probability density function (PDF) and the irreversible probability current as presented in Eq. (6) of Ref. Li et al. (2019). The PDF and the irreversible probability current can be estimated solely from system trajectories without knowledge of applied forces in the overdamped Langevin dynamics. Nevertheless, obtaining them precisely for a high-dimensional system is infeasible in practice, which is called the “curse of dimensionality”.

To overcome these shortcomings of the direct methods, several indirect methods using a thermodynamic “inequality” have been suggested. Here, the EP can be estimated from an ensemble of system trajectories, and the curse of dimensionality can be mitigated by measuring several observable currents only. The indirect methods are based on an inequality in the general form of $\Delta S^{\rm tot}\geq B(\Theta)$ , where the EP bound $B(\Theta)$ is determined by an observable current $\Theta$ . In a certain condition, one can find an optimal observable current $\Theta^{*}$ , which yields $B(\Theta^{*})=\Delta S^{\rm tot}$ . Then, the EP can be accurately estimated by measuring $\Theta^{*}$ .

Regarding the above indirect methods, there exist two representative inequalities. The first inequality is in TUR form Barato and Seifert (2015); Gingrich et al. (2016); Horowitz and Gingrich (2017); Dechant and Sasa (2018a); Hasegawa and Vu (2019), where the EP is bounded by the relative fluctuation of a certain observable current. To access a tighter bound of this TUR, multidimensional TUR Dechant (2019) and Monte Carlo methods Li et al. (2019) have been developed. However, TURs depend on the nature of the system dynamics; e.g, the TUR must be modified when a time-dependent protocol is involved Koyuk and Seifert (2020) or when a system follows underdamped Langevin dynamics Lee et al. (2019); Hasegawa and Van Vu (2019); Lee et al. (2021). Thus, EP estimation based on TURs is not universal. Moreover, if we use a TUR for an underdamped system, EP estimation is not possible from only system trajectories, but rather needs full knowledge and full controllability of the applied forces Lee et al. (2019); Hasegawa and Van Vu (2019); Lee et al. (2021). Thus, no proper method via TUR exists for estimating the EP solely from system trajectories for underdamped dynamics.

The second inequality for the indirect method is the Donsker–Varadhan inequality Donsker and Varadhan (2020). Recently, a machine learning technique named NEEP (neural estimator for entropy production) Kim et al. (2020); Otsubo et al. (2020); Kim et al. (2022) utilized this inequality as an optimized function for a given neural network. Though this technique yields a reliable result in overdamped Langevin systems, a high computational cost is required for a process with a time-dependent protocol since the parameters of the neural network should be reoptimized every single time. Otherwise, this machine learning technique has also been applied to underdamped Langevin dynamics; however, it has difficulty in estimating the EP accurately for large inertia Kim et al. (2022).

In this study, we propose a unified and computationally efficient method to estimate the EP by using the entropic bound (EB) inequality introduced by Dechant and Sasa Dechant and Sasa (2018b). Inspired by the multidimensional TUR Dechant (2019), we use multiple observable currents to obtain the optimal EB for the EP. Thus, we call this the “multidimensional entropic bound” (MEB). The MEB is universal in the sense that it provides a unified platform to estimate the EP for both overdamped and underdamped systems regardless of the time dependence of the driving protocol. The MEB can estimate the EP from system-trajectory information when an irreversible force is absent, a common experimental setup. When an irreversible force is involved, additional information about the force is required to estimate the EP. For an underdamped system, supplementary information on the velocity relaxation time, which can be determined experimentally, is necessary to estimate the EP.

This paper is organized as follows. In Sec. II, we derive the formulae for the MEB and describe the EP estimation process using it. In Sec. III, we explain the relation between the MEB and various TUR bounds. In Sec. IV, we apply the MEB to three systems with time-dependent driving forces that can be realized in experiments using optical tweezers. We conclude the paper in Sec. V.

II Multidimensional Entropic Bound

The EB is the inequality between the EP and an observable current Dechant and Sasa (2018b). As this bound holds for both overdamped and underdamped Langevin systems with an arbitrary time-dependent protocol, it can be a good starting point to obtain a unified and efficient EP estimator applicable to both overdamped and underdamped Langevin dynamics. In this section, we introduce the multidimensional entropic bound (MEB) estimator by incorporating multiple observable currents systematically into the EB estimator.

II.1 Derivation of the integral and the rate EB

Here, we consider an $M$ -dimensional Langevin system with a state vector ${\bm{q}}(t)=(q_{1},\cdots,q_{M})^{\textsf{T}}$ , where T denotes the transpose of a matrix, described by the following equation of motion:

\displaystyle\dot{\bm{q}}(t)=\bm{A}(\bm{q}(t),t)+\sqrt{2\mathsf{B}(\bm{q}(t),t)}\bullet\bm{\xi}(t),

(1)

where $\bm{A}=(A_{1},\cdots,A_{M})^{\textsf{T}}$ is a time-dependent drift force, $\mathsf{B}$ is a positive-definite symmetric $M\times M$ diffusion matrix, and $\bm{\xi}=(\xi_{1},\cdots,\xi_{M})^{\textsf{T}}$ is a Gaussian white noise satisfying $\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$ for $i,j\in\{1,\cdots,M\}$ . The symbol $\bullet$ in Eq. (1) represents the Itô product. From now on, we sometimes drop the arguments of functions for simplicity.

A component of $\bm{q}$ can be an odd-parity variable such as a velocity under time-reversal operation. The time reversal of a state $\bm{q}$ is denoted by $\tilde{\bm{q}}=(\tilde{q}_{1},\cdots,\tilde{q}_{M})^{\textsf{T}}$ with $\tilde{q}_{i}=\epsilon_{i}q_{i}$ , where $\epsilon_{i}=1$ for an even-parity variable and $\epsilon_{i}=-1$ otherwise. The drift force can be divided into reversible and irreversible parts as $\bm{A}(\bm{q},t)=\bm{A}^{\rm rev}(\bm{q},t)+\bm{A}^{\rm irr}(\bm{q},t)$ with Risken (1991)

	$\displaystyle\bm{A}^{\rm rev}(\bm{q},t)$	$\displaystyle\equiv\frac{1}{2}\left[\bm{A}(\bm{q},t)-\bm{\epsilon}\odot\bm{A}^{\dagger}(\bm{\epsilon}\odot\bm{q},t)\right],$
	$\displaystyle\bm{A}^{\rm irr}(\bm{q},t)$	$\displaystyle\equiv\frac{1}{2}\left[\bm{A}(\bm{q},t)+\bm{\epsilon}\odot\bm{A}^{\dagger}(\bm{\epsilon}\odot\bm{q},t)\right],$		(2)

where $\bm{\epsilon}=(\epsilon_{1},\cdots,\epsilon_{M})^{\textsf{T}}$ , $\odot$ denotes the element-wise product, i.e., $\bm{a}\odot\bm{b}=(\cdots,a_{i}b_{i},\cdots)^{\textsf{T}}$ , and $\dagger$ is an operation changing the sign of the odd-parity parameters.

The Fokker–Planck equation associated with Eq. (1) is

\displaystyle\partial_{t}P(\bm{q},t)=-\bm{\nabla}[\bm{J}^{\rm rev}(\bm{q},t)+\bm{J}^{\rm irr}(\bm{q},t)]~{},

(3)

with the PDF $P(\bm{q},t)$ . The reversible current $\bm{J}^{\rm rev}(\bm{q},t)$ and the irreversible current $\bm{J}^{\rm irr}(\bm{q},t)$ are defined as

	$\displaystyle J_{i}^{\rm rev}(\bm{q},t)$	$\displaystyle\equiv A_{i}^{\rm rev}(\bm{q},t)P(\bm{q},t),$		(4)
	$\displaystyle J_{i}^{\rm irr}(\bm{q},t)$	$\displaystyle\equiv A_{i}^{\rm irr}(\bm{q},t)P(\bm{q},t)-\sum_{j}\partial_{q_{j}}\left[\mathsf{B}_{ij}(\bm{q},t)P(\bm{q},t)\right],$		(5)

with $\mathsf{B}(\bm{\epsilon}\odot\bm{q},t)=\mathsf{B}(\bm{q},t)$ . Note that for an overdamped Langevin system with even-parity variables only, $\bm{A}^{\rm rev}(\bm{q},t)$ vanishes, and thus $\bm{J}^{\rm rev}(\bm{q},t)=0$ and the total current coincides with $\bm{J}^{\rm irr}(\bm{q},t)$ . As dissipation originates from the irreversible current, the EP is determined only by $\bm{J}^{\rm irr}(\bm{q},t)$ . Therefore, the total EP rate $\sigma^{\rm tot}$ is given by Spinney and Ford (2012); Kwon et al. (2016); Dechant and Sasa (2018b)

\displaystyle\sigma^{\rm tot}(t)\equiv\int d\bm{q}\frac{{\bm{J}}^{\rm irr}(\bm{q},t)^{\textsf{T}}\mathsf{B}(\bm{q},t)^{-1}{\bm{J}}^{\rm irr}(\bm{q},t)}{P(\bm{q},t)}.

(6)

Hereafter, we use the $k_{B}=1$ unit. Note that for underdamped Langevin systems, the matrix is not directly invertible since $\mathsf{B}_{ij}=0$ when the component index $i$ or $j$ denotes a positional variable. For such an index $i$ , we first set $\mathsf{B}_{ii}=b$ ( $b>0$ ) and $\mathsf{B}_{ij}=0$ ( $i\neq j$ ), then take the inverse of $\mathsf{B}$ and calculate ${\bm{J}}^{\rm irr}(\bm{q},t)^{\textsf{T}}\mathsf{B}(\bm{q},t)^{-1}{\bm{J}}^{\rm irr}(\bm{q},t)$ in Eq. (6), and finally take the $b\rightarrow 0$ limit. Since $J_{i}^{\rm irr}(\bm{q},t)\propto b$ , this limit leads to ${J}_{i}^{\rm irr}(\bm{q},t)^{2}\mathsf{B}_{ii}(\bm{q},t)^{-1}\sim b\rightarrow 0$ . For underdamped systems, this procedure amounts to writing B and ${\bm{J}}^{\rm irr}$ in terms of velocity-variable components only.

In this study, we consider the following form of an averaged observable current generated by the irreversible current during time $\tau$ :

\displaystyle\langle\Theta(\tau)\rangle=\int^{\tau}_{0}dt\int d\bm{q}~{}\bm{\Lambda}(\bm{q},t)^{\textsf{T}}\bm{J}^{\rm irr}(\bm{q},t),

(7)

where $\bm{\Lambda}(\bm{q},t)=(\Lambda_{1},\cdots,\Lambda_{M})^{\textsf{T}}$ is a weight vector of the irreversible current for a given observable. Then, the averaged current rate at time $t$ is given as

\displaystyle\langle\dot{\Theta}(t)\rangle=\int d\bm{q}~{}\bm{\Lambda}(\bm{q},t)^{\textsf{T}}\bm{J}^{\rm irr}(\bm{q},t)~{}.

(8)

The EB in an integral form can be derived from Eq. (7) as follows:

	$\displaystyle\langle\Theta(\tau)\rangle$
	$\displaystyle=\int^{\tau}_{0}dt\int d\bm{q}P(\bm{q},t)^{\frac{1}{2}}\bm{\Lambda}(\bm{q},t)^{\textsf{T}}\mathsf{B}(\bm{q},t)^{\frac{1}{2}}\frac{\mathsf{B}(\bm{q},t)^{-\frac{1}{2}}\bm{J}^{\rm irr}(\bm{q},t)}{P(\bm{q},t)^{\frac{1}{2}}}$
	$\displaystyle\leq\sqrt{\int^{\tau}_{0}dt\left\langle\bm{\Lambda}^{\textsf{T}}\mathsf{B}\bm{\Lambda}\right\rangle_{\bm{q}}}\sqrt{\Delta S^{\rm tot}(\tau)},$		(9)

where $\langle\cdots\rangle_{\bm{q}}=\int d\bm{q}\cdots P(\bm{q},t)$ and the total EP $\Delta S^{\rm tot}(\tau)=\int_{0}^{\tau}dt~{}\sigma^{\rm tot}(t)$ . The Cauchy–Schwartz inequality is used for the last inequality of Eq. (9). Hence, the total EP is bounded in an integral form as

	$\displaystyle\Delta S^{\rm tot}(\tau)\geq\frac{\langle\Theta(\tau)\rangle^{2}}{\int^{\tau}_{0}ds\left\langle\bm{\Lambda}^{\textsf{T}}\mathsf{B}\bm{\Lambda}\right\rangle_{\bm{q}}}\equiv$	$\displaystyle\Delta S^{\rm EB}(\Theta,\tau).$
		(integral EB)		(10)

Similarly, the EB in a rate form can also be obtained from Eq. (8) as

\displaystyle\sigma^{\rm tot}(t)\geq\frac{\langle\dot{\Theta}(t)\rangle^{2}}{\left\langle\bm{\Lambda}^{\textsf{T}}\mathsf{B}\bm{\Lambda}\right\rangle_{\bm{q}}}\equiv\sigma^{\rm EB}(\Theta,t).~{}~{}\textrm{(rate EB)}

(11)

The equality of the integral EB is satisfied when the weight vector has the following form:

\displaystyle\bm{\Lambda}^{\rm e}(\bm{q},t)=c\frac{\mathsf{B}(\bm{q},t)^{-1}\bm{J}^{\rm irr}(\bm{q},t)}{P(\bm{q},t)}~{}~{}\textrm{(for the integral EB)},

(12)

where $c$ is an arbitrary constant that is independent of $\bm{q}$ and $t$ . This can be easily checked by inserting Eq. (12) into Eq. (10). The weight vector in this case corresponds to the observable current proportional to the total EP, i.e., $\langle\Theta(\tau)\rangle=c\Delta S^{\rm tot}(\tau)$ . Similarly, we find the equality condition for the rate EB as

\displaystyle\bm{\Lambda}^{\rm e}(\bm{q},t)=c(t)\frac{\mathsf{B}(\bm{q},t)^{-1}\bm{J}^{\rm irr}(\bm{q},t)}{P(\bm{q},t)}~{}~{}\textrm{(for the rate EB)},

(13)

where $c(t)$ is an arbitrary time-dependent function that is independent of $\bm{q}$ . This weight vector corresponds to the observable current rate as $\langle\dot{\Theta}(t)\rangle=c(t)\sigma^{\rm tot}(t)$ . Note that $c$ and $c(t)$ can be arbitrary, and thus we may choose $c$ and $c(t)$ freely in order to simplify the measurement of an observable current. A relevant example is presented in Sec. IV.1.

II.2 Derivation of the integral and the rate MEB

With the knowledge of the functional form of $\bm{\Lambda}^{\rm e}(\bm{q},t)$ , one may obtain the tight EP bound. However, except for very simple examples, it is impossible to identify $\bm{\Lambda}^{\rm e}(\bm{q},t)$ without knowing all driving and interaction forces. Instead, we measure multiple observable currents to access a tighter bound, thereby systematically approaching the total EP. Our MEB method is analogous to the multidimensional TUR Dechant (2019) but is more general in the sense that it can be applicable to wider classes of Langevin dynamics.

In this method, a linear combination of multiple weight vectors is adopted to approximate $\bm{\Lambda}^{\rm e}(\bm{q},t)$ . The linear combination of $\ell$ weight vectors $\{\Lambda_{i,1},\cdots,\Lambda_{i,\ell}\}$ for the $i$ -th component is written as

\displaystyle\Lambda_{i}^{(\ell)}(\bm{q},t)=\sum_{\alpha=1}^{\ell}k_{i,\alpha}\Lambda_{i,\alpha}(\bm{q},t),

(14)

where $k_{i,\alpha}$ is the coefficient for $\Lambda_{i,\alpha}(\bm{q},t)$ and is independent of $\bm{q}$ and $t$ . From now on, we will consider the case $\mathsf{B}_{ij}=B_{i}\delta_{ij}$ for simplicity. Even when the diffusion matrix has off-diagonal elements, we can always diagonalize the diffusion matrix by using a proper transformation of the coordinate if the full information of $\mathsf{B}_{ij}$ is given. Thus, we can still set $\mathsf{B}_{ij}=B_{i}\delta_{ij}$ on the transformed coordinate. In cases where it is difficult to obtain the information of $\mathsf{B}_{ij}$ , and thus not possible to find the proper transformation, we cannot use the following MEB in component-wise form. However, even in such cases, we can still derive the MEB in “component-combined” form as we show in Appendix C. The observable in Eq. (7) can be divided into the sum of its components as $\langle\Theta(\tau)\rangle=\sum_{i=1}^{M}\langle\Theta_{i}(\tau)\rangle$ , where $\langle\Theta_{i}(\tau)\rangle=\int^{\tau}_{0}dt\int d\bm{q}~{}\Lambda_{i}(\bm{q},t)J_{i}^{\rm irr}(\bm{q},t)$ . With the $i$ -th component current $\langle\Theta_{i}(\tau)\rangle$ , we derive the component-wise EB as

\displaystyle\Delta S_{i}(\tau)\geq\frac{\langle\Theta_{i}(\tau)\rangle^{2}}{\int^{\tau}_{0}dt\left\langle\Lambda_{i}B_{i}\Lambda_{i}\right\rangle_{\bm{q}}}.~{}~{}\textrm{($i$-th integral EB)},

(15)

where $\Delta S_{i}(\tau)=\int_{0}^{\tau}dt\sigma_{i}(t)$ with the $i$ -th component EP rate $\sigma_{i}(t)=\int d\bm{q}~{}B_{i}(\bm{q},t)^{-1}J_{i}^{\rm irr}(\bm{q},t)^{2}/P(\bm{q},t)$ . Thus, $\Delta S^{\rm tot}(\tau)=\sum_{i}\Delta S_{i}(\tau)$ and $\sigma^{\rm tot}(t)=\sum_{i}\sigma_{i}(t)$ . By substituting Eq. (14) into Eq. (15), we have

\displaystyle\Delta S_{i}(\tau)\geq\frac{\left\{\bm{k}_{i}^{\textsf{T}}{\langle\bm{\Theta}}_{i}^{(\ell)}(\tau)\rangle\right\}^{2}}{\bm{k}_{i}^{\textsf{T}}\mathsf{L}_{i}^{(\ell)}(\tau)\bm{k}_{i}}\equiv\Delta\hat{S}_{i}^{(\ell)}(\bm{k}_{i}),

(16)

where $\bm{k}_{i}=(k_{i,1},\cdots,k_{i,\ell})^{\textsf{T}}$ , and the vector $\langle\bm{\Theta}_{i}^{(\ell)}(\tau)\rangle=(\langle\Theta_{i,1}(\tau)\rangle,\cdots,\langle\Theta_{i,\ell}(\tau)\rangle)^{\textsf{T}}$ and the $\ell\times\ell$ matrix $\mathsf{L}_{i}(\tau)$ are defined as

	$\displaystyle\langle\Theta_{i,\alpha}(\tau)\rangle\equiv$	$\displaystyle\int^{\tau}_{0}dt\int d\bm{q}~{}\Lambda_{i,\alpha}(\bm{q},t)J_{i}^{\rm irr}(\bm{q},t)~{}~{}\textrm{and}$		(17)
	$\displaystyle\left(\mathsf{L}_{i}^{(\ell)}(\tau)\right)_{\alpha,\beta}\equiv$	$\displaystyle\int^{\tau}_{0}dt\left\langle\Lambda_{i,\alpha}(\bm{q},t)B_{i}(\bm{q},t)\Lambda_{i,\beta}(\bm{q},t)\right\rangle_{\bm{q}},$		(18)

respectively. Note that $\mathsf{L}_{i}^{(\ell)}(t)$ is a positive definite matrix since $\bm{z}^{\textsf{T}}\mathsf{L}_{i}^{(\ell)}(\tau)\bm{z}=\int^{\tau}_{0}dt\int d\bm{q}\langle\left(\sum_{\alpha}\sqrt{B_{i}}\Lambda_{i,\alpha}z_{\alpha}\right)^{2}\rangle_{\bm{q}}>0$ for an arbitrary $\bm{z}$ .

The bound $\Delta\hat{S}_{i}^{(\ell)}(\bm{k}_{i})$ in Eq. (16) is a function of $\bm{k}_{i}$ ; thus, the tightest bound can be written as $\Delta\hat{S}_{i}^{(\ell)}(\bm{k}_{i}^{*})$ , where $\bm{k}_{i}^{*}$ is the optimal vector maximizing the bound. The optimal vector is obtained by solving the following equation:

	$\displaystyle\partial_{k_{i,\alpha}}\Delta\hat{S}_{i}^{(\ell)}(\bm{k}_{i})=0$
	$\displaystyle=\frac{2\bm{k}_{i}^{\textsf{T}}\langle\bm{\Theta}_{i}^{(\ell)}\rangle\left\{\bm{k}_{i}^{\textsf{T}}{\mathsf{L}}_{i}^{(\ell)}\bm{k}_{i}\cdot\langle\Theta_{i,\alpha}\rangle-\bm{k}_{i}^{\textsf{T}}\langle\bm{\Theta}_{i}^{(\ell)}\rangle\cdot({\mathsf{L}}_{i}^{(\ell)}\bm{k}_{i})_{\alpha}\right\}}{(\bm{k}_{i}^{\textsf{T}}\mathsf{L}_{i}^{(\ell)}\bm{k}_{i})^{2}}.$		(19)

We can easily check that the numerator vanishes with the choice of $\bm{k}_{i}^{*}=\mathsf{(}L_{i}^{(\ell)})^{-1}\cdot\langle\bm{\Theta}^{(\ell)}_{i}(\tau)\rangle$ . By plugging $\bm{k}_{i}^{*}$ into Eq. (16), we find the component-wise MEB as

\displaystyle\Delta S_{i}(\tau)\geq\langle\bm{\Theta}_{i}^{(\ell)}(\tau)\rangle^{\textsf{T}}\mathsf{(}L_{i}^{(\ell)}(\tau))^{-1}\langle\bm{\Theta}_{i}^{(\ell)}(\tau)\rangle\equiv\Delta S_{i}^{{\rm MEB}(\ell)}(\tau).

(20)

By summing over all components, we finally obtain our main result, namely MEB in integral form, as follows:

	$\displaystyle\Delta S^{\rm tot}(\tau)\geq$	$\displaystyle\sum_{i=1}^{M}\Delta S_{i}^{{\rm MEB}(\ell)}(\tau)$
		$\displaystyle\equiv\Delta S^{{\rm MEB}(\ell)}(\tau).~{}~{}~{}\textrm{(integral MEB)}$		(21)

We can also derive the MEB in rate form. The derivation of the rate MEB is essentially the same as that of the integral MEB. It starts from the component-wise rate EB as

\displaystyle\sigma_{i}(t)\geq\frac{\langle\dot{\Theta}_{i}(t)\rangle^{2}}{\left\langle\Lambda_{i}^{\textsf{T}}\mathsf{B}_{i}\Lambda_{i}\right\rangle_{\bm{q}}}~{}.~{}~{}~{}~{}~{}\textrm{($i$-th rate EB)}

(22)

In this case, it is usually sufficient to choose a time-independent basis as

\displaystyle\Lambda_{i}^{(\ell)}(\bm{q},t)=\sum_{\alpha=1}^{\ell}k_{i,\alpha}(t)\Lambda_{i,\alpha}(\bm{q}),

(23)

where the time-dependence is encoded in the coefficients instead of in $\Lambda_{i,\alpha}(\bm{q})$ , as the equality condition in Eq. (13) also allows a time-dependent overall multiplicative coefficient. Following the same derivation procedure as in Eqs. (16) $-$ (21), we finally obtain

$\displaystyle\sigma^{\rm tot}(t)\geq$	$\displaystyle\sum_{i=1}^{M}\langle\dot{\bm{\Theta}}_{i}^{(\ell)}(t)\rangle^{\textsf{T}}\left(\dot{\mathsf{L}}_{i}^{(\ell)}(t)\right)^{-1}\langle\dot{\bm{\Theta}}_{i}^{(\ell)}(t)\rangle$
$\displaystyle=$	$\displaystyle\sum_{i=1}^{M}\sigma_{i}^{{\rm MEB}(\ell)}(t)$
$\displaystyle\equiv$	$\displaystyle\sigma^{{\rm MEB}(\ell)}(t),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{(rate MEB)}$	(24)

where the vector $\langle\dot{\bm{\Theta}}_{i}^{(\ell)}(t)\rangle=(\langle\dot{\Theta}_{i,1}(t)\rangle,\cdots,\langle\dot{\Theta}_{i,\ell}(t)\rangle)^{\textsf{T}}$ and the $\ell\times\ell$ matrix $\dot{\mathsf{L}}_{i}^{(\ell)}(\tau)$ are defined as

	$\displaystyle\langle\dot{\Theta}_{i,\alpha}(t)\rangle\equiv$	$\displaystyle\int d\bm{q}~{}\Lambda_{i,\alpha}(\bm{q})J_{i}^{\rm irr}(\bm{q},t)~{}~{}\textrm{and}$		(25)
	$\displaystyle\left(\dot{\mathsf{L}}_{i}^{(\ell)}(t)\right)_{\alpha,\beta}\equiv$	$\displaystyle\left\langle\Lambda_{i,\alpha}(\bm{q})B_{i}(\bm{q},t)\Lambda_{i,\beta}(\bm{q})\right\rangle_{\bm{q}}.$		(26)

The total EP during a finite time $\tau$ can be evaluated by integrating $\sigma^{{\rm MEB}(\ell)}(t)$ over time.

The weight vectors for the rate MEB are not time-dependent, so we need a lower number of weight vectors to approximate $\Lambda_{i}^{\rm e}(\bm{q},t)$ compared to the integral MEB where time-dependent weight vectors are necessary. Practically, too many weight vectors can overfit all the fluctuations originating from a finite number of trajectories, sometimes giving rise to an undesirable result. Thus, the rate MEB is usually preferable in estimating the EP for a system driven by a time-dependent protocol.

The MEBs in Eqs. (21) and (24) are the maximum bounds for a given finite number of observables. If we add one more observable to the existing $\ell$ observables, the MEB becomes tighter, i.e.,

	$\displaystyle\Delta S_{i}^{{\rm MEB}(\ell+1)}(\tau)-\Delta S_{i}^{{\rm MEB}(\ell)}(\tau)$	$\displaystyle\geq 0,$
	$\displaystyle\sigma_{i}^{{\rm MEB}(\ell+1)}(t)-\sigma_{i}^{{\rm MEB}(\ell)}(t)$	$\displaystyle\geq 0.$		(27)

The proof of Eq. (27) is basically the same as that presented in Ref. Dechant and Sasa (2021). To be self-contained, we include the proof in Appendix A. As we increase $\ell$ , the MEB estimator also increases and eventually saturates to the maximum value, i.e., $\Delta S_{i}(\tau)$ or $\sigma_{i}(t)$ . It can often saturate even at finite $\ell=\ell^{\rm sat}$ for simple systems. Therefore, by observing the saturation regime in a plot of the EP estimator versus $\ell$ , we can accurately estimate the total EP (see Sec. IV) without resorting to a time-consuming optimization procedure.

There exists no limitation for choosing a set of $\ell$ weight vectors as long as they are linearly independent of each other. In this study, we adopt a Gaussian weight vector set Van Vu et al. (2020) for numerical verification of the rate MEB method in Sec. IV. The first weight vector is a Gaussian function, the width of which corresponds to the difference between the maximum and minimum state values. The second and third weight vectors are Gaussian functions with a width half that of the first one, and so on. This represents one way to add the Gaussian basis uniformly for a given state-variable range, which yields an accurate and reliable EP estimation as shown in Sec. IV. The mathematical expression for the weight vector set is as follows:

\displaystyle\{\Lambda_{i,\alpha}\}=\left\{\cdots,\exp{[-(q_{i}-a_{i,\alpha})^{2}/2b^{2}_{i,\alpha}]},\cdots\right\}.

(28)

In Eq. (28), $a_{i,\alpha}$ and $b_{i,\alpha}$ are given as

	$\displaystyle\{a_{i,\alpha}\}_{\alpha\leq\ell}=$	$\displaystyle\{\frac{1}{2}q_{i}^{\rm min}+\frac{1}{2}q_{i}^{\rm max},\frac{2}{3}q_{i}^{\rm min}+\frac{1}{3}q_{i}^{\rm max},$
		$\displaystyle\,\,\,\frac{1}{3}q_{i}^{\rm min}+\frac{2}{3}q_{i}^{\rm max},\frac{3}{4}q_{i}^{\rm min}+\frac{1}{4}q_{i}^{\rm max},...\},$
	$\displaystyle\{b_{i,\alpha}\}_{\alpha\leq\ell}=$	$\displaystyle\{\Delta q_{i},\frac{1}{2}\Delta q_{i},\frac{1}{2}\Delta q_{i},\frac{1}{3}\Delta q_{i},\frac{1}{3}\Delta q_{i},\frac{1}{3}\Delta q_{i},...\},$

where $q_{i}^{\rm max}$ ( $q_{i}^{\rm min}$ ) is the maximum (minimum) value of $q_{i}$ in a given trajectory ensemble and $\Delta q_{i}=q_{i}^{\rm max}-q_{i}^{\rm min}$ . Explicit forms of $a_{i,\alpha}$ and $b_{i,\alpha}$ are given as $a_{i,\alpha}=q_{i}^{\rm min}+(u(\alpha)+1)^{-1}[\alpha-\frac{u(\alpha)(u(\alpha)-1)}{2}]\Delta q_{i}$ and $b_{i,\alpha}=\Delta q_{i}/u(\alpha)$ for $\alpha\geq 1$ with $u(n)=\lfloor\frac{1+\sqrt{8\alpha-7}}{2}\rfloor$ . Here, $\lfloor x\rfloor=\tilde{m}$ , where $\tilde{m}$ is an integer satisfying $\tilde{m}\leq x<\tilde{m}+1$ . We note that the Gaussian weight vector set usually yields reliable estimation results compared to other sets. For example, the polynomial basis set $\{q_{i},q_{i}^{2},\cdots,q_{i}^{\ell}\}$ typically overestimates the EP with large fluctuations when $\ell$ is large, since the polynomial term with a large exponent is highly sensitive to rare-event data.

II.3 EP estimation procedure via MEB

In this section, we describe how to estimate the EP with the MEB from an ensemble of system trajectories. Here we consider both overdamped and underdamped Langevin dynamics. For an overdamped system, the system states consist of only position variables, i.e., $\bm{q}(t)=\bm{x}(t)=(x_{1},\cdots,x_{M})$ , and the Langevin equation is written as

\displaystyle\dot{x}_{i}(t)=\frac{1}{\gamma}F_{i}(\bm{q}(t),t)+\sqrt{2B_{i}(\bm{q},t)}\bullet\xi_{i}(t),

(29)

where $F_{i}$ is a force applied to the $x_{i}$ component. The reversible current $J_{i}^{\rm rev}=0$ , while the irreversible current is given as

\displaystyle J_{i}^{\rm irr}(\bm{q},t)=\left(\frac{1}{\gamma}F_{i}(\bm{q},t)-\partial_{x_{i}}B_{i}(\bm{q},t)\right)P(\bm{q},t).

(30)

In the case of underdamped dynamics, the system states consist of both position and velocity variables, i.e., $\bm{q}(t)=(x_{1},\cdots,x_{N},v_{1},\cdots,v_{N})$ with $M=2N$ , and the Langevin equation is written as

	$\displaystyle\dot{x}_{i}=v_{i}$
	$\displaystyle\dot{v}_{i}(t)=\frac{1}{m}F_{i}(\bm{q}(t),t)-\frac{\gamma}{m}v+\sqrt{2B_{i}(\bm{q},t)}\bullet\xi_{i}(t).$		(31)

As done in Eq. (2), the external force $F_{i}$ can be divided into reversible and irreversible parts as $F_{i}=F_{i}^{\rm rev}+F_{i}^{\rm irr}$ with

	$\displaystyle F_{i}^{\rm rev}(\bm{q},t)$	$\displaystyle\equiv\frac{1}{2}\left[F_{i}(\bm{q},t)+F_{i}^{\dagger}(\bm{\epsilon}\odot\bm{q},t)\right],$
	$\displaystyle F_{i}^{\rm irr}(\bm{q},t)$	$\displaystyle\equiv\frac{1}{2}\left[F_{i}(\bm{q},t)-F_{i}^{\dagger}(\bm{\epsilon}\odot\bm{q},t)\right].$		(32)

Then, the irreversible currents are given as

	$\displaystyle J_{x_{i}}^{\rm irr}(\bm{q},t)$	$\displaystyle=0,$
	$\displaystyle J_{v_{i}}^{\rm irr}(\bm{q},t)$	$\displaystyle=\left(\frac{1}{m}F_{i}^{\rm irr}(\bm{q},t)-\frac{\gamma}{m}v_{i}-\partial_{v_{i}}B_{i}(\bm{q},t)\right)P(\bm{q},t)~{},$		(33)

while the reversible currents $J_{x_{i}}^{\rm rev}(\bm{q},t)=v_{i}P(\bm{q},t)$ and $J_{v_{i}}^{\rm rev}(\bm{q},t)=\frac{1}{m}F_{i}^{\rm rev}(\bm{q},t)P(\bm{q},t)$ .

We focus on the rate MEB in the following discussions, but note that the procedure for the integral MEB is essentially the same.

II.3.1 Determination of $B_{i}$ and $\textsf{L}_{i}$

For an overdamped dynamics described by Eq. (29), the diffusivity $B_{i}$ is determined from the average of short-time mean square displacements as

\displaystyle B_{i}(\bm{q},t)=\lim_{\delta t\rightarrow 0}\frac{\langle\delta x_{i}(t)^{2}\rangle_{(\bm{q},t)}}{\delta t},~{}~{}\textrm{(overdamped)}

(34)

where $\delta x_{i}(t)\equiv x_{i}(t+\delta t)-x_{i}(t)$ and $\langle\cdots\rangle_{(\bm{q},t)}$ denotes the average over the trajectory ensemble at position $\bm{q}$ and time $t$ . When $B_{i}$ is independent of position and time, all short-time trajectories can be utilized for estimating the diffusivity. For an underdamped dynamics, $B_{i}$ can be estimated from the ensemble of velocity trajectories as

\displaystyle B_{i}(\bm{q},t)=\lim_{\delta t\rightarrow 0}\frac{\langle\delta v_{i}(t)^{2}\rangle_{(\bm{q},t)}}{\delta t},~{}~{}\textrm{(underdamped)}

(35)

where $\delta v_{i}(t)\equiv v_{i}(t+\delta t)-v_{i}(t)$ . With these estimations for $B_{i}$ , we calculate $(\dot{\textsf{L}}_{i}^{(\ell)}(t))_{\alpha,\beta}$ from Eq. (26). Note that $\lim_{\delta t\rightarrow 0}\frac{\langle\delta x_{i}(t)^{2}\rangle_{(\bm{q},t)}}{\delta t}=0$ in the underdamped case. In the case where the estimation of $B_{i}(\bm{q},t)$ requires heavy computation, instead of evaluating the diffusivity, we can directly estimate the $\left(\dot{\mathsf{L}}_{i}^{(\ell)}(t)\right)_{\alpha,\beta}$ matrix from the average of the products of two infinitesimal observable changes as

	$\displaystyle\left(\dot{\mathsf{L}}_{i}^{(\ell)}(t)\right)_{\alpha,\beta}$	$\displaystyle=\left\langle\Lambda_{i,\alpha}(\bm{q})B_{i}(\bm{q},t)\Lambda_{i,\beta}(\bm{q})\right\rangle_{\bm{q}}$
		$\displaystyle\approx\langle\Lambda_{i,\alpha}(\bm{q})\circ\delta q_{i}\cdot\Lambda_{i,\beta}(\bm{q})\circ\delta q_{i}\rangle_{\bm{q}}/(2\delta t).$		(36)

When the diffusion matrix has non-zero off-diagonal terms, we first have to directly evaluate ${\mathsf{B}}_{ij}=\lim_{\delta t\rightarrow 0}\langle\delta q_{i}\delta q_{j}\rangle_{(\bm{q},t)}/\delta t$ and find the proper transformation to diagonalize the matrix.

II.3.2 Measurement of an observable current

Now we describe how to measure the observable current rate in Eq. (25) for both overdamped and underdamped dynamics. First, in the overdamped dynamics with $\bm{q}(t)=(x_{1},\cdots,x_{M})$ , the $i$ -th component of an observable current rate can be measured by averaging $\Lambda_{i,\alpha}(\bm{q}(t))\circ\dot{x}_{i}(t)$ over the ensemble of system trajectories as

\displaystyle\langle\dot{\Theta}_{i,\alpha}(t)\rangle=\langle\Lambda_{i,\alpha}(\bm{q}(t))\circ\dot{x}_{i}(t)\rangle_{(\bm{q},t)},

(37)

where $\circ$ denotes the Stratonovich product. This can be checked from the fact that $\langle H(\bm{q},t)\circ\dot{x}_{i}\rangle_{(\bm{q},t)}=\int d\bm{q}~{}H(\bm{q},t)J_{i}^{\rm irr}(\bm{q},t)$ with an arbitrary function $H(\bm{q},t)$ . For $\Lambda_{i,\alpha}(\bm{q})=1$ , the observable current is the displacement in the $x_{i}$ direction as $\Theta_{i,\alpha}(\tau)=x_{i}(\tau)-x_{i}(0)$ .

In the underdamped dynamics with $\bm{q}(t)=(x_{1},\cdots,x_{N},v_{1},\cdots,v_{N})$ , by plugging Eq. (33) into Eq. (25), we have

	$\displaystyle\langle\dot{\Theta}_{i,\alpha}(t)\rangle=$	$\displaystyle-\frac{\gamma}{m}\langle\Lambda_{i,\alpha}(\bm{q})v_{i}\rangle_{\bm{q}}+\frac{1}{m}\langle\Lambda_{i,\alpha}(\bm{q})F_{i}^{\rm irr}(\bm{q},t)\rangle_{\bm{q}}$
		$\displaystyle+B_{i}(\bm{q},t)\langle\partial_{v_{i}}\Lambda_{i,\alpha}(\bm{q})\rangle_{\bm{q}}~{}.$		(38)

When $F_{i}^{\rm irr}=0$ and $\Lambda_{i,\alpha}$ has no explicit velocity dependence, the observable current is proportional to $\langle\Lambda_{i,\alpha}(\bm{q})v_{i}\rangle_{\bm{q}}$ , similar to that of the overdamped system, Eq. (37). However, extra information $\gamma/m$ is required to evaluate Eq. (38) in the underdamped case. This constant can be determined experimentally by measuring the velocity relaxation time Lee and Park (2018). Otherwise, when $F_{i}^{\rm irr}=0$ and $\Lambda_{i,\alpha}$ has an explicit velocity dependence, extra calculation of the last term in Eq. (38) is necessary. For the most general case with $F_{i}^{\rm irr}\neq 0$ (velocity-dependent force), $\langle\dot{\Theta}_{i,\alpha}(t)\rangle$ cannot be determined solely by system trajectories, but rather concrete information on the force is necessary.

II.3.3 EP Estimation

Utilizing numerical data for $\dot{\textsf{L}}_{i}^{(\ell)}$ and $\langle\dot{\Theta}_{i,\alpha}(t)\rangle$ obtained in Secs. II.3.1 and II.3.2, one can evaluate $\sigma^{{\rm MEB}(\ell)}(t)$ from Eq. (24) and then obtain $\Sigma^{{\rm MEB}(\ell)}(\tau)\equiv\int_{0}^{\tau}dt~{}\sigma^{{\rm MEB}(\ell)}(t)$ for each $\ell=1,2,3,\cdots$ . As proved in Sec. II.2, $\Sigma^{{\rm MEB}(\ell)}(\tau)$ is an increasing function of $\ell$ and saturates to the maximum value at some $\ell^{\rm sat}$ . This saturation indicates that $\Lambda_{i}^{\ell^{\rm sat}}$ coincides with $\Lambda_{i}^{\rm e}(\bm{q},t)$ , thus satisfying the equality of the EB. Therefore, the total EP corresponds to the MEB estimator at $\ell=\ell^{\rm sat}$ , i.e., $\Delta S^{\rm tot}(\tau)=\Sigma^{{\rm MEB}(\ell^{\rm sat})}(\tau)$ .

III Relation between MEB and TUR

In this section, we discuss the relation between MEB and TURs. We first consider a one-dimensional (1D) overdamped Langevin dynamics in the steady state (without a time-dependent protocol) as described by Eq. (29). The total EP $\Delta S^{\rm tot}(t,t^{\prime})$ during a small time segment between $t$ and $t^{\prime}=t+\delta t$ and the corresponding accumulated current $\Theta(t,t^{\prime})$ are written as

	$\displaystyle\Delta S^{\rm tot}(t,t^{\prime})\coloneqq\Delta S^{\rm tot}(t^{\prime})-\Delta S^{\rm tot}(t)=\int_{t}^{t^{\prime}}dt~{}\sigma^{\rm tot}(t),$
	$\displaystyle\Theta(t,t^{\prime})\coloneqq\Theta(t^{\prime})-\Theta(t)=\int_{t}^{t^{\prime}}dt~{}\Lambda(x(t),t)\circ\dot{x}(t).$		(39)

Then, the TUR is given by

\displaystyle\Delta S^{\rm tot}(t,t^{\prime})\geq\frac{2\langle\Theta(t,t^{\prime})\rangle^{2}}{\langle\Delta\Theta(t,t^{\prime})^{2}\rangle},

(40)

where $\langle\Theta(t,t^{\prime})\rangle$ and $\Delta\Theta(t,t^{\prime})$ are

	$\displaystyle\langle\Theta(t,t^{\prime})\rangle$	$\displaystyle=\int_{t}^{t^{\prime}}dt\int dx~{}\Lambda(x,t)J^{\rm irr}(x,t),$
	$\displaystyle\Delta\Theta(t,t^{\prime})$	$\displaystyle=\Theta(t,t^{\prime})-\langle\Theta(t,t^{\prime})\rangle.$		(41)

In the short-time limit $\delta t\rightarrow 0$ , $\langle\Theta(t,t^{\prime})\rangle=\langle\dot{\Theta}(t)\rangle\delta t$ and $\langle\Delta\Theta(t,t^{\prime})^{2}\rangle=\langle(\Lambda\circ\delta x-\langle\dot{\Theta}(t)\rangle\delta t)^{2}\rangle=2\langle\Lambda B\Lambda\rangle_{x}\delta t+O(\delta t^{2})$ . With these short-time forms, Eq. (40) becomes identical to the rate EB equation (11). Therefore, the previous EP estimation methods using the multidimensional TUR Dechant (2019); Van Vu et al. (2020) are identical to our MEB method for a 1D overdamped Langevin system in the short-time limit. For a higher-dimensional process, it is not possible to write the component-wise TUR, such as $\Delta S_{i}(\tau)\geq 2\frac{\langle\Theta_{i}(\tau)\rangle^{2}}{\langle\Delta\Theta_{i}(\tau)^{2}\rangle}$ , and consequently we cannot obtain the component-wise bound for the EP from the TUR. In this sense, the MEB provides more detailed information on the EP, compared to the TUR method.

We note that other modified TURs with an arbitrary time-dependent protocol or an arbitrary initial state do not approach the rate EB in the short-time limit even in one dimension. As an example, consider a 1D overdamped Langevin system driven by a time-dependent protocol. The modified TUR for this process with an arbitrary protocol was introduced by Koyuk and Seifert Koyuk and Seifert (2020) as

\displaystyle\Delta S^{\rm tot}(t,t^{\prime})\geq~{}

\displaystyle\frac{2\left[\hat{h}(t^{\prime})\langle\Theta(t,t^{\prime})\rangle\right]^{2}}{\langle\Delta\Theta(t,t^{\prime})^{2}\rangle}\equiv\Delta S^{\rm KS}(t,t^{\prime}),

(42)

where $\hat{h}(t)=t\partial_{t}-\omega\partial_{\omega}$ and $\omega$ denotes the protocol speed. In the $\delta t\rightarrow 0$ limit, the numerator of $\Delta S^{\rm KS}(t,t^{\prime})$ becomes $2\{(1-\omega\partial_{\omega})\langle\dot{\Theta}(t)\rangle\}^{2}\delta t^{2}$ , and thus this TUR is written as

\displaystyle\sigma^{\rm tot}(t)\geq\frac{2\{(1-\omega\partial_{\omega})\langle\dot{\Theta}(t)\rangle\}^{2}}{\langle\Lambda B\Lambda\rangle_{x}}\equiv\sigma^{KS}(t),

(43)

which is different from the rate EB in general. Note that we have to evaluate the sub-leading-order contribution when the numerator vanishes in Eq. (43). Experimental estimation of $\sigma^{KS}(t)$ is a very laborious task since we need a sufficiently large ensemble of trajectories, slightly perturbed with respect to $\omega$ at time $t$ for every single $t$ . Therefore, compared to this modified TUR method, MEB is a much more efficient approach to correctly estimate the EP of a system driven by a time-dependent protocol.

In addition, short-time TURs for underdamped dynamics do not correspond to the rate EB either. The TUR for a 1D underdamped system with a time-dependent protocol and the observable current $\langle\dot{\Theta}(t)\rangle=\langle\Lambda(x,t)v\rangle_{x}$ can be written as Lee et al. (2021)

\displaystyle\Delta S^{\rm tot}(t,t^{\prime})\geq~{}

\displaystyle\frac{2\left[\hat{h}_{\rm u}(t^{\prime})\langle\Theta(t,t^{\prime})\rangle\right]^{2}}{\langle\Delta\Theta(t,t^{\prime})^{2}\rangle}-I(t),

(44)

where $\hat{h}_{\rm u}(t)=t\partial_{t}-s\partial_{s}-r\partial_{r}-\omega\partial_{\omega}$ and $I(t)$ is the initial-state dependent term, defined as $I(t)=2\langle(1+\hat{h}_{\rm u}^{\prime}\ln P(x,t))^{2}\rangle$ with $\hat{h}_{\rm u}^{\prime}=x\partial_{x}-s\partial_{s}-r\partial_{r}-\omega\partial_{\omega}$ . Here, $s$ and $r$ are scaling parameters for force and position, respectively. In the $\delta t\rightarrow 0$ limit, since $\Theta(t,t^{\prime})=\Lambda v\delta t$ , $\langle\Delta\Theta(t,t^{\prime})^{2}\rangle=\langle(\Lambda v-\langle\dot{\Theta}(t)\rangle)^{2}\rangle\delta t^{2}$ , which is not $O(\delta t)$ Fischer et al. (2020). Thus, in the $\delta t\rightarrow 0$ limit, Eq. (44) becomes

\displaystyle\sigma^{\rm tot}(t)\delta t\geq\frac{2\left[(1-s\partial_{s}-r\partial_{r}-\omega\partial_{\omega})\langle\dot{\Theta}(t)\rangle\right]^{2}}{\langle(\Lambda v-\langle\dot{\Theta}(t)\rangle)^{2}\rangle}-I(t),

(45)

which is also different from the rate EB. For evaluating Eq. (45) experimentally, slight scalings of all forces and position variables are necessary, which demands full knowledge and full controllability of all forces. Thus, it is clear that the underdamped TUR is not practically useful to estimate the EP for a complicated system, such as complex biological systems where such detailed information is not available.

We conclude that the MEB is a unified tool that enables the efficient estimation of the EP from a trajectory ensemble for an overdamped or underdamped Langevin process without an irreversible force. Finally, it is interesting to note that the integral MEB can be tight when we choose the optimal observable current, a feature that no finite-time TUR can achieve.

IV Numerical Verification of MEB

In this section, we estimate the EP of three physical systems driven by time-dependent protocols via the MEB method. All these systems can be realized experimentally using optical tweezers. The first example is a dragged Brownian particle, the second is a pulled harmonic chain, and the last is an RNA unfolding process. We also compare the MEB results to those of other well-established EP measurement methods, namely the direct method utilizing $\Delta S^{\rm tot}=\Delta S^{\rm sys}+Q/T$ and a machine learning technique (NEEP) Kim et al. (2020); Otsubo et al. (2020); Kim et al. (2022).

IV.1 Brownian particle dragged by optical tweezers

We consider a 1D Brownian particle dragged by optical tweezers as illustrated in Fig. 1. The center of the harmonic potential of the tweezers is initially ( $t\leq 0$ ) located at the origin and moves with a constant speed $\omega$ for $t>0$ . Then the corresponding overdamped Langevin equation for the position $x(t)$ is written as

\displaystyle\dot{x}(t)=-\mu k(x(t)-\lambda(t))+\sqrt{2B}\xi(t),

(46)

where $\lambda(t)=\omega t$ is a time-dependent protocol, $\mu$ is the mobility, $k$ is the spring constant of the harmonic potential, and $B=\mu T$ with the environmental temperature $T$ . The initial state at $t=0$ is set as the equilibrium state. As the driving force is linear in position, we can solve the equation of motion analytically. The procedure for deriving the analytic solutions is presented in Appendix B.

Refer to caption — Figure 1: Plot for the rate EP estimator $\tilde{\sigma}$ normalized with respect to the total EP rate $\sigma^{\rm tot}(t)$ as a function of time $t$ . The green solid and red dotted line denote the MEB results of the overdamped and underdamped dynamics, respectively. The orange dashed line represents the result of the modified TUR by Koyuk and Seifert, $\sigma^{\rm KS}(t)/\sigma^{\rm tot}(t)$ . The parameters used for this plot are $k=\mu=\omega=T=1$ . The inset shows a schematic of the Brownian particle dragged by optical tweezers. Note that these plots are obtained from the analytic expressions.

We measure the displacement of the particle, that is, we take the weight vector $\Lambda(x,t)=1$ . With this observable current, we evaluate the rate MEB estimator $\sigma^{\rm MEB}(t)$ for each time $t$ analytically and plot the results in Fig. 1. Note that the normalized estimator is defined by $\tilde{\sigma}\equiv\sigma^{\rm MEB}(t)/\sigma^{\rm tot}(t)$ , which turns out to be unity, i.e., the estimated EP exactly matches the true EP. This is a rather surprising result, as we use only one current (displacement). In fact, one can analytically find the tight weight factor $\Lambda^{\rm e}(x,t)$ in Eq. (13) with the help of the exact solution in Eq. (65) as

\displaystyle\Lambda^{\rm e}(x,t)=c(t)\frac{J^{\rm irr}(x,t)}{BP(x,t)}=c(t)\frac{\omega}{\mu T}(1-e^{-t/\tau_{\mu}}),

(47)

where $\tau_{\mu}=1/(\mu k)$ . Note that $\Lambda^{\rm e}(x,t)$ is $x$ -independent. Thus, by choosing the arbitrary $c(t)$ to cancel the $t$ -dependence exactly in Eq. (47), one can easily see that the unity weight vector $\Lambda^{\rm e}=1$ also satisfies the equality condition of the rate EB.

For the purpose of comparison, we also plot the ratio of the modified TUR by Koyuk and Seifert Koyuk and Seifert (2020) to the EP rate in Fig. 1. For evaluating the ratio of this example, the sub-leading-order terms that we neglected for deriving Eq. (43) are necessary since the numerator in Eq. (43) vanishes in this example, so we calculate $\Delta S^{KS}(t+\delta t,t)/\Delta S(t+\delta t,t)$ with small $\delta t=0.001$ . This approximated ratio coincides with the result in Ref. Koyuk and Seifert (2020). The modified TUR deviates largely from the correct one for small $t$ . This confirms that our MEB method outperforms the modified TUR method for this simple case.

We also consider the same process in the underdamped version. The corresponding underdamped Langevin equation is written as

	$\displaystyle\dot{x}(t)=$	$\displaystyle v(t)$
	$\displaystyle\dot{v}(t)=$	$\displaystyle-\frac{1}{m\mu}v(t)-\frac{k}{m}(x(t)-\lambda(t))+\sqrt{2B}\xi(t),$		(48)

where $\lambda(t)=\omega t$ and $B=T/(\mu m^{2})$ . The initial state is also set as the equilibrium state. The analytic solution of this equation is also available via similar procedure for solving Eq. (46). The derivation is presented in Appendix B. From Eq. (71), the tight weight vector is obtained as

\displaystyle\bm{\Lambda}^{\rm e}(x,v,t)=c(t)\frac{\bm{J}^{\rm irr}(x,v,t)}{BP(x,v,t)}=c(t)\begin{pmatrix}0\\ -\frac{m\langle v(t)\rangle}{T}\end{pmatrix},

(49)

where $\langle v(t)\rangle$ is evaluated by taking the time derivative of $\langle x(t)\rangle$ in Eq. (69). We find that $\bm{\Lambda}^{\rm e}(x,v,t)$ depends only on time but not position, like in the overdamped case. Therefore, the unit weight vector again provides the EP exactly. The analytic result of $\tilde{\sigma}=\sigma^{\rm MEB}(t)/\sigma^{\rm tot}(t)$ for this underdamped dynamics is also plotted in Fig. 1, which confirms the exact estimation of the EP from the rate MEB by measuring only the displacement.

IV.2 Harmonic chain pulled by optical tweezers

The next example is an $M$ -bead harmonic chain dragged by optical tweezers as illustrated in Fig. 2(a). The harmonic potential of the optical tweezers is exerted on the rightmost particle of the chain, and the leftmost spring clings to the wall. Here, we consider an overdamped Langevin dynamics described by

\displaystyle\dot{\bm{x}}(t)=-\mu\mathsf{K}\bm{x}(t)+\mu k\bm{\lambda}(t)+\sqrt{2\mathsf{B}}\bm{\xi}(t),

(50)

where $K_{ij}=2k\delta_{i,j}-k(\delta_{i+1,j}+\delta_{i-1,j})$ , $\lambda_{i}=\omega t\delta_{M,i}$ , and $B_{ij}=\mu T\delta_{i,j}$ with $i,j\in\{1,...,M\}$ . We can solve Eq. (50) in a similar way used for the dragged Brownian particle. The derivation details are presented in Appendix B.

For validating the MEB estimator, we generate $10^{5}$ trajectories of $\tau=1$ by solving Eq. (50) numerically with $M=5$ and $T=1$ via the 2nd-order stochastic differential equation integrator. We set the time resolution $dt$ to $0.01$ . The initial state of the chain is in equilibrium, with the center of the harmonic potential being located at the origin. From the trajectories, we estimate $\Delta S_{i}$ , the EP for the $i$ -th bead, by using the rate MEB estimator with the $\ell=4$ Gaussian weight vector set. The estimated data are plotted in the inset of Fig. 2(b). As the figure shows, the estimated EP of each bead perfectly matches the analytic result. We plot the total EP by adding all these $\Delta S_{i}$ in Fig. 2(b). For comparison, we also estimate the total EP with the NEEP machine learning technique Kim et al. (2020), which coincides with the MEB result precisely. The detailed procedure for the NEEP calculation is explained in Appendix E. Both methods exactly estimate the total EP within $0.5\%$ error.

Figure 2(c) is a plot of the total EP at $t=1$ against the number of weight vectors. Surprisingly, the EP estimated by the MEB with only one weight vector is already very close to the analytic result. This is due to the fact that a constant weight vector results in the exact EP value in this system, as explained Appendix B. As a Gaussian function with a broad width can be approximated as a constant, the EP can be approximately estimated solely with the broadest Gaussian function. The EP for $\ell>1$ saturates to the analytic result as expected in Sec. II.2. In Fig. 2(c), we also plot the result of the NEEP calculation, which is also close to the analytic result.

To test the performance of the MEB method for a high-dimensional system, we perform the same simulation with $M=100$ . The estimated EP from the MEB method is $14.4$ , which is only $3.5\%$ distant from the analytic result, $13.9$ . We note that this accurate estimation is infeasible for the direct (plug-in) method for such a high-dimensional system.

IV.3 RNA unfolding process

The final example is an RNA unfolding process, which involves a nonlinear potential and thus an analytic treatment is not possible. A typical experimental setup consists of a single RNA hairpin molecule whose terminals are connected to DNA handles that are controlled by two optical tweezers, as illustrated in Fig. 3(a). By moving the center of the rightmost optical tweezers, a pulling force is exerted on the rightmost particle and the RNA is unfolded. For the RNA hairpin P5GA Kim et al. (2012), a pulling force that amounts to $14.7$ pN yields equal probabilities for folded and unfolded states. The governing equation of motion in this case is

\displaystyle\dot{x}(t)=f_{14.7}(x(t))+\omega t+\sqrt{2D}\xi(t),

(51)

where $x(t)$ is the distance between the two ends of the RNA at time $t$ . The force function $f_{14.7}$ is estimated from coarse-grained molecular dynamics simulation data when the RNA is pulled by a $14.7$ pN force Kim et al. (2012), where a polynomial function of degree $10$ is employed to fit the force. The reflection boundary condition is imposed at $x_{\rm min}=1.01$ nm and $x_{\rm max}=9.07\,\text{nm}$ , as distances larger than $x_{\rm max}$ and smaller than $x_{\rm min}$ were not found in the simulation Kim et al. (2012). We set the initial condition as the equilibrium state at room temperature $300$ K, which is an ordinary experimental setup. During the process time $\tau=7.19~{}\text{ns}$ , the pulling force linearly increases up to $19.7\,\text{pN}$ with a constant rate $\omega=\frac{5.0}{7.19}$ pN/ns. We generate $10^{5}$ trajectories from the simulations. The time resolution $dt$ is set to 79.1 ps. The initial and final distributions at $t=0$ and $\tau$ are presented in Fig. 3(a), respectively.

We estimate the total EP by evaluating the rate MEB from the trajectory ensemble. Here we use the Gaussian weight vector set from Eq. (28) for evaluating the rate MEB estimator. Figure 3(b) shows a plot of the estimated EP as a function of time for $\ell=4$ . As the analytic expression of the EP for this system is not available, we evaluate the EP using other numerical methods to check the validity of the MEB method. First, we employ the NEEP using the same trajectory ensemble and present the result in Fig. 3(b). We find that the MEB and NEEP results coincide with each other precisely. Second, we apply the direct method using the following equality:

\displaystyle\Delta S^{\rm tot}(\tau)=\Delta S^{\rm sys}(\tau)+\frac{1}{T}\int_{0}^{\tau}\left\langle(f_{14.7}(x)+\omega t)\circ\dot{x}\right\rangle dt,

(52)

where $\Delta S^{\rm sys}(\tau)=\langle-\ln{p(x,\tau)+\ln{p(x,0)}}\rangle$ . The integration over the process time of the last term in Eq. (52) denotes the dissipated heat during the process. The initial and the final PDFs can be estimated from the trajectory ensemble. This task becomes much harder with increasing system dimension. The estimated total EP from the direct method is denoted as the red solid line in Fig. 3(b), which matches the MEB and NEEP results very well. We note that the computational cost of the MEB method is lower than that of the NEEP method; it takes 3 s for the MEB method with $4$ weight vectors, while it takes 60 s for the NEEP process including the learning time ¹¹1Intel i9-11900K is used for MEB evaluation. The same CPU and an RTX 3090 are used for NEEP calculation.. Here, we do not take into account the time required to find the proper hyperparameters for the NEEP.

Figure 3(c) shows the way to determine the proper number of weight vectors $\ell$ . From a given trajectory ensemble, we estimate the total EP by using the MEB estimator; the estimated EP as a function of $\ell$ for this RNA unfolding process is plotted in Fig. 3(c). For $\ell<3$ , the estimated EP increases as $\ell$ increases, which indicates that no combinations of two Gaussian weight vectors, Eq. (23), are sufficiently close to the optimal weight vector, Eq. (13). The estimated EP saturates to a certain value for $\ell\geq 3$ , which indicates that the estimator is now sufficiently close to the optimal one. Thus, accurate EP estimation can be obtained by choosing $\ell\geq 3$ . We also examine the dependence of the time resolution of an experiment and a limited number of trajectories on the EP in Appendix D.

V Conclusion and Discussion

In this study, we suggested an EP estimator, named MEB, by applying multidimensional observable currents to the entropic bound. The MEB provides a unified way to estimate the EP for both overdamped and underdamped Langevin dynamics regardless of the time dependence of the protocol. The MEB estimator can be obtained in either integral or rate form. The tight EP bound is always achievable for any finite-time processes via both the integral and the rate MEBs, whereas it is possible for TURs only in the short-time limit. From numerical simulations, we confirmed that the MEB estimates the EP with high accuracy from an ensemble of system trajectories of overdamped Langevin systems. For an underdamped system with an irreversible force, information about the force is additionally required to estimate the EP. Moreover, extra information on the relaxation time is necessary for underdamped systems. Therefore, a precise estimation of the EP may be possible via MEB even for various complicated physical processes, in particular biological systems.

In future research, it will be interesting to develop a method to estimate the stochastic EP at the level of a single trajectory for general Langevin dynamics, rather than the averaged EP over an ensemble of trajectories. Moreover, the extension of the EP estimation to an open quantum system will be another intriguing problem. The quantum TUR recently proposed in Ref. Van Vu and Saito (2022) could be a good candidate for an estimator of the EP of an open quantum system, if one can measure the coherent-effect term in the formulation. It is also worthwhile to mention the recently proposed method for directly inferring the stochastic differential equations from a given trajectory ensemble Brückner et al. (2020); Frishman and Ronceray (2020). It would be interesting to compare the accuracy and efficiency of EP estimation between MEB and the inferring method.

Acknowledgements.

This research was supported by an NRF Grant No. 2017R1D1A1B06035497 (H.P.) and KIAS Individual Grants Nos. PG081801 (S.L.), PG074002 (J.-M.P.), CG076002 (W.K.K.), QP013601 (H.P.), and PG064901 (J.S.L.) at the Korea Institute for Advanced Study. D.-K.K. was supported by the Institute for Basic Science of Korea (IBS-R029-C2).

Appendix A Derivation of Eq. (27)

Here, we focus on the derivation of the integral MEB. The derivation for the rate MEB is essentially the same as that of the integral MEB. ${\mathsf{L}}_{i}^{(\ell+1)}$ can be expressed as the following block matrix form:

\displaystyle{\mathsf{L}}_{i}^{(\ell+1)}=\begin{bmatrix}{\mathsf{L}}_{i}^{(\ell)}&\bm{b}\\ \bm{b}^{\mathsf{T}}&h\end{bmatrix}

(53)

where $\bm{b}^{\mathsf{T}}=[({\mathsf{L}}_{i})_{\ell+1,1},\cdots,({\mathsf{L}}_{i})_{\ell+1,\ell}]$ and $h=({\mathsf{L}}_{i})_{\ell+1,\ell+1}=\int^{\tau}_{0}dt\langle\Lambda_{i,\ell+1}\mathsf{B}_{i}\Lambda_{i,\ell+1}\rangle_{\bm{q}}$ . From the Schur complement, the determinant of the block matrix $L_{i}^{(\ell+1)}$ in Eq. (53) is given by

\displaystyle{\rm{det}}({\mathsf{L}}_{i}^{(\ell+1)})={\rm{det}}({\mathsf{L}}_{i}^{(\ell)})\left[h-\bm{b}^{\mathsf{T}}({\mathsf{L}}_{i}^{(\ell)})^{-1}\bm{b}\right].

(54)

The determinant of $\mathsf{L}_{i}^{(\ell)}$ for any $\ell$ is positive since it is a positive-definite matrix. This implies that the term $h-\bm{b}^{\mathsf{T}}({\mathsf{L}}^{(l)})^{-1}\bm{b}$ in Eq. (54) is also positive. Moreover, via the following inverse block matrix formula,

\displaystyle\begin{bmatrix}\mathsf{A}&\mathsf{B}\\ \mathsf{C}&\mathsf{D}\end{bmatrix}^{-1}=

\displaystyle\begin{bmatrix}\mathsf{A}^{-1}+\mathsf{A}^{-1}\mathsf{B}(\mathsf{D}-\mathsf{C}\mathsf{A}^{-1}\mathsf{B})^{-1}\mathsf{C}\mathsf{A}^{-1}&-\mathsf{A}^{-1}\mathsf{B}(\mathsf{D}-\mathsf{C}\mathsf{A}^{-1}\mathsf{B})^{-1}\\ -(\mathsf{D}-\mathsf{C}\mathsf{A}^{-1}\mathsf{B})^{-1}\mathsf{C}\mathsf{A}^{-1}&(\mathsf{D}-\mathsf{C}\mathsf{A}^{-1}\mathsf{B})^{-1}\end{bmatrix},

(55)

the inverse matrix of $\mathsf{L}_{i}^{(\ell+1)}$ can be expressed as

\displaystyle({\mathsf{L}}_{i}^{(\ell+1)})^{-1}=

\displaystyle\begin{bmatrix}({\mathsf{L}}_{i}^{(\ell)})^{-1}&0\\ 0&0\end{bmatrix}+(h-\bm{b}^{\mathsf{T}}({\mathsf{L}}_{i}^{(\ell)})^{-1}\bm{b})^{-1}\bm{d}\bm{d}^{\mathsf{T}},

(56)

where $\bm{d}^{\mathsf{T}}\equiv(-\bm{b}^{\mathsf{T}}({\mathsf{L}}_{i}^{(\ell)})^{-1},1)$ . Using Eq. (56), we can prove Eq. (27) as follows:

$\displaystyle\Delta S_{i}^{{\rm MEB}(\ell+1)}=$	$\displaystyle\langle{\bm{\Theta}}_{i}^{(\ell+1)}\rangle^{\mathsf{T}}({\mathsf{L}}_{i}^{(\ell+1)})^{-1}\langle{\bm{\Theta}}_{i}^{(\ell+1)}\rangle$
$\displaystyle=$	$\displaystyle\Delta S_{i}^{{\rm MEB}(\ell)}+(h-{\bm{b}}^{\mathsf{T}}({\mathsf{L}}_{i}^{(\ell)})^{-1}{\bm{b}})(\bm{d}^{\mathsf{T}}\langle\bm{\Theta}_{i}^{(\ell+1)}\rangle)^{2}$
$\displaystyle\geq$	$\displaystyle\Delta S_{i}^{{\rm MEB}(\ell)}.$	(57)

The positiveness of $h-{\bm{b}}^{\mathsf{T}}({\mathsf{L}}_{i}^{(\ell)})^{-1}{\bm{b}}$ is used for showing the last inequality of Eq. (57).

Appendix B Analytic solutions of a dragged Brownian particle and pulled harmonic chain by optical tweezers

We consider a Brownian particle dragged by optical tweezers of which dynamics is governed by the following equation Koyuk and Seifert (2020); van Zon and Cohen (2003):

\displaystyle\dot{x}(t)=-\mu k(x(t)-\lambda(t))+\sqrt{2B}\xi(t),

(58)

where $B=\mu T$ and $\lambda(t)$ is an arbitrary time-dependent protocol with the condition $\lambda(0)=0$ . The initial state is set as the equilibrium distribution of Eq. (58) with $\lambda(0)=0$ , and thus, $\langle x(0)\rangle=0$ . To obtain the analytic solution of Eq. (58), we decompose $x(t)$ into the deterministic part $\langle x(t)\rangle$ and the stochastic part $X(t)\equiv x(t)-\langle x(t)\rangle$ . Taking the average of both sides of Eq. (58) leads to an equation for the deterministic part as

\displaystyle\langle\dot{x}(t)\rangle=-\mu k(\langle x(t)\rangle-\lambda(t)).

(59)

Then, the solution of $\langle x(t)\rangle$ is given by

	$\displaystyle\langle x(t)\rangle=$	$\displaystyle e^{-t/\tau_{\mu}}\langle x(0)\rangle+\tau_{\mu}^{-1}\int^{t}_{0}dt^{\prime}e^{-(t-t^{\prime})/\tau_{\mu}}\lambda(t^{\prime})$
	$\displaystyle=$	$\displaystyle\lambda(t)-\int^{t}_{0}dt^{\prime}e^{-(t-t^{\prime})/\tau_{\mu}}\dot{\lambda}(t^{\prime}),$		(60)

where $\tau_{\mu}=(\mu k)^{-1}$ is a characteristic relaxation time. The equation for the stochastic component $X(t)$ can be obtained by simply substituting $X(t)+\langle x(t)\rangle$ for $x(t)$ in Eq. (58) as

\displaystyle\dot{X}(t)=-\tau_{\mu}^{-1}X(t)+\sqrt{2B}\xi(t).

(61)

Since the initial state is in equilibrium, the distribution of $X(t)$ does not change in time. Therefore, the distribution for all time is given by the equilibrium distribution as

\displaystyle P(X,t)=\sqrt{\frac{\beta k}{2\pi}}\,e^{-\frac{k}{2}\beta X^{2}}.

(62)

By substituting $x-\langle x(t)\rangle$ for $X$ in Eq. (62), we have

\displaystyle P(x,t)=\sqrt{\frac{\beta k}{2\pi}}\,e^{-\frac{k}{2}\beta(x-\langle x(t)\rangle)^{2}}.

(63)

Using Eqs. (5) and (60), the irreversible current is

$\displaystyle J^{\rm irr}(x,t)=$	$\displaystyle\left[-\mu k(x-\lambda(t))-B\partial_{x}\right]P(x,t)$
$\displaystyle=$	$\displaystyle\mu k[\lambda(t)-\langle x(t)\rangle]P(x,t)$
$\displaystyle=$	$\displaystyle\mu k\int^{t}_{0}dt^{\prime}e^{-(t-t^{\prime})/\tau_{\mu}}\dot{\lambda}(t^{\prime})P(x,t).$	(64)

When $\lambda(t)=\omega t$ as in Sec. IV.1, the irreversible current and EP rate are

$\displaystyle J^{\rm irr}(x,t)=$	$\displaystyle\omega[1-e^{-t/\tau_{\mu}}]P(x,t),$	(65)
$\displaystyle\sigma^{\rm tot}(t)=$	$\displaystyle\int dx\frac{{J}^{\rm irr}(x,t)^{2}}{BP(x,t)}$
$\displaystyle=$	$\displaystyle\frac{\omega^{2}\beta}{\mu}(1-e^{-t/\tau_{\mu}})^{2}.$	(66)

The derivation procedure for the underdamped Langevin equation (48) is essentially the same as that of the overdamped equation. By decomposing $x(t)$ into $\langle x(t)\rangle$ and $X(t)=x(t)-\langle x(t)\rangle$ and $v(t)$ into $\langle v(t)\rangle$ and $V(t)=v(t)-\langle x(t)\rangle$ , we have

	$\displaystyle\frac{d^{2}}{dt^{2}}\langle x(t)\rangle$	$\displaystyle=-\frac{1}{m\mu}\frac{d}{dt}\langle x(t)\rangle-\frac{k}{m}(\langle x(t)\rangle-\omega t),$		(67)
	$\displaystyle\frac{d}{dt}V(t)$	$\displaystyle=-\frac{1}{m\mu}V(t)-\frac{k}{m}X(t)+\sqrt{2B}\xi(t).$		(68)

For this underdamped case, $B=T\mu^{-1}m^{-2}$ . The second-order differential equation (67) can be solved with the boundary conditions $\langle x(0)\rangle=0$ and $\langle v(0)\rangle=d/dt\langle x(t)\rangle|_{t=0}=0$ . The result is

\displaystyle\langle x(t)\rangle=C_{+}e^{\alpha_{+}t}+C_{-}e^{\alpha_{-}t}+\omega t-\frac{\omega}{\mu k},

(69)

where $\alpha_{\pm}=-1/(2m\mu)\pm\sqrt{1/(2m\mu)^{2}-k/m}$ and $C_{\pm}=\mp\omega(\alpha_{\mp}/(\mu k)+1)/(\alpha_{+}-\alpha_{-})$ . As the initial state is in equilibrium, the distribution of $X(t)$ and $V(t)$ in Eq. (68) for all time is the following equilibrium distribution:

\displaystyle P(X,V,t)=\sqrt{\frac{\beta k}{2\pi}}\sqrt{\frac{\beta m}{2\pi}}\exp\left[-\frac{\beta}{2}(kX^{2}+mV^{2})\right].

(70)

Therefore, $P(x,v,t)$ is given by substituting $x-\langle x(t)\rangle$ for $X$ and $v-\langle v(t)\rangle$ for $V$ in Eq. (70), as was done in Eq. (63). From Eq. (5), the irreversible current is written as

\displaystyle J_{x}^{\rm irr}(x,v,t)=0,~{}~{}~{}J_{v}^{\rm irr}(x,v,t)=-\frac{\langle v(t)\rangle}{m\mu}P(x,v,t).

(71)

Finally, the EP rate is evaluated as

\displaystyle\sigma^{\rm tot}(t)=\int dx\frac{{J}^{\rm irr}(x,v,t)^{2}}{BP(x,v,t)}=\frac{\langle v(t)\rangle^{2}}{\mu T}.

(72)

The analytic solution of Eq. (50) can be obtained in a similar way. By decomposing $x_{i}(t)$ into $X_{i}(t)=x_{i}(t)-\langle x_{i}(t)\rangle$ and rearranging the terms of Eq. (50), we have

	$\displaystyle\langle\dot{\bm{x}}(t)\rangle$	$\displaystyle=-\mu\mathsf{K}\bm{\langle}\bm{x}(t)\rangle+\mu k\bm{\lambda}(t),$		(73)
	$\displaystyle\dot{\bm{X}}(t)$	$\displaystyle=-\mu\mathsf{K}\bm{X}(t)+\sqrt{2\mathsf{B}}\bm{\xi}(t).$		(74)

The expression of $\langle x_{i}(t)\rangle$ can be obtained by solving Eq. (73), and it is certain that $\langle x_{i}(t)\rangle$ is a function of time. Since the initial state is the equilibrium state of Eq. (74), the probability density function (PDF) is given by the Boltzmann factor $\exp[-\beta U(\bm{X})]$ , where $U(\bm{X})=\frac{1}{2}\bm{X}^{\textsf{T}}\textsf{K}\bm{X}$ is the potential energy of the harmonic chain. Thus, by substituting $\bm{X}=\bm{x}-\langle\bm{x}\rangle$ into the Boltzmann factor, the PDF is written as

	$\displaystyle P(\bm{x},t)=$	$\displaystyle\frac{e^{-\frac{\beta}{2}(\bm{x}-\langle\bm{x}(t)\rangle)^{\mathsf{T}}\mathsf{K}({\bm{x}}-\langle\bm{x}(t)\rangle)}}{\sqrt{\det(2\pi\textsf{K}^{-1}/\beta)}},$		(75)
	$\displaystyle\bm{J}^{\rm irr}(\bm{x},t)=$	$\displaystyle[\mu k\bm{\lambda}(t)-\mu\mathsf{K}\langle\bm{x}(t)\rangle]P(\bm{x},t).$		(76)

The tight weight vector is then given by

	$\displaystyle\Lambda_{i}^{\rm e}(\bm{x},t)$	$\displaystyle=c_{i}(t)\frac{J_{i}^{\rm irr}(\bm{x},t)}{\mu TP(\bm{x},t)}$
		$\displaystyle=c_{i}(t)\left(\mu k\lambda_{i}(t)-\mu\sum_{j}K_{ij}\langle x_{j}(t)\rangle\right).$		(77)

Note that $J_{i}^{\rm irr}(\bm{x},t)/P(\bm{x},t)$ depends on time but not position. Thus, the MEB estimator evaluated by measuring displacement, i.e., $\Lambda_{i}^{\rm e}(\bm{x},t)=1$ , results in the correct EP.

Appendix C Component-combined MEB

In this section, we present the derivation of the component-combined MEB, which is useful when the diffusion matrix has off-diagonal elements. To this end, similar to Eq. (14), we consider the following linear combination of weight vectors as

\displaystyle\bm{\Lambda}^{(\ell)}(\bm{q},t)=\sum_{\alpha=1}^{\ell}k_{\alpha}\bm{\Lambda}_{\alpha}(\bm{q},t).

(78)

After substituting $\bm{\Lambda}$ in Eq. (10) with $\bm{\Lambda}^{(\ell)}(\bm{q},t)$ in Eq. (78), we have

	$\displaystyle\Delta S^{\rm tot}(t)\geq$	$\displaystyle\frac{\left[\int^{\tau}_{0}dt\int d\bm{q}\bm{\Lambda}^{(\ell)}(\bm{q},t)^{\textsf{T}}\bm{J}^{\rm irr}(\bm{q},t)\right]^{2}}{\int^{\tau}_{0}dt\langle\bm{\Lambda}^{(\ell)}(\bm{q},t)^{\textsf{T}}\mathsf{B}\bm{\Lambda}^{(\ell)}(\bm{q},t)\rangle_{\bm{q}}}$
	$\displaystyle=$	$\displaystyle\frac{(\bm{k}^{\textsf{T}}\langle\bm{\Theta}^{(\ell)}(\tau)\rangle)^{2}}{\bm{k}^{\textsf{T}}\mathsf{L}^{(\ell)}(\tau)\bm{k}}=\Delta{\hat{S}}^{(\ell)}(\bm{k}),$		(79)

where the components of $\langle\bm{\Theta}^{(\ell)}(\tau)\rangle$ and $\mathsf{L}^{(\ell)}(\tau)$ are given as

	$\displaystyle\langle\Theta^{(\ell)}_{\alpha}(\tau)\rangle=$	$\displaystyle\int^{\tau}_{0}dt\int d\bm{q}\bm{\Lambda}_{\alpha}(\bm{q},t)^{\textsf{T}}\bm{J}^{\rm irr}(\bm{q},t),$		(80)
	$\displaystyle(\mathsf{L}^{(\ell)}(\tau))_{\alpha,\beta}=$	$\displaystyle\int^{\tau}_{0}dt\langle\bm{\Lambda}_{\alpha}(\bm{q},t)^{\textsf{T}}\mathsf{B}\bm{\Lambda}_{\beta}(\bm{q},t)\rangle_{\bm{q}}.$		(81)

Note that $\mathsf{L}^{(\ell)}(\tau)$ is a positive definite matrix since $\bm{z}\mathsf{L}^{(\ell)}(\tau)\bm{z}=\int^{\tau}_{0}dt\langle||\sum_{\alpha}\mathsf{B}^{1/2}\bm{\Lambda}_{\alpha}(\bm{q},t)z_{\alpha}||^{2}\rangle_{\bm{q}}>0$ for an arbitrary non-zero vector $\bm{z}$ , positive definite matrix $\mathsf{B}$ , and non-zero vector $\bm{\Lambda}_{\alpha}$ . Then, the optimal condition for $\Delta{\hat{S}}^{(\ell)}(\bm{k})$ can be written as

\displaystyle\nabla_{\bm{k}}\Delta{\hat{S}}^{(\ell)}(\bm{k})=\frac{2(\bm{k}^{\textsf{T}}\langle\bm{\Theta}^{(\ell)}(\tau)\rangle)[(\bm{k}^{\textsf{T}}\mathsf{L}^{(\ell)}(\tau)\bm{k})\langle\bm{\Theta}^{(\ell)}(\tau)\rangle-(\bm{k}^{\textsf{T}}\langle\bm{\Theta}^{(\ell)}(\tau)\rangle)(\mathsf{L}^{(\ell)}(\tau)\bm{k})]}{[\bm{k}^{\textsf{T}}\mathsf{L}^{(\ell)}(\tau)\bm{k}]^{2}}=0,

(82)

which is similar to Eq. (19). The solution of Eq. (82) is $\bm{k}^{*}=(\mathsf{L}^{(\ell)}(\tau))^{-1}\cdot\langle\bm{\Theta}^{(\ell)}(\tau)\rangle$ . By plugging this $\bm{k}^{*}$ into Eq. (79), we obtain the component-combined MEB as

\displaystyle\Delta S^{\rm tot}\geq

\displaystyle\Delta\hat{S}^{(\ell)}(\bm{k}^{*})=\langle\bm{\Theta}^{(\ell)}(\tau)\rangle(\mathsf{L}^{(\ell)}(\tau))^{-1}\langle\bm{\Theta}^{(\ell)}(\tau)\rangle.

(83)

This is the integral form of the component-combined MEB. Following a similar way, we can derive the rate form of the component-combined MEB as

\displaystyle\sigma^{\rm tot}(t)\geq\langle\dot{\bm{\Theta}}^{(\ell)}(\tau)\rangle^{\textsf{T}}(\dot{\mathsf{L}}^{(\ell)}(\tau))^{-1}\langle\dot{\bm{\Theta}}^{(\ell)}(\tau)\rangle,

(84)

where the components of $\langle\dot{\bm{\Theta}}^{(\ell)}(\tau)\rangle$ and $\dot{\mathsf{L}}^{(\ell)}(\tau)$ are given by

	$\displaystyle\langle\dot{\Theta}^{(\ell)}_{\alpha}(t)\rangle=$	$\displaystyle\int d\bm{q}\bm{\Lambda}_{\alpha}(\bm{q},t)^{\textsf{T}}\bm{J}^{\rm irr}(\bm{q},t),$		(85)
	$\displaystyle(\dot{\mathsf{L}}^{(\ell)}(t))_{\alpha,\beta}=$	$\displaystyle\langle\bm{\Lambda}_{\alpha}(\bm{q},t)^{\textsf{T}}\mathsf{B}\bm{\Lambda}_{\beta}(\bm{q},t)\rangle_{\bm{q}}.$		(86)

Similar to Eq. (36), in the case where the estimation of the diffusion matrix requires a heavy computational cost, we can directly obtain $\dot{\mathsf{L}}^{(\ell)}(t)$ by evaluating

	$\displaystyle(\dot{\mathsf{L}}^{(\ell)}(t))_{\alpha,\beta}=$	$\displaystyle\langle\bm{\Lambda}_{\alpha}(\bm{q},t)^{\textsf{T}}\mathsf{B}\bm{\Lambda}_{\beta}(\bm{q},t)\rangle$
	$\displaystyle=$	$\displaystyle\lim_{\delta t\rightarrow 0}\frac{1}{\delta t}\langle[\bm{\Lambda}_{\alpha}(\bm{q},t)^{\textsf{T}}\circ\delta\bm{q})][\bm{\Lambda}_{\beta}(\bm{q},t)^{\textsf{T}}\circ\delta\bm{q}]\rangle.$		(87)

Then, $({\mathsf{L}}^{(\ell)}(\tau))_{\alpha,\beta}$ can be estimated by integrating Eq. (87) over time from $t=0$ to $t=\tau$ .

Appendix D Effect of limited samples or time resolution on the EP estimation

D.1 Limited number of trajectories

Though we use $10^{5}$ trajectories in our simulation in Sec. IV.3, only several thousand repetitions are usually feasible in real experiments Jun et al. (2014). Thus, in this section we examine the effect of a limited number of trajectories on the estimated EP. To this end, we perform additional simulations of the RNA unfolding process for various trajectory numbers of 1000, 2000, 4000, 8000, and 10000 and estimate the EP using both the MEB and NEEP methods. The results are plotted in Fig. 4. Open and filled circles denote the MEB and NEEP results, respectively. To plot the error bars, we first generated 5 independent data sets for each number of samples, and then evaluated the average and standard deviation for the 5 sets. The red solid line is the estimated EP via the direct method, which is the same line as in Fig. 3(c). The figure shows the tendency that both MEB and NEEP overestimate the EP for small sample sizes. In fact, we recommend using both methods together to cross-check the reliability of the estimated value.

D.2 Limited time resolution

The estimated EP depends on the time resolution $dt$ (time gap between two consecutive data points) of the measurement. Since decreasing the time resolution (increasing $dt$ ) causes a ‘coarse-graining’ of the trajectory data, doing so typically leads to a lower value of the EP. This can also be checked in our simulations; here, we simulate the RNA unfolding process presented in Sec. IV.3 with various $dt$ . Other parameters are the same as those in the previous simulation. The results are plotted in Fig. 5. The NEEP, the MEB, and the direct method show the same decreasing behavior as $dt$ increases. From an extrapolation of the data, we can also estimate the EP value in the $dt\rightarrow 0$ limit. Accordingly, it is important to specify the $dt$ information in a given experiment or simulation.

Appendix E NEEP algorithm

Here we explain the training details of the NEEP and its architecture configurations Kim et al. (2020). We apply the NEEP to one step from $\bm{x}_{t}$ to $\bm{x}_{t+\Delta t}$ . For brevity, we will use the notation $\bm{x}_{t}$ for $\bm{x}(t)$ . The NEEP is designed to maximize the following objective function:

\displaystyle C(\theta)\equiv\left<\Delta{S_{\theta}(\bm{x}_{t+\Delta t},\bm{x}_{t},t)}-e^{-\Delta S_{\theta}(\bm{x}_{t+\Delta t},\bm{x}_{t},t)}\right>

(88)

where $\Delta S_{\theta}$ is an antisymmetric function with respect to the exchange of $\bm{x}_{t}$ and $\bm{x}_{t+\Delta t}$ as

\displaystyle\Delta S_{\theta}(\bm{x}_{t+\Delta t},\bm{x}_{t},t)=h_{\theta}(\bm{x}_{t},\bm{x}_{t+\Delta t},t)-h_{\theta}(\bm{x}_{t+\Delta t},\bm{x}_{t},t).

(89)

In Eq. (89), the function $h_{\theta}$ is the output of a multi-layer perceptron (MLP) and $\theta$ denotes the trainable parameters of the MLP. The MLP has a scalar output unit and three hidden layers of 512 units with the rectified linear unit (ReLU) activation function. It is shown that $\Delta S_{\theta}=\Delta S^{\rm tot}$ with the optimized $\theta^{*}$ in Ref. Kim et al. (2020).

In order to employ the cross-validation method, we split the trajectory data into 20% for the validation set and 80% for the training set. We train the MLP $h_{\theta}$ to maximize Eq. (88) by using the Adam optimizer Kingma and Ba (2014) with learning rate $10^{-4}$ , batch size 4096, and weight decay $5\times 10^{-5}$ . Before feeding the input $(\bm{x}_{t+\Delta t},\bm{x}_{t},t)$ to the MLP, we normalize each element of trajectory data $\bm{x}$ by using the following equation:

\displaystyle x^{(i)}_{t}\leftarrow(x^{(i)}_{t}-\text{mean}[x^{(i)}])/\text{std}[x^{(i)}],

where $x^{(i)}$ indicates the $i$ -th component of $\bm{x}$ , $\text{mean}[x^{(i)}]$ is the mean of $x^{(i)}$ , and $\text{std}[x^{(i)}]$ is the standard deviation of $x^{(i)}$ . We also normalize the time information $t=0\dots\tau$ to be set as $t=-0.5\dots 0.5$ so that mean of the input vector $(\bm{x}_{t+\Delta t},\bm{x}_{t},t)$ becomes a zero vector. The total number of training iterations is $10^{4}$ , and we evaluate $C(\theta)$ values from the validation set per every 500 (50) training iterations for the pulled harmonic chain (RNA unfolding process). The best trained parameter set $\theta^{\star}$ is determined from the case where the NEEP produces the maximum value of $C$ during the training process. The results presented in Sec. IV are those evaluated at the best trained parameter $\theta^{\star}$ over the total trajectory data.

References

Jarzynski (1997) Christopher Jarzynski, “Nonequilibrium equality for free energy differences,” Phys. Rev. Lett. 78, 2690 (1997).
Crooks (1999) Gavin E Crooks, “Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences,” Phys. Rev. E 60, 2721 (1999).
Seifert (2019) Udo Seifert, “From stochastic thermodynamics to thermodynamic inference,” Annu. Rev. Condens. Matter Phys. 10, 171 (2019).
Seifert (2005) Udo Seifert, “Entropy production along a stochastic trajectory and an integral fluctuation theorem,” Phys. Rev. Lett. 95, 040602 (2005).
Collin et al. (2005) D. Collin, F. Ritort, C. Jarzynski, S. B. Smith, I. Tinoco, and C. Bustamante, “Verification of the crooks fluctuation theorem and recovery of rna folding free energies,” Nature 437, 231–234 (2005).
Barato and Seifert (2015) Andre C Barato and Udo Seifert, “Thermodynamic uncertainty relation for biomolecular processes,” Phys. Rev. Lett. 114, 158101 (2015).
Gingrich et al. (2016) Todd R. Gingrich, Jordan M. Horowitz, Nikolay Perunov, and Jeremy L. England, “Dissipation bounds all steady-state current fluctuations,” Phys. Rev. Lett. 116, 120601 (2016).
Horowitz and Gingrich (2017) Jordan M Horowitz and Todd R Gingrich, “Proof of the finite-time thermodynamics uncertainty relation for steady-state currents,” Phys. Rev. E 96, 020103 (2017).
Dechant and Sasa (2018a) Andreas Dechant and Shin-Ichi Sasa, “Current fluctuations and transport efficiency for general langevin systems,” J. Stat. Mech.: Theor. Exp. , 063209 (2018a).
Hasegawa and Vu (2019) Yoshihiko Hasegawa and Tan Van Vu, “Uncertainty relation in stochastic processes: An information inequality approach,” Phys. Rev. E 99, 062126 (2019).
Dechant (2019) Andreas Dechant, “Multidimensional thermodynamic uncertainty relations,” J. Phys. A:Math. Theor. 52, 035001 (2019).
Hasegawa and Van Vu (2019) Yoshihiko Hasegawa and Tan Van Vu, “Fluctuation theorem uncertainty relation,” Phys. Rev. Lett. 123, 110602 (2019).
Koyuk and Seifert (2020) Timur Koyuk and Udo Seifert, “Thermodynamic uncertainty relation for time-dependent driving,” Phys. Rev. Lett. 125, 260604 (2020).
Benenti et al. (2011) Giuliano Benenti, Keiji Saito, and Giulio Casati, “Thermodynamic bounds on efficiency for systems with broken time-reversal symmetry,” Phys. Rev. Lett. 106, 230602 (2011).
Brandner et al. (2013) Kay Brandner, Keiji Saito, and Udo Seifert, “Strong bounds on onsager coefficients and efficiency for three-terminal thermoelectric transport in a magnetic field,” Phys. Rev. Lett. 110, 070603 (2013).
Shiraishi et al. (2016) Naoto Shiraishi, Keiji Saito, and Hal Tasaki, “Universal trade-off relation between power and efficiency for heat engines,” Phys. Rev. Lett. 117, 190601 (2016).
Dechant and Sasa (2018b) Andreas Dechant and Shin-ichi Sasa, “Entropic bounds on currents in langevin systems,” Phys. Rev. E 97, 062101 (2018b).
Lee et al. (2020) Jae Sung Lee, Jong-Min Park, Hyun-Myung Chun, Jaegon Um, and Hyunggyu Park, “Exactly solvable two-terminal heat engine with asymmetric onsager coefficients: Origin of the power-efficiency bound,” Phys. Rev. E 101, 052132 (2020).
Shiraishi et al. (2018) Naoto Shiraishi, Ken Funo, and Keiji Saito, “Speed limit for classical stochastic processes,” Phys. Rev. Lett. 121, 070601 (2018).
Nicholson et al. (2020) Schuyler B. Nicholson, Luis Pedro García-Pintos, Adolfo del Campo, and Jason R. Green, “Time–information uncertainty relations in thermodynamics,” Nature Physics 16, 1211–1215 (2020).
Vo et al. (2020) Van Tuan Vo, Tan Van Vu, and Yoshihiko Hasegawa, “Unified approach to classical speed limit and thermodynamic uncertainty relation,” Phys. Rev. E 102, 062132 (2020).
Falasco and Esposito (2020) Gianmaria Falasco and Massimiliano Esposito, “Dissipation-time uncertainty relation,” Phys. Rev. Lett. 125, 120604 (2020).
Yoshimura and Ito (2021) Kohei Yoshimura and Sosuke Ito, “Thermodynamic uncertainty relation and thermodynamic speed limit in deterministic chemical reaction networks,” Phys. Rev. Lett. 127, 160601 (2021).
Lee et al. (2022) Jae Sung Lee, Sangyun Lee, Hyukjoon Kwon, and Hyunggyu Park, “Speed limit for a highly irreversible process and tight finite-time landauer’s bound,” Phys. Rev. Lett. 129, 120603 (2022).
Roldán and Parrondo (2010) Édgar Roldán and Juan M. R. Parrondo, “Estimating dissipation from single stationary trajectories,” Phys. Rev. Lett. 105, 150607 (2010).
Roldán and Parrondo (2012) Édgar Roldán and Juan MR Parrondo, “Entropy production and kullback-leibler divergence between stationary trajectories of discrete systems,” Phys. Rev. E. 85, 031129 (2012).
Battle et al. (2016) Christopher Battle, Chase P Broedersz, Nikta Fakhri, Veikko F Geyer, Jonathon Howard, Christoph F Schmidt, and Fred C MacKintosh, “Broken detailed balance at mesoscopic scales in active biological systems,” Science 352, 604–607 (2016).
Avinery et al. (2019) Ram Avinery, Micha Kornreich, and Roy Beck, “Universal and accessible entropy estimation using a compression algorithm,” Phys. Rev. Lett. 123, 178102 (2019).
Martiniani et al. (2019) Stefano Martiniani, Paul M. Chaikin, and Dov Levine, “Quantifying hidden order out of equilibrium,” Phys. Rev. X 9, 011031 (2019).
Li et al. (2019) Junang Li, Jordan M. Horowitz, Todd R. Gingrich, and Nikta Fakhri, “Quantifying dissipation using fluctuating currents,” Nat. Commun. 10, 1666 (2019).
Manikandan et al. (2020) Sreekanth K. Manikandan, Deepak Gupta, and Supriya Krishnamurthy, “Inferring entropy production from short experiments,” Phys. Rev. Lett. 124, 120603 (2020).
Martínez et al. (2019) Ignacio A Martínez, Gili Bisker, Jordan M Horowitz, and Juan MR Parrondo, “Inferring broken detailed balance in the absence of observable currents,” Nat. Commun. 10, 3542 (2019).
Kim et al. (2020) Dong-Kyum Kim, Youngkyoung Bae, Sangyun Lee, and Hawoong Jeong, “Learning entropy production via neural networks,” Phys. Rev. Lett. 125, 140604 (2020).
Otsubo et al. (2020) Shun Otsubo, Sosuke Ito, Andreas Dechant, and Takahiro Sagawa, “Estimating entropy production by machine learning of short-time fluctuating currents,” Phys. Rev. E 101, 062106 (2020).
Otsubo et al. (2022) Shun Otsubo, Sreekanth K. Manikandan, Takahiro Sagawa, and Supriya Krishnamurthy, “Estimating time-dependent entropy production from non-equilibrium trajectories,” Communications Physics 5, 11 (2022).
Skinner and Dunkel (2021) Dominic J. Skinner and Jörn Dunkel, “Estimating entropy production from waiting time distributions,” Phys. Rev. Lett. 127, 198101 (2021).
Kim et al. (2022) Dong-Kyum Kim, Sangyun Lee, and Hawoong Jeong, “Estimating entropy production with odd-parity state variables via machine learning,” Phys. Rev. Research 4, 023051 (2022).
van der Meer et al. (2022) Jann van der Meer, Benjamin Ertel, and Udo Seifert, “Thermodynamic inference in partially accessible markov networks: A unifying perspective from transition-based waiting time distributions,” arXiv preprint arXiv:2203.12020 (2022).
Seif et al. (2021) Alireza Seif, Mohammad Hafezi, and Christopher Jarzynski, “Machine learning the thermodynamic arrow of time,” Nature Physics 17, 105–113 (2021).
Seifert (2012) Udo Seifert, “Stochastic thermodynamics, fluctuation theorems and molecular machines,” Rep. Prog. Phys. 75, 126001 (2012).
Sekimoto (2010) Ken Sekimoto, Stochastic energetics (Springer-Verlag, 2010).
Lee et al. (2019) Jae Sung Lee, Jong-Min Park, and Hyunggyu Park, “Thermodynamic uncertainty relation for underdamped langevin systems driven by a velocity-dependent force,” Phys. Rev. E 100, 062132 (2019).
Lee et al. (2021) Jae Sung Lee, Jong-Min Park, and Hyunggyu Park, “Universal form of thermodynamic uncertainty relation for langevin dynamics,” Phys. Rev. E 104, L052102 (2021).
Donsker and Varadhan (2020) M. D. Donsker and S. R. S. Varadhan, “Asymptotic evaluation of certain markov process expectations for large time. iv,” Commun. Pure Appl. Math. 36, 183 (2020).
Risken (1991) H. Risken, The Fokker-Planck Equation: Methods of Solutions and Applications (Springer-Verlag, Berlin, 1991).
Spinney and Ford (2012) Richard E. Spinney and Ian J. Ford, “Entropy production in full phase space for continuous stochastic dynamics,” Phys. Rev. E 85, 051113 (2012).
Kwon et al. (2016) Chulan Kwon, Joonhyun Yeo, Hyun Keun Lee, and Hyunggyu Park, “Unconventional entropy production in the presence of momentum-dependent forces,” J. Korean Phys. Soc. 68, 633 (2016).
Dechant and Sasa (2021) Andreas Dechant and Shin-ichi Sasa, “Improving thermodynamic bounds using correlations,” Phys. Rev. X 11, 041061 (2021).
Van Vu et al. (2020) Tan Van Vu, Van Tuan Vo, and Yoshihiko Hasegawa, “Entropy production estimation with optimal current,” Phys. Rev. E 101, 042138 (2020).
Lee and Park (2018) Jae Sung Lee and Hyunggyu Park, “Additivity of multiple heat reservoirs in the langevin equation,” Phys. Rev. E 97, 062135 (2018).
Fischer et al. (2020) Lukas P. Fischer, Hyun-Myung Chun, and Udo Seifert, “Free diffusion bounds the precision of currents in underdamped dynamics,” Phys. Rev. E 102, 012120 (2020).
Kim et al. (2012) Won Kyu Kim, Changbong Hyeon, and Wokyung Sung, “Weak temporal signals can synchronize and accelerate the transition dynamics of biopolymers under tension,” Proceedings of the National Academy of Sciences 109, 14410–14415 (2012).
Note (1) Intel i9-11900K is used for MEB evaluation. The same CPU and an RTX 3090 are used for NEEP calculation.
Van Vu and Saito (2022) Tan Van Vu and Keiji Saito, “Thermodynamics of precision in markovian open quantum dynamics,” Phys. Rev. Lett. 128, 140602 (2022).
Brückner et al. (2020) David B. Brückner, Pierre Ronceray, and Chase P. Broedersz, “Inferring the dynamics of underdamped stochastic systems,” Phys. Rev. Lett. 125, 058103 (2020).
Frishman and Ronceray (2020) Anna Frishman and Pierre Ronceray, “Learning force fields from stochastic trajectories,” Phys. Rev. X 10, 021009 (2020).
van Zon and Cohen (2003) R. van Zon and E. G. D. Cohen, “Stationary and transient work-fluctuation theorems for a dragged brownian particle,” Phys. Rev. E 67, 046102 (2003).
Jun et al. (2014) Yonggun Jun, Mom čilo Gavrilov, and John Bechhoefer, “High-precision test of landauer’s principle in a feedback trap,” Phys. Rev. Lett. 113, 190601 (2014).
Kingma and Ba (2014) Diederik P. Kingma and Jimmy Ba, “Adam: A method for stochastic optimization,” arXiv:1412.6980 (2014).