Introduction #

Harmonic oscillators are good model for a very wide class of systems. Interaction with an external environment will introduce thermal noise. In this post we will quantify how much thermal noise can affect a harmonic oscillator over a given time period. This is important for experiments where a harmonic oscillator is used to make a precise measurement, or to investigate quantum effects. In such cases we must make sure thermal fluctuations are smaller than whatever we are trying to observe.

This article will be a bit rough, it's mostly me writing down a calculation for later reference than exposition. If you have any questions or if anything is unclear, feel free to send me an email.

Model #

A classical harmonic oscillator driven by thermal noise has equation of motion
[eqOscillatorEOM]: $\ddot{x}(t)+\Gamma\dot{x}(t)+\Omega^2 x(t)=\sqrt{\frac{2\Gamma k_BT}{m}}\xi(t).$
Here $x$ is the position of the oscillator, $m$ its mass, $\Gamma$ its damping rate, and $\Omega$ its resonant frequency. The right-hand side describes thermal flutuations. $k_B$ is Boltzmann's constant, $T$ is temperature, and $m$ is the mass of the oscillator.

Randomness comes from $\xi(t)$, which is a white noise process. Intuitively, function whose value changes randomly with time. It is defined such that over an infinitesimal interval $\mathrm{d}t$, $\xi(t)\mathrm{d}t$ is Gaussian distributed with mean $0$ and variance $\mathrm{d}t$ [J10]. To work with $\xi(t)$, we only need two relations. The first is that it is zero mean:
[eqXiMean]: $\langle\xi(t)\rangle=0,$
where $\langle\cdot\rangle$ is the expected value. Secondly, we have the expected correlation between the noise at different times:
[eqXiCorrelator]: $\langle\xi(t)\xi(t')\rangle=\delta(t-t').$
This says that noise at different times is completely uncorrelated. In other words, knowing the value of $\xi$ at some time $t$, gives you exactly no knowledge about the value of $\xi$ a a different $t'$. Since $t,t'$ could be arbitrarily close together, this requires $\xi$ to be varying infinitely fast. This is an approximation, but it will be valid so long as the thermal fluctuations are much faster than our system frequencies of $\Omega$ and $\Gamma$.

Note that the damping rate $\Gamma$ appears on both sides of [eqOscillatorEOM]. This is a manifestation of the fluctuation-dissipation theorem, which states that any mechanism that can dissipate energy from the system, will also add fluctuations. For example if a particle is moving through water, the particle's motion will be damped due to collision with water molecules. At the same time, random motion of the water molecules will introduce fluctuations via Brownian motion. A consequence of this is that the less damped your system is, the less fluctuations the environment will induce. It is for this reason that in quantum experiments, we try and reduce as much as possible the damping of our system.

Approximate equations of motion #

To understand thermal noise in the oscillator, we need to solve [eqOscillatorEOM]. In this section we will use our intuition for how the solutions should behave, to derive approximate equations of motion that we will be able to solve.

We start by making an ansatz. The solution will be oscillatory, so the dynamics should look like
$x(t)=A(t)\cos\left(\Omega t+\phi(t)\right),$
where $A(t)$, $\phi(t)$ are a time-varying amplitude and phase. Substituting these into [eqOscillatorEOM] gives
$\begin{aligned}&\left[\ddot{A}+\Gamma\dot{A}-A\dot{\phi}(2\Omega+\dot{\phi})\right]\cos\left(\Omega t+\phi\right)\\&-\left[(\Gamma A+2\dot{A})(\Omega+\dot{\phi})+A\ddot{\phi}\right]\sin\left(\Omega t+\phi\right) \\ &= \sqrt{\frac{2\Gamma k_BT}{m}}\xi(t).\end{aligned}$

