Stochastic thermodynamics with SFI

SFI infers the two fields that define an overdamped stochastic system — the drift \(\mathbf{F}(\mathbf{x})\) and the diffusion \(\mathbf{D}(\mathbf{x})\). For a physicist these are more than a model of the dynamics: they carry the system’s thermodynamics. From a fitted model one can reconstruct an energy landscape, exhibit the probability currents that break detailed balance, and measure the rate at which the system dissipates free energy — the entropy production — directly from trajectory data, with error bars and no ground truth.

This page collects the concepts behind these measurements and routes to the executable demos. The worked example is Measuring dissipation: entropy production in a driven optical trap — a driven optical trap whose dissipation is invisible to every static observable.

Throughout, units are natural: \(k_B T\) is absorbed so that energies are in units of \(k_B T\), entropies in units of \(k_B\) (nats), and the mobility is unity so that \(D = k_B T\).

Equilibrium: forces from energies

A system is in equilibrium when its drift derives from a potential, \(\mathbf{F} = -\nabla U\) (constant \(D\); state-dependent diffusion modifies this relation — see Multiplicative noise: equilibrium, spurious drifts, blowtorches). Then the dynamics obeys detailed balance: every trajectory and its time-reverse are equally probable in the steady state, probability currents vanish, and the density is the Boltzmann weight

\[\rho(\mathbf{x}) \;\propto\; e^{-U(\mathbf{x})/D} .\]

Three SFI workflows are natural here:

  • Reconstruct the potential from a generic fit. Fit the force with infer_force_linear(), then integrate (or, for a linear force, symmetrize) the result. The Boltzmann weight of the reconstructed potential can be compared to the observed density — a strong self-consistency check.

  • Fit a gradient basis — still linear. Pick scalar functions \(u_i\) and use their (negative) gradients as the vector basis, \(\mathbf{b}_i = -\nabla u_i\) — one d_x() call on a scalar Basis builds it. The fit stays linear in the coefficients, so infer_force_linear() applies unchanged, the fitted force is exactly conservative, and the potential comes out linearly, \(U = \sum_i c_i u_i\), with error bars inherited from the coefficients. (Drop the constant function: its gradient is the zero feature.)

  • Parametrize by the energy. With the parametric estimator (infer_force()) and JAX automatic differentiation, the model can be written as \(\mathbf{F}_\theta = -\nabla U_\theta\) for any differentiable energy family — the fit returns the energy parameters directly, with error bars. See the energy-recovery section of Measuring dissipation: entropy production in a driven optical trap for a three-parameter example, and Neural-network energy landscape — Müller-Brown potential for a neural-network landscape (Müller–Brown) fitted through exactly this -U.d_x() pattern.

Tip

Fitting \(-\nabla U_\theta\) instead of a generic vector field builds the physics in: the fitted model is exactly conservative, with typically far fewer parameters — and what the data cannot express as a gradient shows up as residual, not as spurious force terms.

Breaking detailed balance: currents

Out of equilibrium, the steady state carries a non-zero probability current \(\mathbf{j} = \rho\,\mathbf{v}\), with mean local (phase-space) velocity

\[\mathbf{v}(\mathbf{x}) \;=\; \mathbf{F}(\mathbf{x}) - D\,\nabla \ln \rho(\mathbf{x}) \qquad\text{(constant } D\text{)}.\]

The velocity field is the order parameter of irreversibility: it vanishes identically in equilibrium and traces the circulation loops of a non-equilibrium steady state. Crucially, the density itself can be blind to it — in the driven optical trap of Measuring dissipation: entropy production in a driven optical trap, \(\rho\) remains exactly Boltzmann while \(\mathbf{v}\) circulates. For a strongly driven example with \(O(1)\) currents, see 2D limit cycle — nonlinear overdamped inference.

The observed-phase-space caveat. SFI sees currents only in the coordinates you record. Hidden degrees of freedom — an unobserved coordinate, a chemical state, fast solvent modes — can hide part (or all) of the circulation: projecting a non-equilibrium process onto fewer coordinates can only decrease the measured irreversibility. Every current-based quantity below is therefore a lower bound on the system’s dissipation, attached to the phase space you observe. This is a feature as much as a limitation: a positive measurement is decisive evidence of broken detailed balance, while a null one only bounds it.

