Sensing with the harmonic oscillator

TL;DR

This paper reviews how harmonic oscillators can be used as sensitive force or acceleration sensors, deriving analytical sensitivity expressions considering various noise sources and frequency jitters, and discusses potential sensitivity improvements and applications.

Contribution

It provides a comprehensive analytical framework for evaluating and enhancing the sensitivity of harmonic oscillator sensors under different noise conditions and frequency fluctuations.

Findings

Derived analytical sensitivity expressions for harmonic oscillator sensors.

Identified conditions for sensitivity improvements beyond standard limits.

Discussed practical applications and examples of oscillator-based sensing.

Abstract

A system obeying the harmonic oscillator equation of motion can be used as a force or proper acceleration sensor. In this short review we derive analytical expressions for the sensitivity of such sensors in a range of different situations, considering noise of thermal and measurement origins and a formalism for dealing with oscillators whose natural frequency $\omega_0$ jitters. A special case where the sensitivity can be improved beyond the standard expressions and some applications with examples are also discussed.

Equations

$$m\ddot{x}+m\Gamma \dot{x}+m{\omega }_0^{2}x=g(t)$$

$$x(t)=h(t)\ast g(t),$$

$$H(\omega )=\frac{1/m}{{\omega }_0^2-{\omega }^2+i\omega {\Gamma }^2}.$$

$$S_{xx}(\omega )=|H(\omega ){|}^2\cdot S_{gg}(\omega ).$$

$$S_{gg}(\omega )=|H(\omega ){|}^{-2}\cdot S_{xx}(\omega ).$$

$$S_{xx}=|H(\omega ){|}^2{S}_{gg}+|H(\omega ){|}^2{S}_{ww}+S_{uu}=|H(\omega ){|}^2{S}_{gg}+|H(\omega ){|}^2{\sigma }_1^2+{\sigma }_2^2,$$

$$|H(\omega ){|}^{-2}\cdot S_{xx}$$

$$=S_{gg}+{\sigma }_1^2+|H(\omega ){|}^{-2}{\sigma }_2^2.$$

$$\text{SNR}=\sqrt{\frac{S_{gg}}{{\sigma }_1^2+|H(\omega ){|}^{-2}{\sigma }_2^2}}.$$

$$\sqrt{S_{gg}}$$

$$=\sqrt{2m\Gamma k_{B}T+{\sigma }_2^2\frac{m^2(4Q^2-1){\omega }_0^4}{4Q^4}}\simeq \sqrt{2m\Gamma k_{B}T+{\sigma }_2^2{m}^2{\omega }_0^2{\Gamma }^2}$$

$$\text{ASD}=\sqrt{2m\Gamma k_{B}T}$$

$$\text{ASD}=\sqrt{{\sigma }_1^2+|H(\omega ){|}^{-2}{\sigma }_2^2}.(\text{N}/\sqrt{\text{Hz}})$$

$$g_{\text{min}}=\sqrt{2\Delta f{S}_{gg}}>\sqrt{4m\Gamma \Delta f{k}_{B}T}.$$

$$\mathbb{E}[\text{ASD}(\omega )]$$

$$={\int }_{-\infty }^{\infty }\sqrt{{\sigma }_1^2+|H(\omega ){|}^{-2}{\sigma }_2^2}\cdot f_{\Omega }(u)\mathrm{d}u.$$

$$\mathbb{E}[\text{ASD}(\omega )]$$

$$\simeq {\int }_{-\infty }^{\infty }\sqrt{{\sigma }_1^2}\cdot f_{\Omega }(u)\mathrm{d}u.={\sigma }_1,$$

$$\sqrt{S_{gg}(\omega )}=\sqrt{\mathbb{E}(|H(\omega ){|}^{-2})\cdot S_{xx}(\omega )},$$

$$m\ddot{x}+m\Gamma \dot{x}+m{\omega }_0^{2}x=-m{k}_d\dot{x},$$

$$H(\omega )=\frac{1/m}{{\omega }_0^2-{\omega }^2+i\omega (\Gamma +k_d{)}^2},$$

$$m(t)=(g(t)+w(t))\ast h(t)+u(t)$$

$$\mathcal{F}\left((g(t)+w(t))\ast h(t)+u(t)\right)=\mathcal{F}\left(g(t)\ast h(t)\right)+\mathcal{F}\left(w(t)\ast h(t)\right)+\mathcal{F}\left(u(t)\right)$$

$$=\pi (\delta (\omega -{\omega }_d)+\delta (\omega +{\omega }_d))\cdot H(\omega )+$$

