Brownian motion in an optical trap

Here we will study the motion of a particle in an optical trap, or for that matter, any classical over-damped particle in a quadratic potential at finite temperature. This should help you to understand a lot of important concepts such as power spectra, correlation functions, thermal equilibrium, and stochastic processes.

Using noise.py, you should be able to fill in most of hw5/harmonic_hints.py to produce a working version harmonic.py. It is important that you do this so you understand better how analysis of noise actually works.

Now run your file harmonic.py. You should see a particle attached to a spring, moving around quite randomly. The represents the solution to the Langevin equation when inertia is neglected.

$γ \dot{x} = - k x + n (t)$ $\gamma {\dot x} = -kx + n(t)$

Here $γ$ $\gamma$ is a drag coefficient of the particle in solution, $k$ $k$ is the spring constant of the trap, and $n (t)$ $n(t)$ is random noise. The random noise is called "white noise" that was discussed in the previous problem. We will talk more about this below.

This equation is not derived but assumed. We'll see that it works very nicely for reasons that'll become clear as we go on. The noise is suppose to be due to solvent applying forces to the particle due to thermal effects. Note that in the simulation $k = γ = 1$ $k=\gamma = 1$ , and the temperature has been set to $k_{B} T = 1$ $k_B T= 1$ as well. Because the system is linear, we only need consider its motion in one dimension at a time, which simplifies the analysis.

1. Discretize the Langevin eqn in time to show that it gives the form shown in the python code under the function "run_spring".

The noise has been taken from a Gaussian (standard normal) distribution, as you can see. The amplitude was set to depend on the temperature as shown as well. It also depends on the "dt", the time step for discretization. We'll try to determine why that is below.

2. For a linear spring at finite temperature, calculate $⟨ x^{2} ⟩$ $\langle x^2 \rangle$ using the equipartition theorem. This result is not derived from the Langevin equation, it is on much firmer foundation than that, namely equilibrium statistical mechanics, i.e. the Gibbs distribution.