Let us consider the right-hand side. We can decompose $\xi(t)$ into quadratures as
[eqXiQuadratures]: $\xi(t)=\frac{1}{\sqrt{2}}\left(\xi_1(t)\sin(\Omega t+\phi)+\xi_2(t)\cos(\Omega t+\phi)\right),$
where $\xi_1,\xi_2$ are independent white noise sources:
$\langle\xi_1(t)\xi_2(t')\rangle=0,$
$\langle\xi_1(t)\rangle=\langle\xi_2(t)\rangle=0,$
$\langle\xi_1(t)\xi_1(t')\rangle=\langle\xi_2(t)\xi_2(t')\rangle=\delta(t-t').$
We can easily verify that both sides of [eqXiQuadratures] have the same first and second moments, and in the following we will only consider the first two moments.

With this we arrive at two coupled equations for the amplitude and phase:
$\ddot{A}(t)+\Gamma\dot{A}(t)-A\dot{\phi}(t)\left(2\Omega+\dot{\phi}(t)\right)=\sqrt{\frac{\Gamma k_BT}{m}}\xi_1(t),$
$\left(\Gamma A(t)+2\dot{A}(t)\right)\left(\Omega+\dot{\phi}(t)\right)+A(t)\ddot{\phi}(t) =\sqrt{\frac{\Gamma k_BT}{m}}\xi_2(t).$

To proceed further we will make some approximations:

1. The thermal relaxation time is slower than the oscillator's frequency: $\Omega\gg\Gamma$. In terms of the mechanical quality factor $Q=\frac{\Omega}{\Gamma}$, this approximation means that we are in the high-Q regime $Q\gg 1$.
2. The envelope $A(t)$ and phase $\phi(t)$ vary on the timescale of $\Gamma$. Then $\dot{A}\sim\mathcal{O}(\Gamma)$ and $\dot{\phi}\sim\mathcal{O}(\Gamma)$ (see big O notation).

Let us first consider the coefficient of cosine. The order of each term on the left-hand side is
$\underbrace{\ddot{A}}_{\mathcal{O}(\Gamma^2)}+\underbrace{\Gamma\dot{A}}_{\mathcal{O}(\Gamma^2)}-\underbrace{2\Omega\dot{\phi}A}_{\mathcal{O}(\Omega\Gamma )}-\underbrace{A\dot{\phi}^2}_{\mathcal{O}(\Gamma^2)}$
The third term is higher order than all others. Keeping just this gives us
[eqPhiDot]: $\dot{\phi}(t)=-\frac{1}{2A}\sqrt{\Gamma \frac{k_BT}{m\Omega^2 }}\xi_1(t).$

Now consider the coefficient of sine. The terms have order
$\underbrace{\Omega\Gamma A}_{\mathcal{O}(\Omega\Gamma)}+\underbrace{\Gamma A\dot{\phi}}_{\mathcal{O}(\Gamma^2)}+\underbrace{2\Omega\dot{A}}_{\mathcal{O}(\Omega\Gamma)}+\underbrace{2\dot{A}\dot{\phi}}_{\mathcal{O}(\Gamma^2)}+\underbrace{A\ddot{\phi}}_{\mathcal{O}(\Gamma^2)}.$
Here there are two terms of equal highest order, which gives us
[eqADot]: $\dot{A}(t)=-\frac{\Gamma}{2} A(t)-\frac{1}{2}\sqrt{\Gamma \frac{k_BT}{m\Omega^2}}\xi_2(t).$

Thus we have found two first-order differential equations, which describe how the amplitude and phase drift due to thermal noise.

Amplitude drift #

We are now in a position to quantify the thermal perturbation to an oscillator. Suppose we observe the system for a time interval $\tau$. Over this time, thermal noise will cause the amplitude to drift by some $\Delta A$. In this section we will study how big this drift is.

Integrating [eqADot] gives us an explicit expression for the change in amplitude:
$\begin{aligned}\Delta A&=\int_0^{\tau}e^{-(\tau -t)\Gamma/2}\frac{1}{2}\sqrt{\Gamma\frac{k_B T}{m\Omega^2}}\xi_1(t)\,\mathrm{d}t, \\ &=\frac{1}{2}\sqrt{\Gamma\frac{k_B T}{m\Omega^2}}\int_0^{\tau}e^{-(\tau-t)\Gamma/2}\xi_1(t)\mathrm{d}t.\end{aligned}$
Since $\xi_1$ is stochastic, this is a random variable. We can use the properties of white noise to understand the distributoin of $\Delta A$.

Let's define
$Z_A(\tau)=\int_0^{\tau}e^{-(\tau-t)\Gamma/2}\xi_A(t)\mathrm{d}t.$
Then
$\Delta A=\frac{1}{2}\sqrt{\Gamma\frac{k_B T}{m\Omega^2}}Z_A(\tau).$
From [eqXiMean] this has zero mean
$\langle Z_A(\tau)\rangle=\int_0^{\tau}e^{-(\tau-t)\Gamma/2}\langle\xi_A(t)\rangle\mathrm{d}t=0,$
so the average change in amplitude is zero: $\langle\Delta A\rangle=0$.

To find the variance we use [eqXiCorrelator]
$\begin{aligned}\langle Z_A(\tau)^2\rangle &=\left\langle\int_0^{\tau}e^{-(\tau-t_1)\Gamma/2}\xi_1(t_1)\mathrm{d}t_1,\int_0^{\tau}e^{-(\tau-t_2)\Gamma/2}\xi_1(t_2)\mathrm{d}t_2\right\rangle \\ &= \iint_0^{\tau}e^{-(\tau-t_1)\Gamma/2}e^{-(\tau-t_2)\Gamma/2}\left\langle\xi_1(t_1)\xi_1(t_2)\right\rangle\mathrm{d}t_1\mathrm{d}t_2, \\ &=\int_0^{\tau}e^{-(\tau-t_1)\Gamma}\mathrm{d}t_1, \\ &=\frac{1-e^{-\Gamma\tau}}{\Gamma}.\end{aligned}$
From this we obtain
$\sqrt{\langle\Delta A^2\rangle}=\frac{1}{2}\sqrt{\Gamma\frac{k_B T}{m\Omega^2}}\sqrt{\frac{1-e^{-\Gamma\tau}}{\Gamma}}.$
This is an explicit expression for the magnitude of thermal drift we can expect over a time interval $\tau$.

Let's first consider the long-time regime $\Gamma\tau\rightarrow\infty$. In this limit we have
$e^{-\Gamma\tau}\approx 0,$
which gives
$\sqrt{\langle\Delta A^2\rangle}\rvert_{\tau\gg 1/\Gamma}=\frac{1}{2}\sqrt{\frac{k_BT}{m\Omega^2}}.$
Thus at long times the drift is independent of $\tau$, depending only on the temperature, and mass and oscillation frequency.

Now consider a short time $\Gamma\tau\rightarrow 0$. Then we can approximate
$e^{-\Gamma\tau}\approx 1-\Gamma\tau,$
which gives
$\sqrt{\langle\Delta A^2\rangle}\rvert_{\tau\ll 1/\Gamma}=\frac{1}{2}\sqrt{\Gamma\tau\frac{k_BT}{m\Omega^2}}.$
So over short times, the amplitude drifts proportional to $\sqrt{\tau}$.

Phase drift #

Now let's study the drift in phase. Integrating [eqPhiDot] gives
$\begin{aligned}\Delta\phi &=\int_0^{\tau}\left(-\frac{1}{2A(t)}\sqrt{\Gamma\frac{k_BT}{m\Omega^2}}\right)\xi_2(t)\,\mathrm{d}t, \\ &= -\frac{1}{2}\sqrt{\Gamma\frac{k_B T}{m\Omega^2}}\int_0^{\tau}\frac{1}{A(t)}\xi_2(t)\mathrm{d}t.\end{aligned}$
Here $A(t)$ is a random variable which is independent of $\xi_2$ (since $A$ is driven by $\xi_1$, an independent noise source).

As before the expected phase drift is zero.
$\langle\Delta\phi\rangle=0$.
To find the variance let us define
$Z_{\phi}(\tau)=\int_0^{\tau}\frac{1}{A(t)}\xi_2(t)\mathrm{d}t.$
Then [eqXiCorrelator] gives
$\begin{aligned}\langle Z(\tau)^2\rangle &= \int_0^{\tau}\left\langle\frac{1}{A(t)^2}\right\rangle\mathrm{d}t.\end{aligned}$
Evaluating this is tricky.

Let's just consider small $\tau$. We can write $A=A_0+\Delta A$, where $A_0$ is the amplitude at the beginning of our time interval. If $\tau$ is small, then $\Delta A\ll A_0$, and we can Taylor expand:
$\begin{aligned}\frac{1}{A^2} &= \frac{1}{A_0^2+2A_0\Delta A+\Delta A^2}, \\ &=\frac{1}{A_0^2}\frac{1}{1+2\Delta A/A_0+\Delta A^2/A_0^2}, \\ &\approx\frac{1}{A_0^2}\frac{1}{1+2\Delta A/A_0}, \\ &\approx \frac{1}{A_0^2}\left(1-2\frac{\Delta A}{A_0}\right) \\ &\approx \frac{1}{A_0^2}. \end{aligned}$
$\left\langle\frac{1}{A^2}\right\rangle\approx \frac{1}{A_0^2}.$
Then at small times, $\Delta\phi$ will increase with $\sqrt{\tau}$:
$\sqrt{\langle\Delta\phi^2\rangle}\rvert_{\tau\ll 1/\Gamma} = \frac{1}{2A_0}\sqrt{\Gamma\tau\frac{k_B T}{m\Omega^2}}.$

Conclusion #

We have found that for time intervals $\tau$ small compared with the thermal relaxation time $1/\Gamma$, both the amplitude and phase of a harmonic oscillator will experience drift proportional to $\sqrt{\tau}$. At large times, the amplitude variance approaches a constant value.

This was all valid in the limit where the oscillation frequency $\Omega$ is fast compared with the termal relaxation time. This is a common regime in sensing experiments, where we generally suppress $\Gamma$ as much as possible.

References #

[J10] Kurt Jacobs (2010) Stochastic processes for physicists: understanding noisy systems Cambridge University Press

---

• Previous: Euler-Lagrange equations for a rigid body with quaternions