$$\mathcal{F}\left(\left(\begin{array}{c}\cos \sqrt{a}t\sin \sqrt{a}t/\sqrt{a}\end{array}\right)\cdot {\int }_0^{t}b{e}^{\frac{\Gamma }{2}(r-t)}\cdot \left(\begin{array}{c}-\sin \sqrt{a}r/\sqrt{a}\\ \cos \sqrt{a}r\end{array}\right)\mathrm{d}{W}_r\right)+\mathcal{F}\left(u(t)\right),$$

$$\mathbb{E}[\mathcal{F}\left((g(t)+w(t))\ast h(t)+u(t)\right)]=\mathbb{E}[\pi (\delta (\omega -{\omega }_d)+\delta (\omega +{\omega }_d))\cdot H(\omega )+$$

$$\mathcal{F}\left(\left(\begin{array}{c}\cos \sqrt{a}t\sin \sqrt{a}t/\sqrt{a}\end{array}\right)\cdot {\int }_0^{t}b{e}^{\frac{\Gamma }{2}(r-t)}\cdot \left(\begin{array}{c}-\sin \sqrt{a}r/\sqrt{a}\\ \cos \sqrt{a}r\end{array}\right)\mathrm{d}{W}_r\right)+\mathcal{F}\left(u(t)\right)]$$

$$=\pi (\delta (\omega -{\omega }_d)+\delta (\omega +{\omega }_d))\cdot H(\omega )+$$

$$\mathcal{F}\left(\left(\begin{array}{c}\cos \sqrt{a}t\sin \sqrt{a}t/\sqrt{a}\end{array}\right)\cdot \mathbb{E}\left[{\int }_0^{t}b{e}^{\frac{\Gamma }{2}(r-t)}\cdot \left(\begin{array}{c}-\sin \sqrt{a}r/\sqrt{a}\\ \cos \sqrt{a}r\end{array}\right)\mathrm{d}{W}_r\right)\right]+\mathcal{F}\left(\mathbb{E}[u(t)]\right)$$

$$=\pi (\delta (\omega -{\omega }_d)+\delta (\omega +{\omega }_d))\cdot H(\omega ),$$

$$m\ddot{x}+m\Gamma \dot{x}+m{\omega }_0^{2}x=0\to m\ddot{x}+m\Gamma (\dot{x}-\dot{y})+m{\omega }_0^2(x-y)=0$$

$$m\ddot{z}+m\Gamma \dot{z}+m{\omega }_0^{2}z=-m\ddot{y}.$$

$$H(\omega )=\frac{-{\omega }^2}{{\omega }_0^2-{\omega }^2+i\omega {\Gamma }^2}.$$

$$Z(\omega )=H(\omega )Y(\omega ){=}_{\omega \gg {\omega }_0}Y(\omega )$$

$$z(t)\simeq y(t).$$

Full text

\SectionNumbersOn

1 Introduction

Linear force and acceleration sensors are often based on systems whose equation of motion is well approximated by a harmonic oscillator. The inertial sensing market is nowadays dominated by micro-electromechanical systems (MEMS)1, devices with sizes ranging between 20 $\mu$m to a mm that can be found in most modern cellphones111In fact, smartphones contain Inertial Measurement Units (IMU), devices incorporating both MEMS Accelerometers and MEMS gyroscopes, which provide linear acceleration and torque measurements respectively., but any harmonic oscillator sensor follows the same working principles. For instance, a massive pendulum around its stable equilibrium position can be used if a large oscillator mass is convenient; if, on the contrary, the mass is required to be small, a levitated nano-particle or a micro-cantilever can be chosen.

Conceptually, force sensors and linear inertial sensors (also known as accelerometers) detect the effect of a driving force as a displacement of the oscillator’s mass; this displacement can be measured by various methods222For instance, modern MEMS accelerometers often consist of little more than a cantilever beam with a proof mass (also known as seismic mass), whose movement is measured as a time dependent capacitance. Other accelerometers work by detecting the electric field due to the strain applied to a piezoelectric crystal., depending on the sensor. In any case, the displacement is created when the driving modifies the equation of motion: the subtle difference is that in force sensors the driving only acts on the oscillator’s mass, whereas in an accelerometer the driving acts on the oscillator’s housing.

In this short review we derive analytical expressions for the sensitivity of such sensors in a range of different situations, considering noise of thermal and measurement origins and a formalism for dealing with oscillators whose natural frequency $\omega_{0}$ jitters. A special case where the sensitivity can be improved beyond the standard expressions and some applications with examples are also discussed.

2 Force sensing

2.1 Principles of force sensing

The harmonic oscillator is a 2nd order constant coefficient linear ordinary differential equation. In the most general case, the equation reads

[TABLE]