Measuring entropy production

The dissipation of an overdamped steady state is the entropy production rate

\[\sigma \;=\; \bigl\langle \mathbf{v} \cdot D^{-1} \mathbf{v} \bigr\rangle \;\ge\; 0 ,\]

in units of \(k_B\) per unit time (heat dissipated into the bath divided by \(k_B T\)). SFI estimates it with compute_entropy_production() (after any linear force fit), by projecting the phase-space velocity onto the force basis — the estimator introduced in Frishman & Ronceray (2020):

\[\Delta\hat S \;=\; 2\, M_v^\top G^{-1} M_v , \qquad M_{v,a} = \sum_t \mathrm{d}\mathbf{X}_t \circ A^{-1} \mathbf{b}_a , \qquad A = 2\bar D ,\]

where \(\circ\) denotes the Stratonovich (mid-point) discretization. The estimate comes with two qualifiers, both reported:

  • a fluctuation bias \(\dot S_{\rm bias} = 2N_b/\tau_N\): with \(N_b\) basis functions and total observation time \(\tau_N\), even a perfectly reversible trajectory yields this much apparent dissipation from noise alone (reported, and also subtracted in the Sdot_debiased / DeltaS_debiased outputs);

  • a statistical error \(\delta\Delta S = \sqrt{2\,\Delta\hat S + (2N_b)^2}\).

The fluctuation bias is nothing exotic — it is the AIC complexity correction in thermodynamic clothing. \(\Delta\hat S\) is four times the plug-in information carried by the projected velocity field (see two integrals below), and a maximized log-likelihood — or information — exceeds its expectation by \(\tfrac12\) per fitted degree of freedom: the same counting that penalizes model complexity in sparsify_force() (criterion="AIC") removes \(2N_b\) from \(\Delta\hat S\). Use the raw value as the conservative plug-in, the debiased one for point estimates — near equilibrium it fluctuates around zero and may come out negative, which is honesty, not pathology.

Together, bias and error set a thermodynamic limit of detection: irreversibility becomes measurable only once the trajectory has dissipated a few \(k_B\) beyond the bias. Near equilibrium this is expensive — currents are linear in the driving, so \(\sigma\) is quadratic: a twice-fainter drive takes four times the data. The detection-limit sweep in Measuring dissipation: entropy production in a driven optical trap shows both effects quantitatively.

Two integrals, two meanings. The same trajectory functional read in the two discretizations measures two different things (see the observables of OverdampedProcess and Physics Reference):

\[I \;=\; \tfrac14 \textstyle\sum_t \mathrm{d}\mathbf{X}_t \cdot D^{-1} \mathbf{F}(\mathbf{x}_t) \qquad\text{(Itô)}\]

is the information the trajectory carries about the force field — the quantity that controls SFI’s error bars (force_information) — while the Stratonovich reading

\[\Delta S \;=\; \textstyle\sum_t \mathrm{d}\mathbf{X}_t \circ D^{-1} \mathbf{F}(\mathbf{x})\]

is the log-ratio of forward to time-reversed path probabilities: the entropy production. The two are tightly related — the entropy production rate is four times the information rate carried by the velocity field, \(\sigma = 4\,\dot I_v\) — so measuring dissipation is exactly accumulating information about the irreversible part of the dynamics.

Practical caveats.

  • The estimator assumes constant (or averaged) diffusion; a warning is emitted when a state-dependent \(D(\mathbf{x})\) field was fitted (see Multiplicative noise: equilibrium, spurious drifts, blowtorches).

  • Measurement (localization) noise is benign at leading order — a parity gift. The Stratonovich moments are odd under time reversal, while independent localization noise is even: the naive \(O(\Lambda/\Delta t)\) noise–increment correlation cancels exactly, telescoping to a boundary term. (This is special to the entropy estimator: the Itô force moments enjoy no such cancellation.) What survives is second order — an \(O(\Lambda)\) smoothing of the basis functions and a noise-inflated variance that raises the effective fluctuation bias — and is not included in the error bar. The exception is noise correlated with the motion: motion blur from finite exposure breaks the parity argument and does bias the estimate. Denoise or subsample for precision, not out of fear of a diverging bias (see Measurement noise and coarse sampling).

