# Harmonic Oscillator subjected to Random Driving Force

I am trying to solve the below problem, but I am unsure about my attempted solution.

## Problem statement

This problem comes from exercise 2 of the notes on Green's function.

## My attempted solution

The equation of motion for harmonic oscillator is

$$\left[\frac{\text{d}^2}{\text{d}t^2}+2\gamma\frac{\text{d}}{\text{d}t}+\omega_0^2\right]x(t)=\frac{f(t)}{m} \tag{1}$$

For over-damped condition $$\gamma > \omega_0$$, the Green's function is given as $$G(t,t')=\Theta(t-t')\frac{1}{\Omega}e^{-\gamma(t-t')}\sinh\left[\Omega(t-t')\right] \tag{2}$$ where $$\Omega=\sqrt{\gamma^2-\omega_0^2}$$ and $$\Theta(t-t')$$ is step function: $$\Theta(t-t')= \begin{cases} 0 &\mbox{if } t

The solution of displacement $$x(t)$$ at time $$t$$ is given by Green's function: $$x(t)=\text{Initial Conditions}+\int_0^t G(t,t')\frac{f(t')}{m}\text{d}t' \tag{4}$$ For large time limit $$t\to\infty$$, the initial conditions are damped out. Therefore, the ensemble average of displacement can be written as: $$\left=\int_0^t G(t,t')\frac{\left}{m}\text{d}t' \tag{5}$$

The mean displacement correlation function $$\left$$ is $$\begin{split} \left =&\int_0^{t_1}\int_0^{t_2} G(t_1,t')G(t_2,t'') \frac{\left}{m^2}\text{d}t'\text{d}t''\\ =&\frac{A}{m^2}\int_0^{t_1} G(t_1,t')G(t_2,t') \text{d}t' \\ =&\frac{A}{4\Omega^2 m^2}\Biggl\{ \frac{1}{2(\gamma-\Omega)}\left[e^{(\gamma-\Omega)|t_1-t_2|}-e^{-(\gamma-\Omega)(t_1+t_2)}\right] \\ &+\frac{1}{2(\gamma+\Omega)}\left[e^{(\gamma+\Omega)|t_1-t_2|}-e^{-(\gamma+\Omega)(t_1+t_2)}\right] \\ &-\frac{1}{2\gamma}\left[e^{\Omega|t_1-t_2|}+e^{-\Omega|t_1-t_2|}\right]\left[e^{\gamma|t_1-t_2|}-e^{-\gamma(t_1+t_2)}\right]\Biggr\} \end{split} \tag{6}$$

The ensemble average of fluctuation should be dependent only on $$\Delta t$$ but not $$t$$, therefore $$\begin{split} \left<\left[x(t+\Delta t)-x(t)\right]\right> =&\left<\left[x(\Delta t)-x(0)\right]\right> =&\left =&\int_0^{\Delta t} G(\Delta t,t')\frac{\left}{m}\text{d}t' \end{split} \tag{7}$$ The ensemble average of squared fluctuation is $$\begin{split} \left<\left[x(t+\Delta t)-x(t)\right]^2\right> =&\left<\left[x(\Delta t)\right]^2\right> \\ =&\int_0^{\Delta t}\int_0^{\Delta t} G(\Delta t,t')G(\Delta t,t'')\frac{\left}{m^2}\text{d}t'\text{d}t'' \\ =&\frac{A}{m^2}\int_0^{\Delta t} \left[G(\Delta t,t')\right]^2 \text{d}t' \\ =&\frac{A}{4\Omega^2 m^2}\left[ \frac{1-e^{-2(\gamma-\Omega)\Delta t}}{2(\gamma-\Omega)}- \frac{1-e^{-2\gamma\Delta t}}{\gamma}+ \frac{1-e^{-2(\gamma+\Omega)\Delta t}}{2(\gamma+\Omega)} \right] \end{split} \tag{8}$$

## My question

It seems that my solution to the average fluctuation $$\left<\left[x(t+\Delta t)-x(t)\right]^2\right>$$ is not proportional to $$\Delta t$$, and I have no idea why. The original lecture notes did not cover anything on Brownian motion or statistical physics, so I just google the notion of $$\left<\cdot\right>$$, and probably misuse it somehow. Please kindly help to comment on any of my conceptual errors and miscalculations. Many thanks!

## Edit 1

Thanks to the hints given by @cyd, the steps starting from equation (6) are corrected as below. Obviously the displacement correlation function should not depend on $$(t_1+t_2)$$.

The displacement correlation function $$\left$$ is $$\begin{split} \left =&\int_{-\infty}^{t_1}\int_{-\infty}^{t_2} G(t_1,t')G(t_2,t'')\frac{\left}{m^2} \text{d}t'\text{d}t'' \\ =&\frac{A}{m^2}\int_{-\infty}^{t_1}\int_{-\infty}^{t_2} G(t_1,t')G(t_2,t'')\delta(t'-t'') \text{d}t'\text{d}t'' \\ =&\frac{A}{m^2}\int_{-\infty}^{t_1} G(t_1,t') \left[\int_{-\infty}^{t_2} G(t_2,t'')\delta(t'-t'') \text{d}t''\right]\text{d}t' \\ =&\frac{A}{m^2}\int_{-\infty}^{\min(t_1,t_2)} G(t_1,t') G(t_2,t')\text{d}t' \\ =&\frac{A}{m^2}\frac{1}{8\gamma\Omega}\left[\frac{e^{-(\gamma-\Omega)|t_1-t_2|}}{\gamma-\Omega}-\frac{e^{-(\gamma+\Omega)|t_1-t_2|}}{\gamma+\Omega}\right] \end{split} \tag{6A}$$ The mean squared displacements are calculated using equation (6A): $$\left<[x(t)]^2\right>=\left<[x(t+\Delta t)]^2\right> =\frac{A}{m^2}\frac{1}{8\gamma\Omega}\left[\frac{1}{\gamma-\Omega}-\frac{1}{\gamma+\Omega}\right] \tag{7.1A}$$ $$\left =\frac{A}{m^2}\frac{1}{8\gamma\Omega}\left[\frac{e^{-(\gamma-\Omega)\Delta t}}{\gamma-\Omega}-\frac{e^{-(\gamma+\Omega)\Delta t}}{\gamma+\Omega}\right] \tag{7.2A}$$ The diffusion is $$\begin{split} \left<[x(t+\Delta t)-x(t)]^2\right> &=\left<[x(t+\Delta t)]^2\right>+\left<[x(t)]^2\right>-2\left \\ &=\frac{A}{m^2}\frac{1}{4\gamma\Omega}\left[\frac{1-e^{-(\gamma-\Omega)\Delta t}}{\gamma-\Omega}-\frac{1-e^{-(\gamma+\Omega)\Delta t}}{\gamma+\Omega}\right] \end{split} \tag{8A}$$ The diffusion in short and long time limits can be obtained by using Taylor expansion: $$\begin{equation*} \left<[x(t+\Delta t)-x(t)]^2\right>=\frac{A}{4\gamma m^2}\times \begin{cases} \Delta t^2 &\text{as } \Delta t\to0 \\ 2/(\gamma^2-\Omega^2) &\text{as } \Delta t\to\infty \end{cases} \tag{9A} \end{equation*}$$

• Taylor expand the exponentials. Keep the expansion to first order in $\Delta t$, since it's small. – d_b Nov 9 '19 at 3:35
• @d_b: have tried Taylor expanding the result in equation (8) using $e^x\approx 1+x$, but that seems to cancel out all the 1st-order terms in $\Delta t$. – York Tsang Nov 9 '19 at 3:49
• Ah I see. Sorry I didn't check that before commenting. – d_b Nov 9 '19 at 4:40
• I would argue the solution is neither Brownian or harmonic. Multiplying $2$ deterministic linear $\sinh()$ solutions together doesn't create Brownian motion. The given $f(t)$ delta function only ensures the $2$ particles collide once. The harmonic eigenfunctions are also solutions of the DE - it just depends on the ratio of $\frac {\omega_0} {\gamma}$. And $<x(t_{1})x(t_{2})>$ isn't a correlation function without a statistical model, i.e., a probability measure - it's just multiplication. – Cinaed Simson Nov 17 '19 at 22:48
• @CinaedSimson: Title and tags are edited to avoid confusion. Thanks for comment! – York Tsang Nov 18 '19 at 5:24

In the second line of Eq. (6), you have to be a bit careful about when the Green's function is zero. The expression should be

$$\frac{A}{m^2}\int^{\textrm{min}(t_1,t_2)} dt'\; G(t_1, t') G(t_2, t')$$

where $$\textrm{min}(t_1,t_2)$$ is the smaller of $$t_1$$ and $$t_2$$.

Eq. (7) doesn't look, right; in this kind of calculations you always want to take the mean of the product of two quantities. Something like this:

$$\langle [x(t+\Delta t) - x(t)]^2\rangle = \langle x(t+\Delta t)^2\rangle + \langle x(t)^2\rangle - 2 \langle x(t+\Delta t) x(t)\rangle$$

As an intermediate step, you should try to calculate the mean square displacement, $$\langle x(t)^2\rangle$$.

# Reply to Edit 1 (updated)

Hmm, you're right, there's something funky going on here which I didn't expect.

I can verify that your newly-derived expression for $$\langle [x(t+\Delta t) - x(t)]^2\rangle$$ is correct, and it indeed doesn't scale linearly for $$\Delta t \rightarrow 0$$ or $$\Delta t \rightarrow \infty$$. But, the interesting thing is that for intermediate times, it's linear:

What's happening here is that the mean squared deviation is a sum of two terms containing exponential decays:

$$\langle [x(t+\Delta t) - x(t)]^2\rangle = \frac{A}{4m^2\gamma\Omega} \left[\frac{1-e^{-\Gamma_1 \Delta t}}{\Gamma_1} + \frac{1-e^{-\Gamma_2 \Delta t}}{\Gamma_2} \right]$$

where

$$\Gamma_1 = \gamma - \Omega,\;\; \Gamma_2 = \gamma + \Omega$$

The first exponential (with decay factor $$\Gamma_1$$) decays slower than the second exponential (which has decay factor $$\Gamma_2 > \Gamma_1$$). So when the second exponential has died away, we end up with

\begin{aligned}\langle [x(t+\Delta t) - x(t)]^2\rangle &\approx \frac{A}{4m^2\gamma\Omega} \left[\frac{1-e^{-\Gamma_1 \Delta t}}{\Gamma_1} + \frac{1}{\Gamma_2} \right] \quad (\textrm{for}\;\, 1/\Gamma_2 \lesssim t \lesssim 1/\Gamma_1) \\ &\approx \frac{A\,\Delta t}{4m^2\gamma\Omega} + \textrm{constant offset} \end{aligned}

For context, this is a thinly-disguised version of a famous problem, the Langevin oscillator. That problem is normally solved not using the Green's function method, but by playing around with different correlation functions and using the equipartition theorem. I'm not sure why there's a discrepancy between the two approaches (e.g., the existence of the offset term); maybe there's some subtle order-of-limits issue at play here. Will have to get back to you on that.

• Thank you for your reply, professor! – York Tsang Nov 15 '19 at 8:55
• The original free particle solution given by Einstein and Langevin is indeed $\left<(\Delta x)^2\right>\propto |\Delta t|$ for all $\Delta t$, but its derivation ignores the effect of particle inertia, or equivalently assumes $\Delta t \gg \frac{1}{2\gamma}$. Later solution by Ornstein–Fürth formula or Ornstein–Uhlenbeck theory is $\left<(\Delta x)^2\right>\propto |\Delta t|^2 \text{ as } |\Delta t|\to 0$, of which the same result can be obtained by using Green's function for free particle. For the harmonic oscillator, its solution reduces to free particle at the limit of $\omega_0\to 0$ – York Tsang Nov 18 '19 at 6:28
• Very nice. I guess I'll have to update that problem. Thanks for actually working through it :-P – cyd Nov 18 '19 at 9:23

I think your definition of the average displacement is missing a normalisation factor of 1/t

• @lux: thank you. Can you point out which equation I need to add a normalization factor $1/t$? like equation (5)? but the unit will not be matched then. – York Tsang Nov 10 '19 at 1:08
• I certainly don't follow why this is a good definition of the average. I thn k 5 needs a 1/t – lux Nov 10 '19 at 1:28