where $m\Gamma\dot{x}$ is a damping force, $m\omega_{0}^{2}x$ is a restoring force and $g(t)$ is an (external) driving force. Force sensing is the detection of these forces $g(t)$ acting upon the oscillator’s mass. Since the harmonic oscillator is a linear time-invariant system (see supplementary), the response of the system to a deterministic (arbitrary) driving function $g(t)$ will be

[TABLE]

where $h(t)$ is the impulse response of the system. By the convolution theorem, $X(\omega)=H(\omega)\cdot G(\omega)$, where the capital letters indicate Fourier transforms. Recall the transfer function of the harmonic oscillator (see supplementary) is

[TABLE]

If the driving function is stochastic, but wide sense stationary333This assumption is important, since the following expressions assume no transient behaviours. When $\Gamma\ll 1$, past transients can contaminate the signal, leading to worse signal to noise ratio (SNR) than expected., this equality still holds in the power spectral density (PSD) sense

[TABLE]

Therefore, given a certain experimentally measurable $x(t)$, to obtain its originating $S_{gg}$ we would need to calculate

[TABLE]

In a more realistic scenario, the oscillator will be driven by thermal noise444As can be seen by the fluctuation-dissipation theorem2. The engineer/scientist should decide whether this Brownian noise can be neglected or should be taken into account, depending on the harmonic oscillator’s mass., a stochastic driving that can be modelled as a white noise $w(t)$ with zero mean and autocorrelation function $R_{w}(\tau)=\sigma_{1}^{2}\delta(\tau)$. There will also be a random noise $u(t)$ originating from our measuring device, that we can consider additive Gaussian white noise (AGWN) (with $R_{u}(\tau)=\sigma_{2}^{2}\delta(\tau)$), summed after the harmonic oscillator frequency response. Therefore, the measured power spectral density of a driving force $g(t)$ is

[TABLE]

and our estimate of $S_{gg}$ will be

[TABLE]

In this expression $S_{gg}$ is deterministic, $\sigma_{1}^{2}$ comes from the thermal noise and is independent of $\omega$ and $|H(\omega)|^{-2}\sigma_{2}^{2}$ has explicit $\omega$ dependency. The signal to noise ratio (SNR), defined as the square root555We use the square root to work with force units, instead of power units. This is, of course, arbitrary. of the ratio of the signal and noise powers, will thus take the expression

[TABLE]

This expression is bounded by $\sqrt{\frac{S_{gg}}{\sigma_{1}^{2}}}$: this is a hard limit that cannot be improved with this system and assumptions. However, the $|H(\omega)|^{-2}\sigma_{2}^{2}$ term can be minimised by using $\omega$ at resonance.

In any case, to get the minimum measurable force $g(t)$ we have to set a limit to what we can detect: usually this limit is $\text{SNR}>1$666Again, this limit is arbitrary: in principle if the system is ergodic and the driving force periodic, signals of any SNR can be detected. From this inequality we can obtain a bound for $g(t)$, which is where the minimum sensitivity expressions come from. Using the fluctuation-dissipation relationship $\sigma_{1}=\sqrt{2m\Gamma k_{B}T}$ and $|H(\omega)|^{2}$ at its maximum (i.e., at resonance) we get

[TABLE]

where after the inequality we assumed resonance, and the last equality is a good approximation for large enough $Q$. Whether we are measuring at resonance or not, if the measuring noise is very small compared to the thermal noise, the second term in the square root may be negligible, in which case we obtain the following expression for the amplitude spectral density (ASD)

[TABLE]

The “minimum resolvable” 3 force or force sensitivity is hence defined as the root mean square of the noise power spectral density (i.e., the ASD is just the square root of the PSD777For reference, check the LISA Pathfinder paper.). Notice that smaller sensitivity values are better, since the benchmark is the smallest signal that can be detected.

Expression (6) is useful to get an intuition of how the sensitivity scales with the oscillator parameters. For example, a smaller mass results in a better sensitivity888Assuming, of course, that noise has no $m$ dependency. For instance, with optically levitated nanoparticles, decreasing the mass increases the measurement noise, since smaller particles scatter less light. Thus, to obtain a signal that is comparable to that of bigger particles, we will need to amplify the measured signal, hence also amplifying measurement noise.. Bear in mind, however, that (6) is just an approximation and the full sensitivity has $\omega$ dependency,

[TABLE]

If we further integrate the power only on an interval of length $\Delta f$ in the Fourier domain999As is done in ref. 4., we obtain

[TABLE]

This expression of $g_{\text{min}}$ is also sometimes known as the force sensitivity of the system. Nevertheless, the previous ASD definition can (and is) also used under the same name. In general, the use of the force sensitivity defined as the ASD (and not as this $g_{\text{min}}$) should always be preferred, since its value doesn’t depend on how long the measurement is, it doesn’t make any assumptions on the bandwidth of the signal and any information of different sensitivities at different $\omega$ values is not lost.

