{'fraction_non_alphanumeric': 0.050491850883879634, 'fraction_numerical': 0.013721006207121855, 'mean_word_length': 4.941341384515851, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Skewed Thermodynamic Geometry and Optimal Free Energy Estimation

Free energy differences are a central quantity of interest in physics, chemistry, and biology. We develop design principles that improve the precision and accuracy of free energy estimators, which has potential applications to screening for targeted drug discovery. Specifically, we derive near-equilibrium approximations for moments of the excess work and find the leading- and next-to-leading order contributions to dissipation are from the Stokes' friction tensor and the supra-Stokes' tensor, with geodesics that respectively minimize the variance and bias of free energy estimators.

Introduction.-Free energy differences determine the equilibrium phases of thermodynamic systems as well as the relative reaction rates and binding affinities of chemical species [1]. Computational and experimental techniques which accurately and precisely predict free energy differences are therefore highly desirable. One important application is in pharmaceutical drug discovery, where computation of free energy differences can aid in the identification and design of ligands for targeted protein binding [2][3][4][5]. Current methods rely primarily on costly and time-consuming experimentation, which can be reduced through screening with efficient computational techniques [1][2][3][4][5][6][7][8].
The free energy difference between two states of a system is typically estimated by measuring the work incurred from a protocol (a dynamic variation in control parameters) that drives the system between controlparameter endpoints corresponding to each state. A unidirectional estimator calculates the free energy difference using only the work done by a protocol driving from initial to final states (forward protocol), while a bidirectional estimator also uses a protocol that drives from final to initial states (reverse protocol). The meanwork estimator [9] equates the mean work with the free energy difference and yields a biased estimate for any non-quasistatic (finite-speed) protocol. The Jarzynski estimator for free energy differences (derived from the Jarzynski equality relating the exponentially averaged work to the free energy change [10]) is unbiased for a large number of samples. The mean-work and Jarzynski estimators can be used as either uni-or bidirectional estimators; however, if bidirectional data is available the maximum log-likelihood estimate is Bennett's acceptance ratio (BAR) [11] which yields (for a large number of samples) the minimum variance of any unbiased estimator [12,13].
Regardless of the estimator, non-quasistatic controlparameter protocols result in excess work (work in excess of the free energy difference), which increases the bias and variance (decreases the accuracy and precision, respectively) of free energy estimates. This excess work (equaling the dissipation), and hence the error, can be reduced by following a different path through control-parameter space and varying the velocity along the path while keep-ing the protocol duration fixed [14][15][16]. Near equilibrium, this can be mapped onto the thermodynamic-geometry framework: for the mean-work estimator, the bias and variance from finite-time protocols is improved by following geodesics of a thermodynamic metric, the forcevariance (FV) metric [17][18][19][20][21]; similarly, for BAR, the variance (but not the bias) is minimized if the protocol follows geodesics of the FV metric [22]. This has been used to improve the precision of calculated binding potentials of mean force [23][24][25][26].
The FV metric only minimizes the variance of the mean-work estimator and BAR if the relaxation time is independent of the control parameters. The frictiontensor metric [27], a generalization of the FV, provides a relatively simple prescription for reducing excess work in more general near-equilibrium processes, where minimum-dissipation protocols follow geodesics of the Riemannian metric induced by the friction tensor. For common unidirectional estimators near equilibrium, the bias and variance are proportional to the first moment (mean) of the excess work, and hence minimumdissipation protocols defined by the friction tensor minimize both the bias and variance.
We show that for bidirectional estimators, the variance is proportional to the sum of the second moments of the excess work from forward and reverse protocols (25a), while the bias is proportional to the difference of the first moments (26a) [28]. With this in mind, we extend the thermodynamic-geometry framework to higher-order moments of the excess work and beyond the leading-order friction-tensor approximation. We find that for bidirectional estimators near equilibrium, minimum-variance protocols follow geodesics of the Riemannian metric induced by the friction tensor, while minimum-bias protocols follow geodesics of a cubic Finsler metric. For the simple model system of a Brownian particle in a quadratic trap with time-dependent stiffness (a breathing harmonic trap), the minimum-variance and minimumbias protocols can improve variance by a factor of 3 − 4 and bias by over a factor of 10 ( Fig. 2).
Derivation.-Consider a system in the canonical ensemble, in thermal equilibrium with a heat bath at temperature T . The probability distribution over microstates x at control parameters λ is π(x|λ) = exp β[F (λ) − E(x, λ)], with energy E(x, λ) and free en- We define the work done in a single realization of an external agent changing the control parameters λ according to a protocol Λ as W ≡ − t 0 dt f iλi (t ), which implies the average excess work (W ex ≡ W − ∆F ) is where we adopt the Einstein summation convention of implied summation over repeated indices. A dot denotes the time derivativeλ i ≡ dλ i /dt, f i ≡ −∂ λi U is the generalized force conjugate to λ i , and δf i ≡ f i − f i eq is the difference from the equilibrium average. Angle brackets · · · Λ denote an average over the nonequilibrium ensemble of system responses to control-parameter protocol Λ. The time derivative of the second moment of the excess-work distribution is For a sufficiently slow protocol, we replace the nonequilibrium average · · · Λ with the equilibrium average · · · λ(t) at fixed control parameters λ(t): where we used the stationarity of the equilibrium average, defined t ≡ t − t, and assumed smooth protocols to expand the control-parameter velocity to zeroth order,λ(t − t ) ≈λ(t). Finally, we assume correlations in the conjugate forces relax quickly relative to the protocol duration and replace the integration bound t with ∞ (Supplemental Material, Finite Bounds [29]), simplifying the approximation to for the friction tensor In analogy with fluid dynamics, this rank-two tensor is the Stokes' friction, since it produces a drag force that depends linearly on velocity. Following parallel arguments, we approximate the third moment of the excess-work distribution as for the rank-three tensor (The factor of three in (6) results from grouping the index permutations {ijk, jik, kij} into one term, using the invariance of the sum under exchange of indices, e.g., ζ (2) ijkλ iλjλk = ζ (2) jikλ iλjλk .) We call the rank-three tensor [Eq. (7)] the supra-Stokes' tensor, as it corresponds to the leading-order correction to dissipation beyond the Stokes' friction (15a).
For fourth and higher moments, Supplemental Material, Finite Bounds [29] shows that the integration bounds must remain finite since the n-time covariance functions do not decay to zero, so there is no clear analogy to frictional drag forces. Nevertheless, we can still approximate them by where we use index notation ν 1 , ν 2 , ν 3 · · · instead of i, j, k, and define integral n-time covariance functions .
This implies that higher-order moments of the excess work are higher order in control-parameter velocity and are therefore smaller for slow protocols. The approximations of (4), (6), and (8) are the leading-order contributions to each moment of the excess work.
To derive the next-order contribution we exploit the connection between time-reversed protocols through the Crooks relation [30,31], which constrains the probability of forward and reverse work measurements as Λ † is the time-reversed protocol, starting at equilibrium in the end state of the forward protocol. Integrating over W ex produces Jarzynski's equality [10] e −βWex Λ = 1 , while first multiplying by W ex then integrating leads to Taylor expanding the exponential in (11) and (12) to third order gives According to (8), near equilibrium the higher-order moments of the excess work are higher order in controlparameter velocity and therefore can be neglected for slow protocols. The leading-order contributions to the sum and difference of the excess work are Near equilibrium, the second and third moments are approximated by (4) and (6), so the time derivatives of the first and second moments are where henceforth we drop the explicit dependence of the friction and control-parameter velocity on λ and t.
To derive an approximation for higher-order moments, we follow parallel arguments. Multiplying (10) by W n ex and integrating over W ex , then either adding or subtracting (−W ex ) n Λ † , the leading-order contributions to the sum and difference for n > 2 are which implies Equations (15a), (15b) and (17) are the central results of this letter, which have applications to designing minimum-dissipation and optimal free energy estimation protocols.
Next-order contribution to excess work.-The approximation of (15a) can improve near-equilibrium estimates of the excess work. The leading-order term is the usual linear-response approximation (Eq. (13) of Ref. [27]) and is positive for all protocols. The Stokes' friction (5) is an autocovariance function, the Hadamard (entry-by-entry) product, •, of the conjugate-force covariance δf i δf j λ(t) and integral relaxation time The Stokes' friction is largest when the conjugate-force fluctuations are largest (large covariance) and most persistent (long relaxation time). In contrast to the positive contribution from the Stokes' friction, the contribution from the supra-Stokes' tensor (second term in (15a)) changes sign under time reversal because it is cubic in control-parameter velocity. The supra-Stokes' tensor (7) is not positive semidefinite, and is the Hadamard product of the unnormalized coskewness δf i δf j δf k λ(t) (related to skewness as covariance is related to variance) and the integral double relaxation time For a protocol with positive velocity, the contribution to dissipation from the supra-Stokes' tensor quantifies the increase (decrease) in excess work from negatively (positively) skewed conjugate-force fluctuation.
The leading-order friction-tensor approximation (first term in (15a)) endows the control-parameter space with a Riemannian metric, such that minimum-work protocols follow geodesics of the Stokes' friction tensor. With the addition of the next-order contribution (7), the excess work can be expressed as for total friction tensor that explicitly depends on the protocol velocity. Minimum-dissipation protocols follow geodesics of the generalized cubic Finsler metric ζ tot ij (Supplemental Material, Finsler Geometry [29]), an extension of the Riemannian metric ζ (1) ij . For small control-parameter velocities, the total friction is positive semi-definite; however, for large control-parameter velocities the approximation breaks down and no longer guarantees positive semidefiniteness, leading to the unphysical possibility of negative excess work.
To illustrate some general properties of the friction, Fig. 1 shows, for the model system of a breathing harmonic trap (Brownian particle in a harmonic trap where the control parameter is the time-dependent stiffness, see Supplemental Material, Breathing Trap [29]), the (a) force variance, (b) Stokes' friction, (c) supra-Stokes' contribution, and (d) total friction. The force variance and Stokes' friction are independent of control-parameter velocity. In general, the force variance differs from the Stokes' friction whenever integral relaxation time depends on the control parameters, as for the breathing harmonic trap. The contribution from the supra-Stokes' tensor (c) is antisymmetric in control-parameter velocity, becoming negative for negative velocity. This antisymmetric contribution skews the total friction (d), which i ζ tot /γ, for initial spring constant ki and damping coefficient γ. The control parameter is the trap stiffness, given in dimensionless form as λ * ≡ k/ki with velocitẏ λ * ≡ γk/k 2 i . Geodesics (solid curves) minimize the magnitude of the corresponding metric's contribution to the excess work, with distinct curves representing different average protocol velocities.
depends on the control-parameter velocity and lacks any symmetry under time reversal.
Geodesics (solid curves in Fig. 1) are protocols that minimize the contribution from the corresponding metric to the excess work. For relatively fast protocols (average scaled velocities ofλ * ≡ k/k i > ∼ 0.5 or < ∼ −0.5, for initial spring constant k i and damping coefficient γ), the velocity is significantly smaller in regions of high friction and larger in regions of low friction (or force variance).
Precision and accuracy of free energy estimates.-We quantify the precision of a free energy estimator by its variance. In the limit of many samples, the expected variance of any unbiased estimator ∆F is bounded by [12] δ ∆F where for simplicity we assume an equal number of forward and reverse work measurements, and the average · · · is over a total of N work measurements. BAR saturates this bound for large N . This bound demonstrates an explicit connection between the minimum variance of free energy estimators and the excess work. The variance is minimized at 0 when the excess-work distribution is a delta function at W ex = 0 (achieved in the quasistatic limit). For any finite-time protocol, the average excess work is positive.
Assuming small excess work, we expand the variance where the second line follows from (15b). Equation (25a) also holds for very few samples, since then BAR is equivalent to the average of the sum of the forward and reverse work measurements [28]. Thus the protocol designed to reduce the variance follows geodesics of ζ (1) ij , and for onedimensional control proceeds at velocityλ ∝ ζ (1) −1/2 .
Unlike the variance, the protocol that maximizes the accuracy (minimum-bias) is different for unidirectional and bidirectional estimators. For unidirectional Jarzynski and mean-work estimators, near equilibrium the minimum-bias protocol is simply the minimumdissipation protocol (protocol that minimizes (15a)) and therefore to leading order is optimized by the same protocol that minimizes (25b). For BAR, for small excess work (or few samples), the bias instead is proportional to the difference between the forward and reverse excess work [28]: where the second line follows from (15a). The protocol designed to reduce the (magnitude of) bias thus follows geodesics of the cubic Finsler metric ζ (2) ijk , simplifying for one-dimensional control toλ ∝ ζ (2) −1/3 .
To illustrate the potential benefit of protocols designed to reduce variance and bias, Fig. 2 shows approximations for the variance (25a) and bias (26a) of designed and naive protocols, as well as their ratio, as a function of protocol duration, for the model system of a breathing harmonic trap (Supplemental Material, Breathing Trap [29]). For a slow protocol (protocol duration longer than the slowest relaxation time, ∆t/τ (1) f > ∼ 1), the designed protocols reduce the variance by a factor of 3 − 4 (Fig. 2c), with the precise protocol (designed to reduce variance) performing the best (smallest ratio). For the bias, the accurate protocol (designed to reduce bias) performs the best (smallest ratio), with reductions by an order of magnitude for the slowest protocols shown (Fig. 2d). For fast protocols (protocol duration shorter than the slowest relaxation time, ∆t/τ (1) f < ∼ 1), the approximations break down, and the naive protocol can outperform the designed in both bias and variance (ratio larger than one). For this system the minimum-variance and minimum-bias protocols achieve similar amounts of bias and variance in all cases, likely due to their similar functional form (Supplemental Material, [29] Fig. 1).  (1) f = γ/(2k f )), for the breathing harmonic trap. Ratios of designed and naive (c) variances (δ ∆F ) 2 des / (δ ∆F ) 2 naive and (d) biases δ ∆F des / δ ∆F naive. Subscript 'naive' denotes the constantvelocity protocol, and 'des' the designed protocols, which are the force-variance-optimized 'force variance', the minimumvariance 'accurate', and the minimum-bias 'precise' protocols. Dash-dotted lines denote the limits for short (black) and long duration (colors). Final control-parameter value (spring constant) is k f /ki = 1/16.
Discussion.-We have developed near-equilibrium approximations for moments of the excess work, (15a), (15b), and (17), that incorporate time-reversal symmetric and antisymmetric contributions. The antisymmetric contribution to the first moment (15a) yields the nextorder contribution beyond linear response, which asymmetrically skews the Riemannian metric (5) into a generalized cubic Finsler metric (23). The Stokes' friction tensor (5) controls the leading-order contribution to the variance of both unidirectional and bidirectional free energy estimators (15b), such that minimum-variance protocols follow geodesics of the Stokes' friction. For unidirectional estimators, these same protocols also minimize the bias; however, for bidirectional estimators such as BAR the leading-order contribution to the bias is from the supra-Stokes' tensor (26b), and therefore minimum-bias protocols instead follow geodesics of the supra-Stokes' tensor. From these near-equilibrium approximations, we design protocols that increase the precision and accuracy of standard nonequilibrium free energy estimators (Fig. 2).
The addition of the supra-Stokes' tensor has several physical implications. Notably, it accounts for time-reversal asymmetric dissipation, which arises from skewed conjugate-force fluctuations. Since the supra-Stokes' contribution is antisymmetric, it cancels out for equal numbers of forward and reverse protocols, which has implications for the design principles of molecular machines: many molecular machines (e.g., kinesin walking toward a microtubule's plus end [32] or ATP synthase synthesizing ATP [33]) achieve directed behavior; the supra-Stokes' tensor quantifies the leading-order energetic cost of directed operation compared to a coequal forward and reverse process.
In some cases the minimum-work protocols can be calculated exactly [34,35], or minimum-variance protocols solved numerically [36,37], which would yield more accurate and precise estimators; however, to date these optimizations are limited to simple systems and do not permit straightforward generalization. The present formalism can be applied to more general settings than exact calculations and makes fewer approximations than Shenfeld et al. [22], who only examine the force-variance metric and do not consider minimizing bias.
Here we considered continuous protocols, but free energies are often estimated from sampling discrete states [38][39][40][41]. Additionally, for discrete state sampling in replica-exchange [42] or parallel-tempering [43] simulations it is important to consider the acceptance probability between states [44][45][46][47][48][49][50][51][52]. The force-variance metric has been applied in this context [22], and the linear-response framework of Ref. 27 has been generalized to discrete control [53]. As Ref. 22 has shown, spacing the states along geodesics of the force variance not only reduces the variance of the free energy estimator but also increases the acceptance probability. We expect that analogous principles hold for designing minimum-variance and minimumbias protocols for discrete state sampling using the milder approximations of our approach.

BREATHING HARMONIC TRAP
A simple non-trivial system for designing and testing minimum-variance and minimum-bias protocols is the breathing harmonic trap, since its correlation functions can be calculated analytically and it has a non-Gaussian work distribution. Consider a colloidal particle in a harmonic trap with variable stiffness. The particle position x obeys the overdamped Langevin equation, in a trap of time-dependent strength k, damping coefficient γ, environmental temperature T , Boltzmann's constant k B , and Gaussian white noise η. The control parameter is the trap strength λ = k, so the conjugate force is f ≡ −∂U/∂k = − 1 2 x 2 . The joint probability distribution of the particle position and work obeys [1,2] γ ∂p(x, w, t) ∂t From this, we can solve for any moment of either the work or position distribution, subject to an initial equilibrium condition p(x, w, t = 0) = π(x|k i )δ(w) = [2π/(βk i )] −1/2 exp{−βk i x 2 /2}δ(w), with δ(w) the Dirac delta function. The force variance (18), the Stokes' (5) and the supra-Stokes' coefficients (7), and the third integral covariance function (9) are respectively. In (S3d), we neglect terms of order t −n for n ≥ 1. Since we only consider one-dimensional control, we can analytically determine the minimum-variance (precise) and minimum-bias (accurate) protocols that minimize Eq. (15b) and (26b), respectively. The force-variance-optimized protocol proceeds at velocityk ∝ k, minimumvariance protocol atk ∝ k 3/2 , and minimum-bias protocol atk ∝ k 5/3 . Figure 1 plots these protocols alongside the naive (constant-velocity) protocol. As a simple test of the approximations of Eqs. (15a), (15b), and (17), Fig. 2 plots the sum and difference of the first three moments of the excess work from forward and reverse protocols. We choose the forward protocol to be a decrease in the control parameter (decreasing k), and therefore the reverse increases the control parameter (increasing k). In all cases the exact calculations agree with the approximation for slow protocols (large ∆t/τ (1) f for slowest relaxation time τ (1) f = γ/(2k f )), and overestimate for fast protocols (small ∆t/τ (1) f ). For large ∆t/τ (1) f , the first and second moments are approximated by (15a) and (15b) so the sums (n = 1 and 2 (+)) are both approximated by the Stokes' friction (5) and decrease as 1/∆t, while the differences (n = 1 and 2 (−)) are proportional to the supra-Stokes' tensor (7) and decrease as 1/(∆t) 2 . The third moment is approximated by (17), so the sum (n = 3 (+)) is approximated by the third integral covariance (9) and decays as 1/(∆t) 2 (due to the term linear in t in (S3d)) while the difference (n = 3 (−)) is proportional to the supra-Stokes tensor and decays as 1/(∆t) 2 .  (1) f = γ/(2k f )). Solid lines are the near-equilibrium approximations, given by Eqs. (15a), (15b), and (17). Dashed (n = 1), dotted (n = 2), and dash-dotted (n = 3) curves show exact results. The protocol k(t) is linear with k f /ki = 1/2. The coefficients α + 1 = 1/2, α − 1 = 2, α + 2 = 1/4, α − 2 = 1/3, α + 3 = 1/2 and α − 3 = 1/4 are chosen such that any moment approximated by the same friction coefficient collapses onto a single curve. In all cases the exact calculations agree with the approximation for large ∆t/τ (1) f .

FINITE INTEGRATION BOUNDS
The infinite integration bound on the friction assumes that correlations relax quickly relative to the protocol speed. For the first two moments, despite the finite integration bound generally yielding a more accurate approximation, the infinite bound considerably simplifies the approximation and allows for straightforward protocol optimization.
The effect of finite integration bounds is significant for two calculations: approximation of fourth-and higherorder moments of the excess work, and next-order approximations (including both leading-and next-to-leading-order contributions) for the moments. In the former, one must treat the bound as finite since the n-time covariance functions do not decay to zero for some subspace of large time arguments.
In more detail, consider the four-time covariance (kurtosis) δf i (0)δf j (t 2 )δf k (t 3 )δf (t 4 ) λ(t) . When any one of the four times {0, t 2 , t 2 , t 4 } significantly differs from the others, the conjugate forces decorrelate, and the kurtosis decays to zero; however, when t 3 ∼ t 4 t 2 ∼ 0, any variables separated by significant time decorrelate, and the kurtosis approaches δf i (0)δf j (t 2 ) λ(t) δf k (t 3 )δf (t 4 ) λ(t) . This limit represents a plane in the (t 1 , t 2 , t 3 ) parameter space where the kurtosis asymptotes to a finite value even for large time arguments. Integrating the above and multiplying by three (to account for permutations of the indices 1, 2, 3, 4) yields the integral four-time covariance in the limit t 3 ∼ t 4 t 2 ∼ 0, Since this is the only case that remains finite as t → ∞, the second line is the approximation (to highest order in t) for the integral four-time covariance. Parallel arguments hold for the higher-order moments. When approximating the average excess work (22) with both the Stokes' and supra-Stokes' tensors (23), finite integration bounds on the Stokes' friction may be necessary. For finite integration bounds on the Stokes' friction we replace (23) with For systems with weakly skewed conjugate-force fluctuation, the contribution from the finite bound on the first term can be comparable in magnitude to the contribution from the supra-Stokes' tensor. Figures 3 shows the first and second moments of the excess work for forward and reverse protocols, compared to different approximations. For the mean excess work, the supra-Stokes' tensor contributes with opposite sign from the Stokes' friction for forward (decreasing k) protocols and with same sign for reverse (increasing k). Since the Stokes'friction approximation overestimates the excess work in both cases (a,c), adding the supra-Stokes' tensor reduces the accuracy of the approximation for the reverse excess work. This effect is an artifact of the infinite integral bound; indeed, if the bound is kept finite (C (1) rather than ζ (1) ), the first integral autocovariance overestimates for forward protocols and underestimates for reverse protocols, and adding the supra-Stokes' tensor improves the approximation in both cases. Finally, we note that the Stokes' friction does not always overestimate for both forward and reverse protocols. For the second moment in (b,d), the Stokes' friction overestimates for forward protocols and underestimates for reverse protocols, and the supra-Stokes' tensor improves the approximation when either the Stokes' friction or first integral covariance are used.
Optimization of the (finite-bound) integral n-time covariance functions is not in general analytically tractable, but still could be done numerically.

FINSLER GEOMETRY
The addition of the supra-Stokes' tensor (7) comes at the cost of a more complex geometric structure. In contrast to Riemannian geometry, Finsler geometry is not restricted to a quadratic norm. In general, the inner products are not characterized by points, but rather by points and directions. Despite this, several useful concepts from Riemannian geometry (notably curvature, length, and geodesics) generalize. There are therefore standard procedures to find the geodesics despite the more complex landscapes induced by the generalized cubic Finsler metric (23).