Optical trap Langevin calibration (source code)

= Optical trap Langevin calibration
{c}
{title2=$\zeta\dot x=-kx+\xi(t),\quad f_c=k/(2\pi\zeta)$}

Add the random thermal <force> to the deterministic <optical trap relaxation time> model. In one stationary coordinate of a <harmonic optical trap>, the overdamped <Langevin equation> is
$$
\zeta\dot x=-kx+\xi(t),\qquad
\langle\xi(t)\rangle=0,\qquad
\langle\xi(t)\xi(t')\rangle=2k_BT\zeta\,\delta(t-t').
$$
Here $\xi$ is <Gaussian white noise> understood through its integrals, $\delta$ is the <Dirac delta function>, and $\zeta$ is the <linear friction coefficient>. The covariance coefficient has units $\mathrm{N^2\,s}$; multiplying by $\delta(t-t')$, of units $\mathrm{s^{-1}}$, gives squared <force>. This <fluctuation-dissipation relation for a Langevin particle> makes the stationary <variance> $\sigma^2=k_BT/k$, consistent with the <equipartition theorem>. The solution is an <Ornstein-Uhlenbeck process> with <autocovariance function>
$$
C_x(t)=\sigma^2e^{-|t|/\tau},\qquad
\tau=\frac{\zeta}{k}.
$$
Indeed, at positive lag the noise has zero correlation with the initial position, so $\zeta C_x'=-kC_x$; its initial value is $\sigma^2$. This supplies both a decay time and an amplitude, making a stronger calibration test than a position histogram alone.

Specify the <Fourier transform> convention before using a <power spectral density>. For <frequency> $f$ in <hertz>, define $S_x^{(2)}(f)=\int_{-\infty}^{\infty}C_x(t)e^{-2\pi ift}\,dt$ and the one-sided density $S_x^{(1)}(f)=2S_x^{(2)}(f)$ for $f>0$. Integration of the exponential gives
$$
S_x^{(1)}(f)=\frac{4k_BT\zeta}{k^2+(2\pi f\zeta)^2}
=\frac{k_BT}{\pi^2\zeta(f^2+f_c^2)},\qquad
f_c=\frac{k}{2\pi\zeta},\qquad
\int_0^\infty S_x^{(1)}(f)\,df=\frac{k_BT}{k}.
$$
The position <power spectral density> has units $\mathrm{m^2/Hz}$; its integral is a <variance>, not a <standard deviation>. At $f=f_c$ it is half its zero-frequency plateau. Angular <frequency> conventions require a corresponding $2\pi$ conversion in both axis and integration measure. If $x=\beta(V-V_0)$ with position sensitivity $\beta$ in $\mathrm{m/V}$, then $S_V^{(1)}=S_x^{(1)}/\beta^2$ in $\mathrm{V^2/Hz}$ before independent detector noise is added.

For an original invented calibration example, choose $a=0.500\,\mathrm{\mu m}$, <dynamic viscosity> $\eta=1.00\times10^{-3}\,\mathrm{Pa\,s}$, $T=300\,\mathrm K$, $k=0.0200\,\mathrm{pN/nm}=2.00\times10^{-5}\,\mathrm{N/m}$, and $\beta=100\,\mathrm{nm/V}$. The unbounded no-slip <Stokes drag law> gives
$$
\zeta=9.42478\times10^{-9}\,\mathrm{N\,s/m},\quad
\tau=0.471239\,\mathrm{ms},\quad
f_c=337.737\,\mathrm{Hz},\quad
\sigma=14.391\,\mathrm{nm}.
$$
The corresponding voltage <standard deviation> is $0.14391\,\mathrm V$, and the position <power spectral density> plateau is $0.390369\,\mathrm{nm^2/Hz}$, or $3.90369\times10^{-5}\,\mathrm{V^2/Hz}$ in the detector. Recovering both $k=k_BT/\sigma^2$ and $k=2\pi\zeta f_c$ is an internal model check; actual agreement also tests the assumed <temperature>, drag and detector response.

\Image[/optical-tweezers-synthetic-calibration.png]
{title=Synthetic <harmonic optical trap> fluctuations, <power spectral density> and <optical trap drag calibration>}
{description=Original model calculations and seeded synthetic <Ornstein-Uhlenbeck process> samples; explicitly not measured data. The $10\,\mathrm s$ record is sampled at $20\,\mathrm{kHz}$ with no detector noise or exposure averaging. The dashed spectrum includes sampling aliases; the red continuous spectrum integrates to the exact <equipartition theorem> variance. The drag line is an exact steady-model prediction.}
{source=optical-tweezers-synthetic-calibration.py}

The generator uses the exact stationary update $x_{j+1}=\rho x_j+\sigma\sqrt{1-\rho^2}\,Z_j$, with $\rho=e^{-\Delta t/\tau}$ and independent standard <normal random variables> $Z_j$. For a demeaned block of $N$ samples and its unnormalized <discrete Fourier transform> $X_j$, the two-sided <periodogram> is $\Delta t|X_j|^2/N$; doubling positive bins except the zero and Nyquist endpoints gives a one-sided density. Summing those bins times $\Delta f=1/(N\Delta t)$ returns the block mean square by <discrete Parseval identity>. Independent detector noise, finite exposure and <aliasing> would alter a measured spectrum and must enter a fit before interpreting $f_c$.