It is worth pointing out that if one assumes that measurement noise can be neglected, the sensitivity is the same for all frequencies: it doesn’t matter if one measures at resonance or out of it101010In fact, measuring at resonance may be a bad idea if the driving signal has a big bandwidth, since the phase response changes a lot around the resonance peak.. This happens because the oscillator gain affects thermal noise and driving signal equally and the SNR stays constant. Therefore, for a fixed $m$, the obvious knobs for reducing the ASD are decreasing the temperature $T$ and/or $\Gamma$.

2.2 The case of a harmonic oscillator with jitter

If the response function of the harmonic oscillator jitters (i.e., is not constant with time), we can model the sensitivity in probabilistic terms. Assume the jittering is due to a stochastic natural frequency of the oscillator $\Omega$111111In contrast with the previous $\omega_{0}$., that now fluctuates with time, but has a certain stationary, well-known probability density function $f_{\Omega}(u)$. Then, the value of the ASD from (7) will also be stochastic, but we can still calculate its statistical moments. In particular, the expected value of the sensitivity will be

[TABLE]

Clearly this value depends on $\omega$, since the expected sensitivity will not be the same for every frequency. However, if, as before, measurement noise can be neglected, then

[TABLE]

recovering the expression of (6). Intuitively, this means that if the natural frequency $\Omega$ jitters only in a region where thermal noise is still dominant (for instance the standard deviation of the natural frequency, $\sigma_{\Omega}$, is less than $\Gamma/2$), the sensitivity is not affected, even if the force amplitude fluctuates due to the jitter. However, the actual value of the measured force will be affected if one naively uses equation (3) without taking into account that $H(\omega)$ is stochastic. To calculate the real expected force, we will need to, again, find the expected value of $H(\omega)$

[TABLE]

which can be thought of as an effective response function obtained by averaging the instantaneous responses at different times.

Real oscillators, and especially small ones (which are more susceptible to perturbations), will always have some jitter. The importance of the effect depends on the parameters of the oscillator: for instance, the effect can be neglected when the damping is large and the oscillator’s width is much broader than the jitter (i.e., $\Gamma\gg\sigma_{\Omega}$). However, for strongly underdamped systems – which have a very narrow response peak –, jitter needs to be taken into account. Beyond the probabilistic treatment that we just described, it is advisable to have some way of artificially increasing the damping in the harmonic oscillator. As long as this extra damping does not induce additional noise (i.e., it is a cold damping5), this has several advantages:

• •

Since the peak is broadened, the relative effect of jitter is reduced.

• •

As long as thermal noise is still considerably larger than measurement noise, the SNR (i.e., the sensitivity) will not be affected.

• •

It will reduce the effect of non-linearities in the oscillator. Low dissipation systems with narrow resonances are prone to large oscillation amplitudes: in this case, non-linearities can dominate121212This is a classical example..

• •

It will also reduce the memory of the system. Recall that, the lower the damping, the longer the memory of the impulse response. In practical terms, this means that if no additional damping is included we will keep measuring signal remnants from some past event long after the driving has stopped.

One possibility to implement this damping is with a feedback force $-mk_{d}\dot{x}$. The equation of motion becomes

[TABLE]

resulting in a new transfer function

[TABLE]

which has a bigger effective damping than before, while the thermal noise will still have spectral density $\sqrt{2k_{B}Tm\Gamma}$. Real feedback systems will never be completely noise-free, but for small gains the effect of this noise is negligible (see supplementary material of ref. 5).

2.3 Sensing a perfectly sinusoidal force

Assume now the force $g(t)$ we want to measure is a sinusoid, with a constant phase relation to a controlled reference signal (i.e., it is a perfect sinusoid). In this case, the sensitivity to this signal can be greatly increased if instead of working with PSDs one averages the Fourier transform measurements (note that this gives $(\text{N}/\sqrt{\text{Hz}})$ units directly). The argument is described in what follows: the measured signal $m(t)$ will take the expression

[TABLE]

where, as before, $g(t)$ is the driving force, $w(t)$ is a white thermal noise and $u(t)$ is an AWGN from the measuring device. By assumption, $g(t)=\cos w_{d}t$. Then, if we take the Fourier transform of the previous expression we get

[TABLE]

where we have used the properties of the Fourier transform, in the second summand we have substituted by the solution of a thermally driven harmonic oscillator (see the supplementary material for the derivation), $a=\omega_{0}^{2}-\Gamma^{2}/4$ and $W_{r}$ is a Wiener process parametrized by the time $r$. Finally, taking expected values and applying Fubini’s theorem

