Harmonic Oscillator in Heat Bath: Exact simulation of time-lapse-recorded data, exact analytical benchmark statistics

The stochastic dynamics of the damped harmonic oscillator in a heat bath is simulated with an algorithm that is exact for time steps of arbitrary size. Exact analytical results are given for correlation functions and power spectra in the form they acquire when computed from experimental time-lapse recordings. Three applications are discussed: (i) Effects of finite sampling-rate and -time, described exactly here, are similar for other stochastic dynamical systems-e.g. motile micro-organisms and their time-lapse recorded trajectories. (ii) The same statistics is satisfied by any experimental system to the extent it is interpreted as a damped harmonic oscillator at finite temperature-such as an AFM cantilever. (iii) Three other models of fundamental interest are limiting cases of the damped harmonic oscillator at finite temperature; it consequently bridges their differences and describes effects of finite sampling rate and sampling time for these models as well. Finally, we give a brief discussion of nondimensionalization.

💡 Research Summary

The paper presents an exact numerical scheme for simulating the stochastic dynamics of a damped harmonic oscillator immersed in a thermal bath, together with closed‑form analytical expressions for the statistical measures that are directly comparable to experimentally recorded time‑lapse data. Starting from the Langevin equation (m\ddot x+\gamma\dot x+\kappa x=F_{\rm therm}(t)) with white‑noise thermal force, the authors rewrite the second‑order equation as a pair of first‑order linear stochastic differential equations. By introducing the 2×2 matrix (M) that encodes the deterministic dynamics, they derive the exact discrete‑time propagator (e^{-M\Delta t}) for any time step (\Delta t). This propagator can be expressed analytically in terms of cosine and sine functions of the damped natural frequency (\omega=\sqrt{\kappa/m-\gamma^{2}/(4m^{2})}) and the relaxation time (\tau=m/\gamma).

The update rule for the state vector ((x_{j},v_{j})) reads
((x_{j+1},v_{j+1})^{T}=e^{-M\Delta t}(x_{j},v_{j})^{T}+(\Delta x_{j},\Delta v_{j})^{T}),
where the stochastic increments (\Delta x_{j}) and (\Delta v_{j}) are jointly Gaussian with a covariance matrix that the authors compute explicitly. The matrix elements (\sigma_{xx},\sigma_{vv},\sigma_{xv}) depend on (\Delta t), (\omega) and (\tau); in the critically damped limit ((\omega=0)) they reduce to simple exponential forms.

Using this exact discretization, the authors obtain analytical expressions for the mean‑squared displacement (MSD), the position and velocity autocorrelation functions, and the cross‑correlation. The MSD is shown to be ( {\rm MSD}(t)=2\langle x^{2}\rangle-2\langle x(t)x(0)\rangle), which is simply the continuous‑time result evaluated at discrete sampling times. The autocorrelation functions decay as (\exp(-t/2\tau)) while oscillating at frequency (\omega) when the system is under‑damped; for over‑damped conditions the oscillatory part disappears and the decay becomes a sum of two exponentials. The paper also discusses the practical issue of estimating velocity from finite‑difference positions (“secant velocity”). It demonstrates that the secant‑velocity autocorrelation deviates markedly from the true instantaneous‑velocity correlation at short lags, a point that experimentalists must consider when interpreting high‑frequency data.

Power spectral densities (PSDs) for position and velocity are derived in a normalization that keeps total power constant as the measurement time grows. The position PSD is
(P_{x}(f)=\frac{D}{2\pi^{2}}\frac{(2\pi\tau)^{2}}{(f^{2}-f_{0}^{2})^{2}+f^{2}/\tau^{2}}),
with (f_{0}=\omega_{0}/2\pi) the undamped resonance frequency and (D=k_{B}T/\gamma) the diffusion coefficient. The velocity PSD follows from the relation (\tilde v= i2\pi f \tilde x). Simulations confirm that the theoretical PSDs match the numerically generated spectra across under‑, critically, and over‑damped regimes, showing the characteristic Lorentzian peak for under‑damped motion and a low‑frequency (1/f) tail for over‑damped motion.

The authors further show that three well‑known limiting cases—(i) vanishing mass (over‑damped Brownian particle in a trap), (ii) vanishing spring constant (free Brownian particle with inertia), and (iii) both mass and spring constant zero (Einstein‑Smoluchowski diffusion)—are obtained by taking appropriate limits of the general solution. This unifying perspective clarifies how finite‑sampling effects propagate across these models.

A brief discussion of nondimensionalization follows, introducing natural scales for time ((\tau)), length ((\sqrt{k_{B}T/\kappa})), and energy ((k_{B}T)). By rescaling, the number of independent parameters is reduced, facilitating comparison between different experimental systems such as AFM cantilevers, optical traps, and cell‑migration trajectories.

In conclusion, the paper provides a rigorously exact algorithm for simulating damped harmonic motion at any time step, together with analytical benchmark formulas for MSD, autocorrelations, and PSDs that incorporate finite‑sampling and finite‑measurement‑time effects. These results are directly applicable to a broad range of experimental contexts where the system can be modeled as a harmonic oscillator coupled to a thermal bath, and they lay the groundwork for extending exact treatments to more complex, possibly nonlinear, stochastic dynamics.