3. Calculate the probability of finding the particle at position $x$ $x$ . Again, you can get this from the Gibbs distribution (probability $\propto \exp (- e n e r g y / (k_{B} T))$ $\propto \exp(-energy/(k_B T))$ .

4. Measure the same distribution as in 3 by running the simulation. You can change the code to set graphics=0 to turn off the graphics. Note that run_spring returns a list containing the position of the particle as a function of time. This is just what you need to get a histogram of positions. This is done in "compare_to_histogram" using the scipy function "hist". "hist" returns three things, the second of which, "bins" you need to use. It's an array containing the points on the x axis of the plot. You can also compare the simulation results with your answer to question 2, as the standard deviation is printed out and called "sigma". Note the code compares the simulation result with the theoretical result you should have calculated in question 3. Again, make sure you turn off the graphics flag or your simulation will take forever to complete! Does your value of sigma agree with what you computed in 2? Try changing the temperature to 0.25 to make sure it works in that case as well.

Now you should have complete agreement with the Gibbs distribution, that is equilibrium statistical properties, and the simulation of the Langevin eqn. We now turn to understanding the dynamics of this better.

First understand what happens to the Langevin equation if at t=0, we turn off the noise term and start the particle at position $x_{0}$ $x_0$ . You can now solve for the position of the particle as a function of time. Now consider what happens if we add noise. Since the equation is linear, we get the same solution except an extra random part added to it.

5. So now consider starting the particle at position $x = x_{0}$ $x = x_0$ and watching the trajectory. Consider doing this many times and averaging over all trajectories. We have a deterministic zero temperature part, that's always the same, and a random part. What happens when we average over many trajectories. From this calculate $⟨ x (t) ⟩$ $\langle x(t)\rangle$ with the condition that we always start at $t = 0$ $t=0$ and the same point $x = x_{0}$ $x = x_0$ .

Now we use "Onsager's Regression Hypothesis" (which one can really prove in this case). The time evolution of the system when we force it to start at $x_{0}$ $x_0$ is statistically identical to what it would have been if it got to $x = x_{0}$ $x = x_0$ by itself. So if we know $⟨ x (t) ⟩$ $\langle x(t)\rangle$ with $x (0)$ $x(0)$ given, we know $⟨ x (0) x (t) ⟩$ $\langle x(0) x(t)\rangle$ with $x (0)$ $x(0)$ given. What is it?

6. Now $x (0)$ $x(0)$ is not given but is taken from the Gibbs distribution. Use your answer from 2 to calculate $⟨ x (0) x (t) ⟩$ $\langle x(0) x(t)\rangle$ where here we are averaging over all possible $x (0)$ $x(0)$ as well. As discussed in problem 1, this is called the "autocorrelation function" sometimes just the "correlation function". We'll call it $C (t)$ $C(t)$ .

Note that $C (t) = ⟨ x (0) x (t) ⟩ = ⟨ x (- t) x (0) ⟩ = ⟨ x (0) x (- t) ⟩ = C (- t)$ $C(t) = \langle x(0) x(t)\rangle = \langle x(-t) x(0)\rangle = \langle x(0) x(-t)\rangle = C(-t)$ . So if we know $C$ $C$ for $t$ $t$ > 0 we also know it for $t < 0$ $t < 0$ . This is important in 8.

We have the positions of the particle as a function of time in the list "a". We can Fourier Transform this and square it. As discussed in problem 1, this is called the "power spectrum". We can now use the Wiener Khinchin Theorem Links to an external site. as discussed in problem 1.

7. As in problem 1, we can compute the correlation function, look at the function "compute_correlation". The only difference is that the noise has been generated a different way, but the way it's analyzed is the same. We apply the above theorem to compute it. Then this is compared with the analytical solution. Are the fitted values found in that function agreement with what you'd predict?

8. Compute the power spectrum. This is a little tricky because it will be quite noisy if not smoothed. To smooth it out over some window in frequency, see "compute_smoothed_power_spect". Here we're convolving the spectrum with a Gaussian filter using the convolution theorem. Use the Wiener-Khinchin Theorem to Fourier Transform your analytic result to obtain the power spectrum analytically. How does this analytic solution compare with the numerical results?

Important note about C(t). We've obtained C(t) for t>0. Make sure you figure it out correctly for t<0. Think about the physics of it. If it's blowing up in an unphysical way, it's not correct.

9. (Optional) Now we can answer the question we posed at the beginning: how to choose the amplitude of the noise term. To do this:

(a) Fourier Transform the Langevin equation. (Remember what happens to the Fourier Transform of derivatives). Then solve for

{\hat{x}}_{ω}

${\hat x}_\omega$ in terms of

{\hat{n}}_{ω}

${\hat n}_\omega$ ,

k

$k$ , and

γ

$\gamma$ .

Now we described the noise that we use in this equation, which is "white". That means that its correlation function $⟨ n (t) n (t^{'}) ⟩ = c δ (t - t^{'})$ $\langle n(t) n(t')\rangle = c \delta(t-t')$ . That is, it's correlation function is a delta function. At the moment $c$ $c$ is not determined and we're about to determine it. Note the by the Wiener-Khinchin theorem, $⟨ | {\hat{n}}_{ω} |^{2} ⟩ = c T_{m e a s}$ $\langle |{\hat n}_\omega|^2\rangle = c T_{meas}$ is frequency independent (which is why it's called white, in analogy with the spectrum of light).

(b) Now square (a) to get

⟨ | {\hat{x}}_{ω} |^{2} ⟩

$\langle |{\hat x}_\omega|^2\rangle$ . The right hand side should have an

⟨ | {\hat{n}}_{ω} |^{2} ⟩

$\langle |{\hat n}_\omega|^2\rangle$ which you can write as

c

$c$ . Compare this form with the form you derived in 8. This then gives a formula for

c

$c$ . Compare this with the amplitude of noise that is used in the code and see if they agree.

10. Get parameters from a real experiment on an optical trap to estimate the power spectrum that you'd expect to see in a real system. For example, you could choose $1 μ m$ $1 \mu m$ diameter polystyrene bead. and a spring constant of 0.5 pN/nm. What is the variance in the bead position that you would expect? If you didn't do 9, you can get the power spectrum here Links to an external site..

11. In real experiments, the fluctuations in position are often large compared to the positions that are useful to measure. Give a method for how the noise in the answer can be reduced to a level where interesting behavior can be observed. Give an example of one experiment that employs such a technique.

12. Give an example of an experiment probing conformations of a single molecule, where you can think of the particle as being in a potential that no longer has a single minimum and the noise plays an important role in understanding the molecule's internal states.