[TABLE]

since the expected values of an Ito integral and $u(t)$ are zero. If we are perfectly rigorous, Fubini’s theorem can’t be applied with delta distributions; nonetheless, the calculation can be repeated in complete analogy with sinc functions instead. Sinc functions appear as the Fourier transforms of finite rectangular windows, and are unavoidable in actual measurements131313Since every real measurement will be a finite time measurement. Therefore, the measured signal of an arbitrary $f(t)$ will in fact be $f(t)\cdot\Pi(t)$ ($\Pi(t)$ being a rectangular window), and by the convolution theorem the Fourier transform of the measurement will be $F(\omega)*\text{sinc}(\omega)$. This avoids the infinite values from the delta distribution..

Thus, defining the SNR as the ratio between expected driving signal amplitude and expected noise amplitude in the frequency domain, we find that the value goes to infinity. A more accurate analysis should compare not the ratio of the expected values but the ratio of the signal and noise random variables themselves (which, to a good approximation, should follow a Cauchy distribution141414The Cauchy distribution appears as the ratio of two normally distributed random variables and has undefined (i.e., going to infinity) statistical moments.). However, this section is enough to show that if more information about the driving force is known (in this case, the fact that the force is sinusoidal), we don’t need to restrict ourselves to the sensitivity described in eq. (7) and more intelligent approaches might exist.

3 Accelerometers

In the last section we saw the expressions for the sensitivity and response of the harmonic oscillator when subjected to an external driving. In a situation where the driving acts on the housing of the oscillator instead of on the oscillator itself, the equation of motion for the harmonic oscillator is modified to

[TABLE]

where $y(t)$ is the motion of the housing. By rewriting the equation in terms of $z\equiv x-y$, which is the quantity that we will measure, we obtain the equation

[TABLE]

Therefore, the transfer function for $y(t)$, as compared with a regular harmonic oscillator, will take the modified expression

[TABLE]

Note that, for large values of $\omega$, $H(\omega)\simeq 1$. Therefore, for frequencies significantly above resonance, the accelerometer has a flat response and behaves approximately as a dirac delta:

[TABLE]

so

[TABLE]

In this range of operation, the oscillator behaves as a seismometer: it measures the displacement of the housing (see Fig. 1).

To understand the accelerometer regime, lets assume we have an arbitrary housing motion $y(t)$. The response in terms of the power spectral densities will be

[TABLE]

But since $-\omega^{2}Y(\omega)=A(\omega)$, where $Y(\omega)$ and $A(\omega)$ are the Fourier transforms of $y(t)$ and the acceleration $a(t)=\ddot{y}(t)$, the response to an acceleration $\ddot{y}(t)$ will be

[TABLE]

In other words, the response of $z(t)$ to an acceleration of the housing is, except for the missing $1/m^{2}$ factor, the same as the response of a force acting on the harmonic oscillator, as plotted in Fig. 2. Now, the response of the harmonic oscillator is flat at low frequencies (i.e., below resonance)

[TABLE]

and this is the frequency band where most of the accelerometers work. Like in seismometers, this is a relevant point: working in a regime where the impulse response is approximately a Dirac delta is necessary when the driving force isn’t restricted to a small bandwidth (otherwise, to recover $a(t)$ out of $z(t)$, some sort of deconvolution may be needed, which is an ill-posed problem that should be avoided if possible).

Assuming the housing mass is large enough, the stochastic component of $y(t)$ due to Brownian noise can be neglected. However, Brownian noise will still drive the harmonic oscillator itself151515Note that an important approximation is being done here: we assume that the acceleration leaves the statistical properties of Brownian noise unchanged. For large accelerations this approximation may not hold.

[TABLE]

Therefore, repeating the analysis of the force sensing section we obtain

[TABLE]

where $H(\omega)$ is the response of the harmonic oscillator including the $1/m^{2}$ factor. Thus, neglecting the measurement noise, $\text{SNR}=\sqrt{\frac{S_{aa}m^{2}}{\sigma_{1}^{2}}}$ so

[TABLE]

This means that larger masses will push down the value of the minimum detectable accelerations.

{acknowledgement}

The author acknowledges financial support from the European Research Council through grant QnanoMECA (CoG - 64790), Fundació Privada Cellex and the Spanish Ministry of Economy and Competitiveness through the Severo Ochoa Programme for Centres of Excellence in R$\&$D (SEV-2015-0522), grant FIS2016-80293-R. The author acknowledges useful discussions and related work5 with the levitodynamics team of the Plasmon Nano-optics (Romain Quidant) group, at ICFO.

4 Supplementary material

