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

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 without knowing the details of a given system. MEB is a unified method in the sense that it is applicable to both overdamped and underdamped Langevin dynamics, irrespective of the time dependence of the driving protocol. In addition, 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 [1][2][3][4][5], 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) [6][7][8][9][10][11][12][13], the power-efficiency trade-off relation [14][15][16][17][18], and the speed limit [19][20][21][22][23][24].
Subsequently, experimentally feasible methods for measuring the EP have been actively suggested and discussed [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39].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 [40], 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, ∆S tot .The first method uses the equality ∆S tot = ∆S sys +Q/T , where ∆S sys is the Shannon entropy change of a system and Q is dissipated heat into a reservoir at temperature T [40,41].In experiments, it is difficult to measure the * jslee@kias.re.kr 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 [41].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. [30].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 ∆S tot ≥ B(Θ), where the EP bound B(Θ) is determined by an observable current Θ.In a certain condition, one can find an optimal observable current Θ * , which yields B(Θ * ) = ∆S tot .Then, the EP can be accurately estimated by measuring Θ * .
Regarding the above indirect methods, there exist two representative inequalities.The first inequality is in TUR form [6][7][8][9][10], where the EP is bounded by the relative fluctuation of a certain observable current.To access a tighter bound of this TUR, multidimensional TUR [11] and Monte Carlo methods [30] have been developed.However, TURs depend on the nature of the system dynamics; e.g, the TUR must be modified when a timedependent protocol is involved [13] or when a system follows underdamped Langevin dynamics [12,42,43].Thus, EP estimation based on TURs is not universal.More-arXiv:2207.05961v2[cond-mat.stat-mech]22 Feb 2023 over, 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 [12,42,43].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 [44].Recently, a machine learning technique named NEEP (neural estimator for entropy production) [33,34,37] 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 [37].
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 [17].Inspired by the multidimensional TUR [11], 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 [17].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.
A. Derivation of the integral and the rate EB Here, we consider an M -dimensional Langevin system with a state vector q(t) = (q 1 , • • • , q M ) T , where T denotes the transpose of a matrix, described by the following equation of motion: where The symbol • in Eq. ( 1) represents the Itô product.From now on, we sometimes drop the arguments of functions for simplicity.
A component of q can be an odd-parity variable such as a velocity under time-reversal operation.The time reversal of a state q is denoted by q = (q 1 , • • • , qM ) T with qi = i q i , where i = 1 for an even-parity variable and i = −1 otherwise.The drift force can be divided into reversible and irreversible parts as A(q, t) = A rev (q, t) + A irr (q, t) with [45] A rev (q, t) ≡ and † is an operation changing the sign of the odd-parity parameters.
The Fokker-Planck equation associated with Eq. ( 1) is with the PDF P (q, t).The reversible current J rev (q, t) and the irreversible current J irr (q, t) are defined as J irr i (q, t) ≡ A irr i (q, t)P (q, t) with B( q, t) = B(q, t).Note that for an overdamped Langevin system with even-parity variables only, A rev (q, t) vanishes, and thus J rev (q, t) = 0 and the total current coincides with J irr (q, t).As dissipation originates from the irreversible current, the EP is determined only by J irr (q, t).Therefore, the total EP rate σ tot is given by [17,46,47] σ tot (t) ≡ dq J irr (q, t) T B(q, t) −1 J irr (q, t) P (q, t) .
Hereafter, we use the k B = 1 unit.Note that for underdamped Langevin systems, the matrix is not directly invertible since B ij = 0 when the component index i or j denotes a positional variable.For such an index i, we first set B ii = b (b > 0) and B ij = 0 (i = j), then take the inverse of B and calculate J irr (q, t) T B(q, t) −1 J irr (q, t) in Eq. ( 6), and finally take the b → 0 limit.Since J irr i (q, t) ∝ b, this limit leads to J irr i (q, t) 2 B ii (q, t) −1 ∼ b → 0. For underdamped systems, this procedure amounts to writing B and J 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 τ : where T is a weight vector of the irreversible current for a given observable.Then, the averaged current rate at time t is given as The EB in an integral form can be derived from Eq. ( 7) as follows: B(q, t) − 1 2 J irr (q, t) P (q, t) where • • • q = dq • • • P (q, t) and the total EP ∆S tot (τ ) = τ 0 dt σ 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 (integral EB) (10) Similarly, the EB in a rate form can also be obtained from Eq. ( 8) as The equality of the integral EB is satisfied when the weight vector has the following form: Λ e (q, t) = c B(q, t) −1 J irr (q, t) P (q, t) (for the integral EB), (12) where c is an arbitrary constant that is independent of 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., Θ(τ ) = c∆S tot (τ ).Similarly, we find the equality condition for the rate EB as Λ e (q, t) = c(t) B(q, t) −1 J irr (q, t) P (q, t) (for the rate EB), (13) where c(t) is an arbitrary time-dependent function that is independent of q.This weight vector corresponds to the observable current rate as Θ(t) = c(t)σ 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 A.