Underdamped systems: dissipation from positions only (experimental)

Warning

Experimental. The underdamped estimator is new and its error theory is provisional: the bias and error formulas mirror the overdamped derivation and are pinned empirically (a \(\gamma\Delta t\) calibration sweep), not derived from first principles. Expect the outputs and error bars to be refined.

With inertia the state is \((\mathbf{x}, \mathbf{v})\) and time reversal flips the velocity: \(\mathbf{x}\) is even, \(\mathbf{v}\) is odd. The fitted acceleration field then splits by parity,

\[F^{\pm}(\mathbf{x},\mathbf{v}) = \tfrac12\bigl[F(\mathbf{x},\mathbf{v}) \pm F(\mathbf{x},-\mathbf{v})\bigr] ,\]

into a reversible part \(F^{+}\) (conservative and confining forces) and an irreversible part \(F^{-}\) (friction and velocity-odd driving) — the decomposition that organises stochastic thermodynamics with odd variables (Spinney & Ford 2012; Lee, Kwon & Park 2013). time_reversal_split() returns the two parts of a fitted model, and compute_entropy_production() evaluates the Stratonovich log path-probability ratio

\[\Delta \hat S = \sum_t \bigl(\mathrm{d}\hat v_t \circ - \hat F^{+}\mathrm{d}t\bigr)^\top \bar D^{-1} \hat F^{-} ,\]

with the velocities \(\hat v\) reconstructed from positions (the ULI kinematics) — no velocity measurements are needed. In a steady state its mean is the entropy production rate whenever \(\nabla_v \cdot F^{+} = 0\); the heat dissipated into the bath follows as \(\dot Q = T\,\sigma\) with \(T = \bar D/\hat\gamma\) (Sekimoto 1998). For any linear underdamped model the exact rate is available in closed form from the phase-space Ornstein–Uhlenbeck formula with the parity-aware irreversible drift (see Qian 2001; Godrèche & Luck 2019) — the validation route used by SFI’s test suite.

Underdamped-specific caveats:

  • Resolution prerequisite. Velocity reconstruction requires \(\Delta t \lesssim \tau_v\) (the velocity correlation time); a warning is emitted when \(\hat\gamma\,\Delta t > 0.5\). Beyond that the measurement degrades toward a coarse-grained lower bound.

  • Fluctuation scale. Friction makes the odd sector large even at equilibrium: the log ratio’s contributions cancel in the mean but not in variance, so error bars are set by the odd quadratic \(\langle F^-\!\cdot \bar D^{-1} F^-\rangle\) (reported, with an empirical block-variance estimate).

  • Same-sample plug-in bias. Evaluating the functional with coefficients fitted on the same trajectory produces a positive \(O(1/\tau_N)\) bias; it is removed by cross-fitting (fit on one half, evaluate on the other — the coefficients argument), and it decays with trajectory length.

  • Overdamped limits are singular. With temperature gradients or state-dependent friction, the small-inertia limit carries an entropic anomaly absent from the overdamped description (Celani et al. 2012): analysing inertial data with an overdamped model can be systematically wrong — infer at the underdamped level when inertia is resolvable (the classify_dynamics() verdict helps decide).

Multiplicative noise: equilibrium, spurious drifts, blowtorches

When the noise amplitude depends on the state — multiplicative noise, \(D = D(\mathbf{x})\) — equilibrium intuition breaks in a subtler way. The Itô drift that SFI infers is then not the mechanical force: for a system with mobility \(\mu(\mathbf{x})\) and physical force \(\mathbf{f}\),

\[\mathbf{F}_{\rm It\hat o}(\mathbf{x}) \;=\; \mu(\mathbf{x})\,\mathbf{f}(\mathbf{x}) \;+\; \nabla\!\cdot\! \mathbf{D}(\mathbf{x}) ,\]