4.1 Linear time-invariant systems

A linear constant-coefficient ordinary differential equation, defined by a linear polynomial in the unknown function $x(t)$ and its derivatives, takes the form

[TABLE]

where $f(\ldots)$ is a linear function and we added a non-homogeneous term $g(t)$. If we now Fourier transform both sides of the equation, we get

[TABLE]

where $p(i\omega)$ is the characteristic polynomial of the differential equation and $X(\omega)$ and $G(\omega)$ are the Fourier transforms of $x(t)$ and $g(t)$ respectively. Therefore

[TABLE]

where we defined the frequency response (or transfer function) $H(\omega)=\frac{1}{p(i\omega)}$. By the convolution theorem, if we apply the inverse Fourier transform we obtain

[TABLE]

[TABLE]

Here, $h(t)$ is known as the impulse response of the system, and is used to find $x(t)$ for an arbitrary driving $g(t)$.

4.2 The harmonic oscillator

The harmonic oscillator is a 2nd order constant coefficient linear ODE that can be used as a first approximation or a paradigmatic model of most simple oscillatory systems. In the most general case, starting from Newton’s second law, the equation reads

[TABLE]

where $m\Gamma\dot{x}$ is a damping force, $kx$ is a restoring force161616For now we are ignoring the fact that a damping force leads to a Brownian random force, due to the fluctuation-dissipation theorem. Since this force scales with $\sqrt{m}$, the randomness doesn’t need to be taken into account for large oscillator masses, but will be relevant in micro and nanoparticles. (of arbitrary origin) and $F(t)$ a driving force.

4.2.1 Conserved quantities

When $\Gamma=F(t)=0$, the system is Hamiltonian, and its energy takes the form

[TABLE]

where $\omega_{0}=\sqrt{\frac{k}{m}}$; it can be proved that this is the only conserved quantity of the system (Hamiltonian systems with $2n$ degrees of freedom can have, at most, $n$ constants of motion). The equality between energy and Hamiltonian is justified because $H$ doesn’t explicitly depend on $t$.

The first summand is the kinetic energy, and the second the (restoring force) potential energy. When feedback is introduced to the system, $H\equiv H(t)$, energy conservation cannot be taken for granted, as in general feedback can lead to cooling, heating, chaotic motion, etc.

4.2.2 Deterministic solution

Damped harmonic oscillator: If $F(t)$ is zero, we define the natural frequency $\omega_{0}=\sqrt{\frac{k}{m}}$, the damping ratio

[TABLE]

and the quality factor

[TABLE]

The quality factor can be understood in different (but mostly equivalent) manners. One way, that will be discussed later, is as the number of “coherent” oscillations of the system. A different definition is as the fraction of the energy $E$ stored in the system versus the energy dissipated $\Delta E$ in a period $\tau_{0}$ of the oscillation, as

[TABLE]

The last approximation comes from an order one Taylor expansion, so when $\Gamma$ is small one recovers the previous quality factor definition. Using these parameters, the ODE now takes the form

[TABLE]

The value of the damping ratio $\zeta$ critically determines the behaviour of the system. A damped harmonic oscillator can be

• •

Overdamped, $Q\leq 0.5$. The system exponentially decays to zero without oscillating. The case $Q=0.5$ is usually called critically damped and is the boundary between oscillation and no oscillation.

• •

Underdamped, $Q>0.5$. The system oscillates at $\omega_{1}=\omega_{0}\sqrt{1-Q^{2}/4}$, so the smaller the damping the closer the oscillation frequency to the natural frequency of the oscillator. The sinusoid has an exponential decay of $\lambda=\frac{\omega_{0}}{2Q}=\frac{\Gamma}{2}$.171717The characteristic time is $\tau=1/\lambda$

Driven harmonic oscillator: In the case of a sinusoidal driving force:

[TABLE]

where $F_{0}$ is the driving amplitude and $\omega$ is the driving frequency for a sinusoidal driving mechanism181818This type of system appears in AC driven inductor-capacitor systems and in driven spring systems having internal mechanical resistance or external, like a particle in an optical trap.

The general solution is a sum of a transient term that depends on initial conditions, and a steady state that is independent of initial conditions and depends only on the driving amplitude $F_{0}$, driving frequency $\omega$, undamped angular frequency $\omega_{0}$, and $\Gamma$. The transient solutions are the same as the unforced ($F_{0}=0$) harmonic oscillator and represent the systems response to other events that occurred previously. However, they typically die out rapidly enough that they can be ignored.

The steady-state solution is proportional to the driving force with an induced phase change of $\phi$:

[TABLE]

where

[TABLE]

is the absolute value of the linear response function, and

[TABLE]

