Equation of motion - SDE to be solved

$\ddot{q}(t) + \Gamma_0\dot{q}(t) + \Omega_0^2 q(t) - \dfrac{1}{m} F(t) = 0 $

where q = x, y or z

Where $F(t) = \mathcal{F}_{fluct}(t) + F_{feedback}(t)$

Taken from page 46 of 'Dynamics of optically levitated nanoparticles in high vacuum' - Thesis by Jan Gieseler

Using $\mathcal{F}_{fluct}(t) = \sqrt{2m \Gamma_0 k_B T_0}\dfrac{dW(t)}{dt}$

and $F_{feedback}(t) = \Omega_0 \eta q^2 \dot{q}$

Taken from page 49 of 'Dynamics of optically levitated nanoparticles in high vacuum' - Thesis by Jan Gieseler

we get the following SDE:

$\dfrac{d^2q(t)}{dt^2} + (\Gamma_0 - \Omega_0 \eta q(t)^2)\dfrac{dq(t)}{dt} + \Omega_0^2 q(t) - \sqrt{\dfrac{2\Gamma_0 k_B T_0}{m}} \dfrac{dW(t)}{dt} = 0$

split into 2 first order ODE/SDE s

letting $v = \dfrac{dq}{dt}$

$\dfrac{dv(t)}{dt} + (\Gamma_0 - \Omega_0 \eta q(t)^2)v + \Omega_0^2 q(t) - \sqrt{\dfrac{2\Gamma_0 k_B T_0}{m}} \dfrac{dW(t)}{dt} = 0$

therefore

$\dfrac{dv(t)}{dt} = -(\Gamma_0 - \Omega_0 \eta q(t)^2)v - \Omega_0^2 q(t) + \sqrt{\dfrac{2\Gamma_0 k_B T_0}{m}} \dfrac{dW(t)}{dt} $

$v = \dfrac{dq}{dt}$ therefore $dq = v~dt$

\begin{align} dq&=v\,dt\\ dv&=[-(\Gamma_0-\Omega_0 \eta q(t)^2)v(t) - \Omega_0^2 q(t)]\,dt + \sqrt{\frac{2\Gamma_0 k_B T_0}m}\,dW \end{align}

Apply Milstein Method to solve

Consider the autonomous Itō stochastic differential equation

${\mathrm {d}}X_{t}=a(X_{t})\,{\mathrm {d}}t+b(X_{t})\,{\mathrm {d}}W_{t}$

Taking $X_t = q_t$ for the 1st equation above (i.e. $dq = v~dt$) we get:

$$ a(q_t) = v $$$$ b(q_t) = 0 $$

Taking $X_t = v_t$ for the 2nd equation above (i.e. $dv = ...$) we get:

$$a(v_t) = -(\Gamma_0-\Omega_0\eta q(t)^2)v - \Omega_0^2 q(t)$$

$$b(v_t) = \sqrt{\dfrac{2\Gamma_0 k_B T_0}m}$$

${\displaystyle b'(v_{t})=0}$ therefore the diffusion term does not depend on ${\displaystyle v_{t}}$ , the Milstein's method in this case is therefore equivalent to the Euler–Maruyama method.

We then construct these functions in python:

Using values obtained from fitting to data from a real particle we set the following constant values describing the system. Cooling has been assumed to be off by setting $\eta = 0$.

partition the interval [0, T] into N equal subintervals of width $\Delta t>0$: $ 0=\tau _{0}<\tau _{1}<\dots <\tau _{N}=T{\text{ with }}\tau _{n}:=n\Delta t{\text{ and }}\Delta t={\frac {T}{N}}$

set $Y_{0}=x_{0}$

Generate independent and identically distributed normal random variables with expected value 0 and variance dt

Apply Milstein's method (Euler Maruyama if $b'(Y_{n}) = 0$ as is the case here):

recursively define $Y_{n}$ for $ 1\leq n\leq N $ by

$ Y_{{n+1}}=Y_{n}+a(Y_{n})\Delta t+b(Y_{n})\Delta W_{n}+{\frac {1}{2}}b(Y_{n})b'(Y_{n})\left((\Delta W_{n})^{2}-\Delta t\right)$

Perform this for the 2 first order differential equations:

We now have an array of positions, $v$, and velocities $p$ with time $t$.

Alternatively we can use a derivative-free version of Milsteins method as a two-stage kind-of Runge-Kutta method, documented in wikipedia (https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_method_%28SDE%29) or the original in arxiv.org https://arxiv.org/pdf/1210.0933.pdf.

The form of $F_{feedback}(t)$ is still questionable

On page 49 of 'Dynamics of optically levitated nanoparticles in high vacuum' - Thesis by Jan Gieseler he uses the form: $F_{feedback}(t) = \Omega_0 \eta q^2 \dot{q}$

On page 2 of 'Parametric feeedback cooling of levitated optomechancs in a parabolic mirror trap' Paper by Jamie and Muddassar they use the form: $F_{feedback}(t) = \dfrac{\Omega_0 \eta q^2 \dot{q}}{q_0^2}$ where $q_0$ is the amplitude of the motion: $q(t) = q_0(sin(\omega_0t)$

However it always shows up as a term $\delta \Gamma$ like so:

$\dfrac{d^2q(t)}{dt^2} + (\Gamma_0 + \delta \Gamma)\dfrac{dq(t)}{dt} + \Omega_0^2 q(t) - \sqrt{\dfrac{2\Gamma_0 k_B T_0}{m}} \dfrac{dW(t)}{dt} = 0$

By fitting to data we extract the following 3 parameters:

1) $A = \gamma^2 \dfrac{k_B T_0}{\pi m}\Gamma_0 $

Where:

• $\gamma$ is the conversion factor between Volts and nanometres. This parameterises the amount of light/ number of photons collected from the nanoparticle. With unchanged allignment and the same particle this should remain constant with changes in pressure.
• $m$ is the mass of the particle, a constant
• $T_0$ is the temperature of the environment
• $\Gamma_0$ the damping due to the environment only

2) $\Omega_0$ - the natural frequency at this trapping power

3) $\Gamma$ - the total damping on the system including environment and feedback etc...

By taking a reference save with no cooling we have $\Gamma = \Gamma_0$ and therefore we can extract $A' = \gamma^2 \dfrac{k_B T_0}{\pi m}$. Since $A'$ should be constant with pressure we can therefore extract $\Gamma_0$ at any pressure (if we have a reference save and therefore a value of $A'$) and therefore can extract $\delta \Gamma$, the damping due to cooling, we can then plug this into our SDE instead in order to include cooling in the SDE model.

For any dataset at any pressure we can do:

$\Gamma_0 = \dfrac{A}{A'}$

And then $\delta \Gamma = \Gamma - \Gamma_0$

Using this form and the same derivation as above we arrive at the following form of the 2 1st order differential equations:

\begin{align} dq&=v\,dt\\ dv&=[-(\Gamma_0 + \delta \Gamma)v(t) - \Omega_0^2 v(t)]\,dt + \sqrt{\frac{2\Gamma_0 k_B T_0}m}\,dW \end{align}

values below are taken from a ~1e-2 mbar cooled save