B. Derivation of the integral and the rate MEB
With the knowledge of the functional form of Λ e (q, t), one may obtain the tight EP bound.However, except for very simple examples, it is impossible to identify Λ e (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 [11] 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 Λ e (q, t).The linear combination of weight vectors where k i,α is the coefficient for Λ i,α (q, t) and is independent of q and t.From now on, we will consider the case B ij = B i δ 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 B ij is given.Thus, we can still set B ij = B i δ ij on the transformed coordinate.In cases where it is difficult to obtain the information of 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 "componentcombined" form as we show in Appendix C. The observable in Eq. ( 7) can be divided into the sum of its components as where ∆S i (τ ) = τ 0 dtσ i (t) with the i-th component EP rate σ i (t) = dq B i (q, t) −1 J irr i (q, t) 2 /P (q, t).Thus, ∆S tot (τ ) = i ∆S i (τ ) and σ tot (t) = i σ i (t).By substituting Eq. ( 14) into Eq.( 15), we have where T , and the vector Θ respectively.Note that L ( ) The bound ∆ Ŝ( ) i (k i ) in Eq. ( 16) is a function of k i ; thus, the tightest bound can be written as ∆ where k * i is the optimal vector maximizing the bound.The optimal vector is obtained by solving the following equation: We can easily check that the numerator vanishes with the choice of k * i = (L 16), we find the component-wise MEB as By summing over all components, we finally obtain our main result, namely MEB in integral form, as follows: 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 In this case, it is usually sufficient to choose a timeindependent basis as Λ where the time-dependence is encoded in the coefficients instead of in Λ i,α (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 where the vector Θ( ) The total EP during a finite time τ can be evaluated by integrating σ MEB( ) (t) over time.
The weight vectors for the rate MEB are not timedependent, so we need a lower number of weight vectors to approximate Λ e i (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 observables, the MEB becomes tighter, i.e., ∆S The proof of Eq. ( 27) is basically the same as that presented in Ref. [48].To be self-contained, we include the proof in Appendix A. As we increase , the MEB estimator also increases and eventually saturates to the maximum value, i.e., ∆S i (τ ) or σ i (t).It can often saturate even at finite = sat for simple systems.Therefore, by observing the saturation regime in a plot of the EP estimator versus , 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 weight vectors as long as they are linearly independent of each other.In this study, we adopt a Gaussian weight vector set [49] 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: In Eq. ( 28), a i,α and b i,α are given as where q max i (q min i ) is the maximum (minimum) value of q i in a given trajectory ensemble and ∆q i = q max i − q min i .Explicit forms of a i,α and b i,α are given as a . Here, x = m, where m is an integer satisfying m ≤ x < 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 2 i , • • • , q i } typically overestimates the EP with large fluctuations when is large, since the polynomial term with a large exponent is highly sensitive to rare-event data.

C. 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., q(t) = x(t) = (x 1 , • • • , x M ), and the Langevin equation is written as where F i is a force applied to the x i component.The reversible current J rev i = 0, while the irreversible current is given as In the case of underdamped dynamics, the system states consist of both position and velocity variables, i.e., q(t)

and the Langevin equation is written as
As done in Eq. ( 2), the external force F i can be divided into reversible and irreversible parts as Then, the irreversible currents are given as while the reversible currents J rev xi (q, t) = v i P (q, t) and J rev vi (q, t) = 1 m F rev i (q, t)P (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.

Determination of Bi and Li
For an overdamped dynamics described by Eq. ( 29), the diffusivity B i is determined from the average of shorttime mean square displacements as where δx i (t) ≡ x i (t+δt)−x i (t) and • • • (q,t) denotes the average over the trajectory ensemble at position 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 δt , (underdamped) (35) where δv i (t) ≡ v i (t + δt) − v i (t).With these estimations for B i , we calculate ( L( ) i (t)) α,β from Eq. (26).Note that lim δt→0 δxi(t) 2 (q,t) δt = 0 in the underdamped case.In the case where the estimation of B i (q, t) requires heavy computation, instead of evaluating the diffusivity, we can directly estimate the L( ) i (t) α,β matrix from the average of the products of two infinitesimal observable changes as When the diffusion matrix has non-zero off-diagonal terms, we first have to directly evaluate B ij = lim δt→0 δq i δq j (q,t) /δt and find the proper transformation to diagonalize the matrix.

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 q(t) = (x 1 , • • • , x M ), the i-th component of an observable current rate can be measured by averaging Λ i,α (q(t)) • ẋi (t) over the ensemble of system trajectories as where • denotes the Stratonovich product.This can be checked from the fact that H(q, t) • ẋi (q,t) = dq H(q, t)J irr i (q, t) with an arbitrary function H(q, t).For Λ i,α (q) = 1, the observable current is the displacement in the x i direction as Θ i,α (τ ) = x i (τ ) − x i (0).
In the underdamped dynamics with , by plugging Eq. (33) into Eq.( 25), we have When F irr i = 0 and Λ i,α has no explicit velocity dependence, the observable current is proportional to Λ i,α (q)v i q , similar to that of the overdamped system, Eq. ( 37).However, extra information γ/m is required to evaluate Eq. (38) in the underdamped case.This constant can be determined experimentally by measuring the velocity relaxation time [50].Otherwise, when F irr i = 0 and Λ i,α has an explicit velocity dependence, extra calculation of the last term in Eq. ( 38) is necessary.For the most general case with F irr i = 0 (velocity-dependent force), Θi,α (t) cannot be determined solely by system trajectories, but rather concrete information on the force is necessary.

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 ∆S tot (t, t ) during a small time segment between t and t = t + δt and the corresponding accumulated current Θ(t, t ) are written as Then, the TUR is given by where Θ(t, t ) and ∆Θ(t, t ) are In the short-time limit δt → 0, Θ(t, t ) = Θ(t) δt and ∆Θ(t, t ).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 [11,49] 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 ∆S i (τ ) ≥ 2 Θi(τ ) 2 ∆Θi(τ ) 2 , 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 [13] as ∆S tot (t, t ) ≥ 2 ĥ(t ) Θ(t, t ) 2 ∆Θ(t, t ) 2  ≡ ∆S KS (t, t ), (42) where ĥ(t) = t∂ t −ω∂ ω and ω denotes the protocol speed.In the δt → 0 limit, the numerator of ∆S KS (t, t ) becomes 2{(1 − ω∂ ω ) Θ(t) } 2 δt 2 , and thus this TUR is written as 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 σ KS (t) is a very laborious task since we need a sufficiently large ensemble of trajectories, slightly perturbed with respect to ω 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 timedependent protocol and the observable current Θ(t) = Λ(x, t)v x can be written as [43] ∆S tot (t, t ) ≥ where ĥu (t) = t∂ t − s∂ s − r∂ r − ω∂ ω and I(t) is the initial-state dependent term, defined as I(t) = 2 (1 + ĥ u ln P (x, t)) 2 with ĥ u = x∂ x − s∂ s − r∂ r − ω∂ ω .Here, s and r are scaling parameters for force and position, respectively.In the δt → 0 limit, since Θ(t, t ) = Λvδt, ∆Θ(t, t ) 2 = (Λv − Θ(t) ) 2 δt 2 , which is not O(δt) [51].Thus, in the δt → 0 limit, Eq. ( 44) becomes 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 wellestablished EP measurement methods, namely the direct method utilizing ∆S tot = ∆S sys + Q/T and a machine learning technique (NEEP) [33,34,37].

A. 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 ≤ 0) located at the origin and moves with a constant speed ω for t > 0. Then the corresponding overdamped Langevin equation for the position x(t) is written as where λ(t) = ωt is a time-dependent protocol, µ is the mobility, k is the spring constant of the harmonic potential, and B = µ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. We measure the displacement of the particle, that is, we take the weight vector Λ(x, t) = 1.With this observable current, we evaluate the rate MEB estimator σ MEB (t) for each time t analytically and plot the results in Fig. 1.Note that the normalized estimator is defined by σ ≡ σ MEB (t)/σ 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 Λ e (x, t) in Eq. ( 13) with the help of the exact solution in Eq. (B8) as where τ µ = 1/(µk).Note that Λ e (x, t) is x-independent.Thus, by choosing the arbitrary c(t) to cancel the tdependence exactly in Eq. ( 47), one can easily see that the unity weight vector Λ 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 [13] 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 ∆S KS (t + δt, t)/∆S(t + δt, t) with small δt = 0.001.This approximated ratio coincides with the result in Ref. [13].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 where λ(t) = ωt and B = T /(µ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. (B14), the tight weight vector is obtained as where v(t) is evaluated by taking the time derivative of x(t) in Eq. (B12).We find that Λ 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 σ = σ MEB (t)/σ 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.

B. 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 where , and B ij = µT δ i,j with i, j ∈ {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 τ = 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  52).(c) Estimated EP via MEB at t = 7.19 ns as a function of the number of weight vectors (black).The green dashed line and the red solid line denote the results of the NEEP and the EP obtained from Eq. ( 52), respectively.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 ∆S i , the EP for the i-th bead, by using the rate MEB estimator with the = 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 ∆S i in Fig. 2(b).For comparison, we also estimate the total EP with the NEEP machine learning technique [33], 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 > 1 saturates to the analytic result as expected in Sec.II B. 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.

C. RNA unfolding process
The final example is an RNA unfolding process, which involves a nonlinear potential and thus an analytic treat-ment 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 [52], 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 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 [52], where a polynomial function of degree 10 is employed to fit the force.The reflection boundary condition is imposed at x min = 1.01 nm and x max = 9.07 nm, as distances larger than x max and smaller than x min were not found in the simulation [52].We set the initial condition as the equilibrium state at room temperature 300 K, which is an ordinary experimental setup.During the process time τ = 7.19 ns, the pulling force linearly increases up to 19.7 pN with a constant rate ω = 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 τ 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 = 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: where ∆S sys (τ ) = − ln p(x, τ ) + ln p(x, 0) .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 [53].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 .From a given trajectory ensemble, we estimate the total EP by using the MEB estimator; the estimated EP as a function of for this RNA unfolding process is plotted in Fig. 3(c).For < 3, the estimated EP increases as 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 ≥ 3, which indicates that the estimator is now sufficiently close to the optimal one.Thus, accurate EP estimation can be obtained by choosing ≥ 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 es-timate 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. [54] 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 [55,56].It would be interesting to compare the accuracy and efficiency of EP estimation between MEB and the inferring method.
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 By substituting x − x(t) for X in Eq. (B5), we have Using Eqs. ( 5) and (B3), the irreversible current is When λ(t) = ωt as in Sec.IV A, the irreversible current and EP rate are The derivation procedure for the underdamped Langevin equation ( 48) is essentially the same as that of the overdamped equation.By decomposing x(t) into x(t) and X(t) = x(t) − x(t) and v(t) into v(t) and For this underdamped case, B = T µ −1 m −2 .The second-order differential equation (B10) can be solved with the boundary conditions x(0) = 0 and v(0) = d/dt x(t) | t=0 = 0.The result is where As the initial state is in equilibrium, the distribution of X(t) and V (t) in Eq. (B11) for all time is the following equilibrium distribution: Therefore, P (x, v, t) is given by substituting x − x(t) for X and v − v(t) for V in Eq. (B13), as was done in Eq. (B6).From Eq. ( 5), the irreversible current is written as Finally, the EP rate is evaluated as 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) − x i (t) and rearranging the terms of Eq. ( 50), we have ẋ(t) = −µK x(t) + µkλ(t), (B16) Ẋ(t) = −µKX(t) + √ 2Bξ(t).(B17) The expression of x i (t) can be obtained by solving Eq. (B16), and it is certain that x i (t) is a function of time.Since the initial state is the equilibrium state of Eq. (B17), the probability density function (PDF) is given by the Boltzmann factor exp[−βU (X)], where U (X) = 1 2 X T KX is the potential energy of the harmonic chain.Thus, by substituting X = x − x into the Boltzmann factor, the PDF is written as J irr (x, t) =[µkλ(t) − µK x(t) ]P (x, t). (B19) The tight weight vector is then given by Λ e i (x, t) = c i (t) J irr i (x, t) µT P (x, t) Note that J irr i (x, t)/P (x, t) depends on time but not position.Thus, the MEB estimator evaluated by measuring displacement, i.e., Λ e i (x, t) = 1, results in the correct EP.

FIG. 1 .
FIG.1.Plot for the rate EP estimator σ normalized with respect to the total EP rate σ 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, σ KS (t)/σ tot (t).The parameters used for this plot are k = µ = ω = 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.

FIG. 2 .
FIG. 2. (a) Schematic of the harmonic chain pulled by optical tweezers.(b) Estimated EP via the MEB method (black) and the NEEP method (green) as a function of time t.The inset shows the EP of the i-th bead., ×, +, , and ♦ represent the estimated EPs for x1, x2, x3, x4, and x5 beads, respectively.Solid lines denote the analytic results.Four Gaussian weight vectors are used for the MEB estimation.(c) Plot of ∆S tot at t = 1.0 against the number of weight vectors .The green dashed line represents the NEEP result, while the red solid lines in (b) and (c) denote the analytic results.The parameters used for these plots are k = 5, µ = 1, ω = 5, and T = 1.

FIG. 3 .
FIG. 3. EP estimation results for the RNA unfolding process.(a) Histogram of the distance x between the two ends of P5GA at initial time 0 ns (light blue) and final time 7.19 ns (orange).The inset shows a schematic of the RNA pulled by optical tweezers.(b) EP estimated via MEB with four Gaussian weight vectors (black) and NEEP (green) as a function of time t.The red solid line denotes the results obtained from Eq. (52).(c) Estimated EP via MEB at t = 7.19 ns as a function of the number of weight vectors (black).The green dashed line and the red solid line denote the results of the NEEP and the EP obtained from Eq. (52), respectively.

FIG. 4 .
FIG. 4. Estimated EP of the RNA unfolding process as a function of the number of trajectory samples.Open circles represent the MEB results with 4 Gaussian bases ( = 4).Green filled circles denote the NEEP results.The red solid line is the estimated EP via the direct method, which is the same line as in Fig.3(c).The other parameters are the same as those used to plot Fig.3(c).The error bars represent the standard deviation of the EP estimations from five different trajectory sets.