is the phase of the oscillation relative to the driving force, if the $\arctan(\cdot)$ value is taken to be between -180 degrees and 0 (that is, it represents a phase lag, for both positive and negative values of the $\arctan$ argument).

These last two expressions are obtained from the complex transfer function

[TABLE]

the power transfer function being

[TABLE]

For a particular driving frequency, called the resonance or resonant frequency

[TABLE]

the amplitude (for a given $F_{0}$) is maximum, and

[TABLE]

The last approximation only holds for large $Q$ factors; notice that the ratio between the power transfer function at $\omega_{r}$ and at zero is $H(\omega_{r})^{2}/H(0)^{2}\simeq Q^{2}$.

The resonance effect only occurs when $Q>\sqrt{2}/2$, i.e. for significantly underdamped systems. For strongly underdamped systems the value of the amplitude can become quite large near the resonance frequency.

We can calculate the half width half maximum (HWHM) by imposing $Z_{m}^{2}$ to be half of the value at resonance. The angular frequencies thus obtained are

[TABLE]

and by approximating $\sqrt{Q^{2}-\frac{1}{4}}\simeq Q$, $\omega_{r}\simeq\omega_{0}$ and applying the Taylor series expansion of the square root in $\omega_{0}\sqrt{1\pm\frac{1}{Q}}$ , we obtain

[TABLE]

This shows that increasing $Q$ (or, equivalently, decreasing $\Gamma$) reduces the width of the frequency response of the oscillator. Therefore, a high $Q$ is particularly important when trying to detect frequency shifts.

4.2.3 Stochastic driving forces

If the intrinsic (classical) randomness of the system is taken into account via the fluctuation-dissipation theorem, then the driving force will have a stochastic driving term $F(t)=\sigma\eta(t)$, with $\eta(t)$ a zero mean and unit standard deviation white noise191919This is not exactly true, as white noise doesn’t exist; a proper treatment requires the use of Ito calculus. However, it is still useful to think of it as white noise. and $\sigma=\sqrt{2k_{B}T\gamma}$ (obtained from a fluctuation-dissipation theorem). The full system can be solved in all generality, but it is useful to consider the “overdamped” regime first.

Overdamped regime: It is common to discard the second order term of the stochastic differential equation when it is “small” compared to the other terms202020This can be made rigorous as a perturbative problem, a regular perturbation problem. However, this is not always the case when the parameter is multiplying the highest order term of the equation. See “Singular perturbation problem” for more details. In this case the equation takes the form of an Ornstein-Uhlenbeck process

[TABLE]

with $\gamma$ the damping constant and $\delta=m\omega_{0}^{2}$ the restoring force. In Ito’s notation

[TABLE]

with $a=\frac{\delta}{\gamma}>0$, $b=\frac{\sigma}{\gamma}$.

The solution of this SDE is

[TABLE]

To see how the process will diffuse, we can calculate its variance process with Ito’s isometry:

[TABLE]

which, for very short times (i.e., applying a first order Taylor expansion), scales as

[TABLE]

The autocorrelation can also be obtained by using Ito’s isometry on the last expression

[TABLE]

When $|t-s|=0$, the autocorrelation equals the variance of $X_{t}$,212121Equivalently this can also be written as $\langle X_{t}^{2}\rangle$

[TABLE]

For small $t$, the variance takes the expression

[TABLE]

while for long times

[TABLE]

as is expected from the equipartition theorem. As $t$ and $s$ increase, the second exponential summand of the process autocorrelation becomes arbitrarily small. The remaining part is a function of $\tau\equiv t-s$ only; thus, we can apply the Wiener-Khinchin theorem to get an analytical expression of the power spectrum of the process:

[TABLE]

[TABLE]

where in the last expression we have reintroduced the mass, as it will help in identifying similarities between this and the full second order system power spectrum.

As expected by Parseval’s theorem, the integral of the power spectral density is

[TABLE]

recovering again the value expected by the equipartition theorem.

Observation: The expression of the power spectral density has two clear different regimes: at low frequencies, the $a^{2}$ term in the denominator dominates and the spectrum is almost flat. At large $f$, $a^{2}$ is negligible and the other term dominates. In a log-log scale representation, the power spectrum looks like two straight lines (first an horizontal line and then a decreasing line); the frequency at which the behaviour changes is known as the corner frequency, $\omega_{c}=a=\frac{\delta}{\gamma}=\frac{\omega_{0}^{2}}{\Gamma}$.

Full 2nd order equation: The full equation takes the expression

[TABLE]

Performing a change of variables, $x=e^{\frac{-\gamma}{2m}t}x_{1}$, we get

[TABLE]