and the divergence term — the spurious drift — pushes probability around even with no force at all. Interpreting an inferred drift thermodynamically therefore requires the diffusion field, not just its average — infer it with infer_diffusion_linear() or infer_diffusion() (see Anisotropic diffusion tensor field for tensor-valued fields).

Equilibrium with multiplicative noise. State-dependent noise does not by itself break detailed balance. When the diffusion field reflects a state-dependent mobility at uniform temperature — fluctuation–dissipation, \(\mathbf{D}(\mathbf{x}) = \mu(\mathbf{x})\,k_BT\) — a conservative force \(\mathbf{f} = -\nabla U\) gives (in our \(k_BT = 1\) units)

\[\mathbf{F}_{\rm It\hat o}(\mathbf{x}) \;=\; -\,\mathbf{D}(\mathbf{x})\,\nabla U(\mathbf{x}) \;+\; \nabla\!\cdot\!\mathbf{D}(\mathbf{x}) ,\]

and the stationary density remains exactly Boltzmann, \(\rho \propto e^{-U}\) — the spurious drift is precisely what keeps it so (plug \(\rho = e^{-U}\) into the zero-current condition \(\mathbf{F}\rho = \nabla\!\cdot\!(\mathbf{D}\rho)\) to check). Read as an inference statement, this is a sharp equilibrium test: a fitted pair \((\hat{\mathbf{F}}, \hat{\mathbf{D}})\) is compatible with equilibrium if and only if

\[\boldsymbol{\Phi}(\mathbf{x}) \;=\; \hat{\mathbf{D}}^{-1}\bigl(\hat{\mathbf{F}} - \nabla\!\cdot\!\hat{\mathbf{D}}\bigr)\]

is a gradient field, \(\boldsymbol{\Phi} = -\nabla U\); any rotational part of \(\boldsymbol{\Phi}\) is a genuine detailed-balance violation in the observed coordinates.

Joint \((U, D)\) inference in practice. Recovering the energy landscape under multiplicative noise therefore takes both fields, and SFI composes them rather than adding a dedicated mode:

  1. fit the diffusion field (infer_diffusion_linear() or the noise-robust infer_diffusion());

  2. fit the drift (infer_force_linear() or infer_force());

  3. form \(\boldsymbol{\Phi}\) from the two fits — the divergence \(\nabla\!\cdot\!\hat{\mathbf{D}}\) is one d_x() call away — and test it for gradientness / line-integrate it for \(U\).

Alternatively, build the equilibrium ansatz in: with the fitted \(\hat{\mathbf{D}}(\mathbf{x})\) frozen, write the drift family

\[\mathbf{F}_\theta \;=\; -\,\hat{\mathbf{D}}(\mathbf{x})\,\nabla U_\theta(\mathbf{x}) \;+\; \nabla\!\cdot\!\hat{\mathbf{D}}(\mathbf{x})\]

(autodiff supplies both derivative terms) and fit \(U_\theta\) directly with infer_force() — the state-dependent-\(D\) generalization of the gradient-basis and autodiff-energy workflows of the equilibrium section. Misfit of this constrained family relative to a free drift fit is then itself evidence of broken detailed balance.

Blowtorches. If instead the noise field reflects a genuine temperature gradient, no potential \(U\) exists at all: the steady state is not the naive Boltzmann weight \(e^{-U/\bar D}\), and occupancy is re-shaped by the temperature landscape — the Landauer blowtorch effect, with heat flowing even in a static system. See Multiplicative noise — the Landauer blowtorch for the blowtorch worked example. Entropy production with a state-dependent \(\mathbf{D}\) is where the current estimator stops: it contracts with the constant \(\bar D\) and warns (the natural extension — a trajectory integral against \(\mathbf{D}(\mathbf{x})^{-1}\) — is noted in compute_entropy_production()).

Which tool for which question

Question

Tool / demo

Is this system out of equilibrium?

compute_entropy_production() — positive beyond bias and error bar is decisive; Measuring dissipation: entropy production in a driven optical trap

What drives it — which terms break detailed balance?

sparsify_force() (PASTIS) + symmetric/antisymmetric decomposition of the fit

What is the energy landscape?

Parametric energy fit via autodiff (infer_force()); Neural-network energy landscape — Müller-Brown potential

Where does probability circulate?

stream_field() on the non-gradient part of the fit; 2D limit cycle — nonlinear overdamped inference

Is the noise state-dependent, and what does it do to occupancy?

infer_diffusion_linear(); Multiplicative noise — the Landauer blowtorch

Equilibrium despite multiplicative noise?

Gradientness of \(\hat{\mathbf{D}}^{-1}(\hat{\mathbf{F}} - \nabla\!\cdot\!\hat{\mathbf{D}})\) — see Multiplicative noise: equilibrium, spurious drifts, blowtorches

Driven, or inertial?

classify_dynamics(); Overdamped or underdamped? Classifying dynamics from data

References

Method:

Frishman, A. & Ronceray, P., Learning force fields from stochastic trajectories, Physical Review X 10, 021009 (2020). DOI: 10.1103/PhysRevX.10.021009 — introduces both the inference framework and the entropy-production estimator used here.

Brückner, D. B., Ronceray, P. & Broedersz, C. P., Inferring the dynamics of underdamped stochastic systems, Physical Review Letters 125, 058103 (2020). DOI: 10.1103/PhysRevLett.125.058103 — the velocity-reconstruction engine behind the underdamped estimator.

Stochastic thermodynamics:

Seifert, U., Stochastic thermodynamics, fluctuation theorems and molecular machines, Reports on Progress in Physics 75, 126001 (2012). DOI: 10.1088/0034-4885/75/12/126001

Sekimoto, K., Langevin equation and thermodynamics, Progress of Theoretical Physics Supplement 130, 17 (1998). DOI: 10.1143/PTPS.130.17

Odd-parity (underdamped) entropy production:

Spinney, R. E. & Ford, I. J., Nonequilibrium thermodynamics of stochastic systems with odd and even variables, Physical Review Letters 108, 170603 (2012). DOI: 10.1103/PhysRevLett.108.170603; and Entropy production in full phase space for continuous stochastic dynamics, Physical Review E 85, 051113 (2012). DOI: 10.1103/PhysRevE.85.051113

Lee, H. K., Kwon, C. & Park, H., Fluctuation theorems and entropy production with odd-parity variables, Physical Review Letters 110, 050602 (2013). DOI: 10.1103/PhysRevLett.110.050602

Celani, A., Bo, S., Eichhorn, R. & Aurell, E., Anomalous thermodynamics at the microscale, Physical Review Letters 109, 260603 (2012). DOI: 10.1103/PhysRevLett.109.260603

Exact linear (phase-space Ornstein–Uhlenbeck) benchmarks:

Qian, H., Mathematical formalism for isothermal linear irreversibility, Proceedings of the Royal Society A 457, 2645 (2001). DOI: 10.1098/rspa.2001.0811

Godrèche, C. & Luck, J.-M., Characterising the nonequilibrium stationary states of Ornstein–Uhlenbeck processes, Journal of Physics A 52, 035002 (2019). DOI: 10.1088/1751-8121/aaf190

Broken detailed balance in living and driven systems:

Battle, C. et al., Broken detailed balance at mesoscopic scales in active biological systems, Science 352, 604 (2016). DOI: 10.1126/science.aac8167

Gnesotto, F. S., Mura, F., Gladrow, J. & Broedersz, C. P., Broken detailed balance and non-equilibrium dynamics in living systems: a review, Reports on Progress in Physics 81, 066601 (2018). DOI: 10.1088/1361-6633/aab3ed

Non-conservative forces in optical traps:

Roichman, Y., Sun, B., Stolarski, A. & Grier, D. G., Influence of nonconservative optical forces on the dynamics of optically trapped colloidal spheres: the fountain of probability, Physical Review Letters 101, 128301 (2008). DOI: 10.1103/PhysRevLett.101.128301

Wu, P., Huang, R., Tischer, C., Jonas, A. & Florin, E.-L., Direct measurement of the nonconservative force field generated by optical tweezers, Physical Review Letters 103, 108101 (2009). DOI: 10.1103/PhysRevLett.103.108101