thus eliminating $\dot{x}$. Setting $a=\frac{k}{m}-\frac{\gamma^{2}}{4m^{2}}$, $b=\frac{\sigma}{m}$ and rewriting the equation as a first order linear system with

[TABLE]

where $v_{2}=\dot{x_{2}}$, we get in Ito’s notation

[TABLE]

The solution of a linear homogeneous SDE is

[TABLE]

where $A(t)$ is the (generally vector) coefficient of $X$. For this SDE a fundamental matrix solution of the associated homogeneous noise-free system is

[TABLE]

The determinant of this matrix is 1, so its inverse matrix will be

[TABLE]

and hence we can solve the complete system. We are interested in the first component of $X$, the position (as we will be calculating the PSD of the trajectory of the particle)

[TABLE]

Finally, $x_{1}(t)=e^{\frac{\gamma t}{2m}}x(t)$, so

[TABLE]

We see that, after a transient time, only the term depending on $\mathrm{d}W_{r}$ remains, so the first moment of the process is zero. Now, applying Ito’s isometry as before to calculate the covariance we get

[TABLE]

[TABLE]

This is a quite uninteresting calculation222222The results are verified with Mathematica. I have the file in ./Simulations/Mathematica. As in the one-dimensional Ornstein-Uhlenbeck process, one gets a term which only depends on the difference $|t-s|=\tau$, and another which is multiplied by $e^{-|t+s|}$, that quickly vanishes. Setting $t=s$ we recover the variance,

[TABLE]

where we have defined $\omega_{1}=\sqrt{\omega_{0}^{2}-(\Gamma/2)^{2}}=\sqrt{a}$ and used the normalized damping constant $\Gamma=\frac{\gamma}{m}$. As in the overdamped case,

[TABLE]

as expected from the equipartition theorem. However, more interesting is the behaviour at short times,

[TABLE]

Observation: this value of the variance assumes that $x_{0}$, $v_{0}$, initial conditions of the problem, are known. This may not be the case for a real experiment in the lab. For instance, assume we have a particle in an optical tweezer, and we want to calculate the mean square displacement (MSD),

$\mathbb{E}[(x(t)-x(0))^{2}]=\mathbb{E}[x^{2}(t)+x^{2}(0)-2x(0)x(t)].$

Then we distinguish two cases:

$x(0)$ is known, e.g. $x(0)=0$. Then

$\mathbb{E}[x^{2}(t)+x^{2}(0)-2x(0)x(t)]=\mathbb{E}[x^{2}(t)]=\mathcal{O}(t^{3}).$

$x(0)$ is not known (or is uncontrolled), and we do the ensemble average. Then

$\displaystyle\mathbb{E}[x^{2}(t)+x^{2}(0)-2x(0)x(t)]$ $\displaystyle=\mathbb{E}[x^{2}(t)]+\mathbb{E}[x^{2}(0)]-2\mathbb{E}[x(0)x(t)]$

$\displaystyle=\frac{2k_{B}T}{m\omega_{0}^{2}}-2R(t)=\mathcal{O}(t^{2})$

We see, therefore, that the leading term of the series expansion of the MSD is different in each case. This is an important distinction that is not clear in Toncang’s Li Science paper232323R. Rica contributed to this.

If in the autocorrelation integral we keep only the term with $\tau$ dependency (since the term multiplied by $e^{-|t+s|}$ will quickly decay), after some simplification one gets

[TABLE]

Worth mentioning is the value of the time constant of the exponential factor, $\frac{2}{\Gamma}$: this value is proportional to the time needed for the autocorrelation to get below a given threshold242424Which is arbitrarily set, otherwise known as losing the coherence of the oscillation. A usual criterion is 3 time constants, with the autocorrelation decreasing to below $e^{-3}<5\%$. The number of oscillations of the system during this time is

[TABLE]

Thus, the quality factor can be understood as the number of oscillations of the system under the presence of Brownian noise before the autocorrelation gets below 5$\%$252525The number of “coherent” oscillations.

From the expression of the autocorrelation we see that $R(t,\tau)=R(\tau)$: therefore, the process is wide-sense stationary and the conditions to apply the Wiener-Khinchin theorem are satisfied. The Fourier transform of this autocorrelation function is the power spectral density

[TABLE]

which, after replacing the variables and some rearranging262626See, again, the Mathematica file “mathem_psd_calculations” takes the simpler and more familiar expression

[TABLE]

where we have replaced the unitary ordinary frequency Fourier transform (in terms of $f$) by the non-unitary angular frequency Fourier transform. Thus, it is readily seen that this is the frequency response of a harmonic oscillator to a flat-spectrum (white noise) driving force.

As in the overdamped case,

[TABLE]

as expected by the equipartition theorem.