Skip to content

Commit

Permalink
suggestions: theory
Browse files Browse the repository at this point in the history
  • Loading branch information
rpauszek authored and JoepVanlier committed Dec 9, 2024
1 parent 35cd56c commit 814f1e0
Show file tree
Hide file tree
Showing 7 changed files with 84 additions and 49 deletions.
42 changes: 27 additions & 15 deletions docs/theory/force_calibration/active.rst
Original file line number Diff line number Diff line change
Expand Up @@ -126,9 +126,9 @@ If we define a factor :math:`c` by which the velocity is reduced, we obtain the
.. math::
\begin{align}
R_{d, corrected} & = c R_d\\
R_{f, corrected} & = \frac{1}{c} R_f\\
\kappa_{corrected} & = \frac{1}{c^2}\kappa
R_{d\mathrm{, corrected}} & = c R_d\\
R_{f\mathrm{, corrected}} & = \frac{1}{c} R_f\\
\kappa_\mathrm{corrected} & = \frac{1}{c^2}\kappa
\end{align}
Where :math:`R_d` is the displacement sensitivity, :math:`R_f` is the force sensitivity and :math:`\kappa` is the stiffness.
Expand All @@ -146,7 +146,9 @@ Filling in the maximal velocity we expect during the oscillation, we find the fo
Re = \frac{\rho u L}{\eta} = 2 \pi f A d \frac{\rho}{\eta}
Here :math:`\rho` refers to the fluid density, :math:`u` the characteristic velocity, :math:`L` the characteristic length scale, :math:`\eta` the viscosity, :math:`f` the oscillation frequency, :math:`A` the oscillation amplitude and :math:`d` the bead diameter.
Here :math:`\rho` refers to the fluid density, :math:`u` the characteristic velocity, :math:`L` the
characteristic length scale, :math:`\eta` the viscosity, :math:`f` the oscillation frequency, :math:`A`
the oscillation amplitude and :math:`d` the bead diameter.
For microfluidic flow, this value is typically much smaller than `1`.

In this limit, the Navier-Stokes equation describing fluid flow reduces to the following expressions:
Expand All @@ -160,32 +162,41 @@ In this limit, the Navier-Stokes equation describing fluid flow reduces to the f
Here :math:`\eta` is the viscosity, :math:`p` is the pressure and :math:`v` is the fluid velocity.
Creeping flow is far removed from every day intuition as it equilibrates instantaneously.
The advantage of this is that for sufficiently low frequencies, the correction factor can be based on the correction factor one would obtain for a steady state constant flow.
The advantage of this is that for sufficiently low frequencies, the correction factor can be based on
the correction factor one would obtain for a steady state constant flow.

For two beads aligned in the flow direction, we can use the analytical solution presented in :cite:`stimson1926motion`.
This model uses symmetry considerations to solve the creeping flow problem for two solid spheres moving at a constant velocity parallel to their line of centers.
This model uses symmetry considerations to solve the creeping flow problem for two solid spheres moving
at a constant velocity parallel to their line of centers.
We denote the correction factor obtained from this model as :math:`c_{\|}`.
This correction factor is given by the ratio of the drag coefficient by the drag coefficient one would expected from a single bead in creeping flow (:math:`3 \pi \eta d v`).
For beads aligned perpendicular to the flow direction, we use a model from :cite:`goldman1966slow`, which we denote as :math:`c_{\perp}`.
This correction factor is given by the ratio of the drag coefficient by the drag coefficient one would
expected from a single bead in creeping flow (:math:`3 \pi \eta d v`).
For beads aligned perpendicular to the flow direction, we use a model from :cite:`goldman1966slow`,
which we denote as :math:`c_{\perp}`.

From the derivations in these papers, it follows that the correction factors obtained depend on the bead diameter(s) :math:`d` and distance between the beads :math:`l`.
From the derivations in these papers, it follows that the correction factors obtained depend on the
bead diameter(s) :math:`d` and distance between the beads :math:`l`.
For equally sized beads, this dependency is a function of the ratio of the distance between the beads over the bead diameter.

Considering the linearity of the equations that describe creeping flow :cite:`goldman1966slow`, we can combine the two analytical solutions by decomposing the incoming velocity (in the direction :math:`\vec{e}_{osc}`) into a velocity perpendicular to the bead-to-bead axis :math:`\vec{e}_{\perp}` and a velocity component aligned with the bead-to-bead axis :math:`\vec{e}_{\|}`.
Considering the linearity of the equations that describe creeping flow :cite:`goldman1966slow`, we can
combine the two analytical solutions by decomposing the incoming velocity (in the direction :math:`\vec{e}_{osc}`)
into a velocity perpendicular to the bead-to-bead axis :math:`\vec{e}_{\perp}` and a velocity component
aligned with the bead-to-bead axis :math:`\vec{e}_{\|}`.

.. math::
\begin{align}
v_{\|} & = (\vec{e}_{\|} \cdot\vec{e}_{osc}) c_{\|}\\
v_{\perp} & = (\vec{e}_{\perp} \cdot \vec{e}_{osc}) c_{\perp}
v_{\|} & = (\vec{e}_{\|} \cdot\vec{e}_\mathrm{osc}) c_{\|}\\
v_{\perp} & = (\vec{e}_{\perp} \cdot \vec{e}_\mathrm{osc}) c_{\perp}
\end{align}
This provides us with contributions for each of those axes, but we still need to project this back to the oscillation axis (since this is where we measure our amplitude).
This provides us with contributions for each of those axes, but we still need to project this back
to the oscillation axis (since this is where we measure our amplitude).
We can calculate our desired hydrodynamic correction factor as:

.. math::
c_{total} = v_{\|} (\vec{e}_{\|} \cdot \vec{e}_{osc}) + v_{\perp} (\vec{e}_{\perp} \cdot \vec{e}_{osc})
c_\mathrm{total} = v_{\|} (\vec{e}_{\|} \cdot \vec{e}_\mathrm{osc}) + v_{\perp} (\vec{e}_{\perp} \cdot \vec{e}_\mathrm{osc})
The response of this combined model for equally sized beads can be calculated as follows::

Expand All @@ -201,5 +212,6 @@ The response of this combined model for equally sized beads can be calculated as

.. image:: figures/correction_factor.png

Here, when providing only a horizontal distance recovers the Stimson model :cite:`stimson1926motion`, while a vertical displacement recovers the Goldman model :cite:`goldman1966slow`.
Here, when providing only a horizontal distance recovers the Stimson model :cite:`stimson1926motion`,
while a vertical displacement recovers the Goldman model :cite:`goldman1966slow`.
To find out more about how to use these correction factors, please refer to the :ref:`tutorial<bead_bead_tutorial>`.
20 changes: 13 additions & 7 deletions docs/theory/force_calibration/diode.rst
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,18 @@ In reality, our measurement is affected by two processes:
This second factor depends on the type of measurement device being used.
Typical position sensitive detectors are made of silicon.
Such a detector has a very high bandwidth for visible light (in the MHz range).
Unfortunately, the bandwidth is markedly reduced for the near infra-red light of the trapping laser :cite:`berg2003unintended,berg2006power`.
Unfortunately, the bandwidth is markedly reduced for the near infra-red light of the trapping laser
:cite:`berg2003unintended,berg2006power`.
This makes it less sensitive to changes in signal at high frequencies.

Why is the bandwidth limited?
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The high bandwidth of visible light detection of a silicon photodiode is achieved when incoming photons are absorbed in the so-called depletion layer of the diode.
The high bandwidth of visible light detection of a silicon photodiode is achieved when incoming photons
are absorbed in the so-called depletion layer of the diode.
Unfortunately, silicon has an increased transparency at the near infra-red wavelength of the trapping laser.
The result of this is that light penetrates deeper into the substrate of the diode, where it generates charge carriers in a different region of the diode.
The result of this is that light penetrates deeper into the substrate of the diode, where it generates
charge carriers in a different region of the diode.
These charge carriers then have to diffuse back to the depletion layer, which takes time.
As a result, a fraction of the signal is detected with a much lower bandwidth.

Expand All @@ -41,19 +44,22 @@ This model is characterized by two numbers whose values depend on the incident l
High corner frequencies
^^^^^^^^^^^^^^^^^^^^^^^

In literature, the diode parameters are frequently estimated simultaneously with the calibration data :cite:`berg2003unintended,hansen2006tweezercalib,berg2006power,tolic2006calibration,tolic2004matlab,berg2004power`.
In literature, the diode parameters are frequently estimated simultaneously with the calibration data
:cite:`berg2003unintended,hansen2006tweezercalib,berg2006power,tolic2006calibration,tolic2004matlab,berg2004power`.
Unfortunately, this can cause issues when calibrating at high powers.

Recall that the physical spectrum is characterized by a corner frequency `fc`, and diffusion constant `D`.
The corner frequency depends on the laser power and bead size (smaller beads resulting in higher corner frequencies).
The parasitic filtering effect also has a corner frequency (`f_diode`) and depends on the incident intensity :cite:`berg2003unintended`.

When these two frequencies get close, they cannot be estimated from the power spectrum reliably anymore.
When these two frequencies are similar, they cannot be estimated from the power spectrum reliably anymore.
The reason for this is that the effects that these parameters have on the power spectrum becomes very similar.
When working with small beads or at high laser powers, it is important to verify that the corner frequency `fc` does not approach the frequency of the filtering effect `f_diode`.
When working with small beads or at high laser powers, it is important to verify that the corner frequency `fc`
does not approach the frequency of the filtering effect `f_diode`.

Sometimes, the parameters of this diode have been characterized independently.
In that case, the arguments `fixed_diode` and `fixed_alpha` can be passed to :func:`~lumicks.pylake.calibrate_force()` to fix these parameters to their predetermined values resolving this issue.
In that case, the arguments `fixed_diode` and `fixed_alpha` can be passed to :func:`~lumicks.pylake.calibrate_force()`
to fix these parameters to their predetermined values, resolving this issue.
For more information on how to achieve this with Pylake, please refer to the :ref:`diode calibration tutorial<diode_tutorial>`.

Mathematical background
Expand Down
37 changes: 22 additions & 15 deletions docs/theory/force_calibration/fitting.rst
Original file line number Diff line number Diff line change
Expand Up @@ -56,8 +56,8 @@ frequency :cite:`berg2004power`.

.. math::
\bar{f} &= \frac{1}{n_b} \sum_{f \in block} f\\
\bar{P}_{meas} &= \frac{1}{n_b} \sum_{f \in block} P_{meas}(f)
\bar{f} &= \frac{1}{n_b} \sum_{f \in \mathrm{block}} f\\
\bar{P}_\mathrm{meas} &= \frac{1}{n_b} \sum_{f \in \mathrm{block}} P_\mathrm{meas}(f)
Setting the number of points per block too low results in a bias from insufficient averaging :cite:`berg2004power`.
Insufficient averaging would result in an overestimation of the force response :math:`R_f` and an
Expand All @@ -76,23 +76,24 @@ the spectral data points are approximately Gaussian distributed with standard de
\sigma(\bar{f}) = \frac{P(\bar{f})}{\sqrt{n_b}}
This means that regular weighted least squares (WLS) can take place.
This means that regular weighted least squares (WLS) can be used for fitting.
To ensure unbiased estimates in WLS, the data and squared weights must be uncorrelated.
However, we know that there is a known correlation between these which results
in a known bias in the estimate for the diffusion constant that can be corrected for after
However, there is a known correlation between these which results
in a known bias in the estimate for the diffusion constant that can be corrected after
fitting :cite:`norrelykke2010power`:

.. math::
D_{corrected} = D_{wls} \frac{n_b}{n_b + 1}
D_\mathrm{corrected} = D_\mathrm{wls} \frac{n_b}{n_b + 1}
.. _noise_floor:

Noise floor
^^^^^^^^^^^

When operating at very low powers (and by extension corner frequencies), a noise floor may be visible at high frequencies.
It is important to ensure that the upper limit of the fitting range does *not* include the noise floor as it is not taken into account in the calibration model.
It is important to ensure that the upper limit of the fitting range does *not* include the noise floor
as it is not taken into account in the calibration model.
In the dataset below, we can see the effect of the noise floor::

lk.download_from_doi("10.5281/zenodo.7729823", "test_data")
Expand Down Expand Up @@ -151,33 +152,39 @@ Note that when we have a diode calibration, excluding the noise floor becomes ev

.. image:: figures/noise_floor_fixed_diode.png

The reason the effect of the noise floor on the calibration parameters is more pronounced is because with the fixed diode model, the model is less free to adjust the diode model to mitigate its impact.
The reason the effect of the noise floor on the calibration parameters is more pronounced is because
with the fixed diode model, the model is not free to adjust the diode parameters to mitigate its impact.
As a result, the model uses the corner frequency in an attempt to capture the shape of the noise floor (strongly biasing the result).

Noise peaks
^^^^^^^^^^^

Optical tweezers are precision instruments.
Despite careful determination and elimination of noise sources, it is not always possible to exclude all potential sources of noise.
One downside of weighted least squares estimation, is that it is very sensitive to outliers.
Despite careful determination and elimination of noise sources, it is not always possible to exclude all potential sources of noise,
which can manifest in the power spectrum as spurious peaks.
One downside of weighted least squares estimation, is that it is very sensitive to such outliers.
It is therefore important to either exclude noise peaks from the data prior to fitting or use :ref:`robust fitting<robust_fitting>`.
Noise peaks are always excluded prior to blocking to minimize data loss.
Noise peaks should always be excluded prior to blocking to minimize data loss.

.. _goodness_of_fit:

Goodness of fit
---------------

When working with the Gaussian error model, we can calculate a goodness of fit criterion.
When sufficient blocking has taken place, the sum of squared residuals that is being minimized during the fitting procedure is distributed according to a chi-squared distribution characterized by :math:`N_{\mathit{dof}} = N_{\mathit{data}} - N_{\mathit{free}}` degrees of freedom.
Here :math:`N_{\mathit{data}}` corresponds to the number of data points we fitted (after blocking) and :math:`N_{\mathit{free}}` corresponds to the number of parameters we fitted.
When sufficient blocking has taken place, the sum of squared residuals that is being minimized during
the fitting procedure is distributed according to a chi-squared distribution characterized by
:math:`N_{\mathrm{dof}} = N_{\mathrm{data}} - N_{\mathrm{free}}` degrees of freedom.
Here :math:`N_{\mathrm{data}}` corresponds to the number of data points we fitted (after blocking) and
:math:`N_{\mathrm{free}}` corresponds to the number of parameters we fitted.
We can use the value we obtain to determine how unusual the fit error we obtained is.

.. math::
\mathrm{support} = 100 P(x > \chi_{\mathit{obtained}}^2) = 100 \int_{\chi_{\mathit{obtained}}^2}^{\infty} \chi^2_{N_{\mathit{dof}}}(x) dx
\mathrm{support} = 100 P(x > \chi_{\mathrm{obtained}}^2) = 100 \int_{\chi_{\mathrm{obtained}}^2}^{\infty} \chi^2_{N_{\mathrm{dof}}}(x) dx
The support or backing is the probability that a repetition of the measurement that produced the data we fitted to will, after fitting, produce residuals whose squared sum is greater than the one we initially obtained.
The support or backing is the probability that a repetition of the measurement that produced the data
we fitted will, after fitting, produce residuals whose squared sum is greater than the one we initially obtained.
More informally, it represents the probability that a fit error at least this large should occur by chance.

Support less than 1% warrants investigating the residuals for any trend in the residuals.
2 changes: 1 addition & 1 deletion docs/theory/force_calibration/force_calibration.rst
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ These are the time intervals we associated with the broader, slow movements of t
If we think about a bead moving due to random collisions, then we can expect that the bead will move
more in these longer time intervals. This is why in the power spectrum of free diffusion, we see a
lot more energy concentrated at these low frequencies, while the rapid jiggles at higher frequency
contribute far less. The amplitude of this power spectrum is related to :math:`D`.
contribute far less. The amplitude of this power spectrum is related to :math:`D`.

.. image:: figures/sim_trap_opt.gif
:nbattach:
Expand Down
15 changes: 10 additions & 5 deletions docs/theory/force_calibration/hyco.rst
Original file line number Diff line number Diff line change
Expand Up @@ -11,9 +11,10 @@ While the idealized model of the bead motion is sometimes sufficiently accurate,
it neglects inertial and hydrodynamical effects of the fluid and bead(s).

The frictional forces applied by the viscous environment to the bead are proportional to the bead's velocity relative to the fluid.
The idealized model is based on the assumption that the bead's relative velocity is constant.
There is no dynamical change of the fluid motion around the bead.
In reality, when the bead moves through the fluid, the frictional force between the bead and the fluid depends on the past motion, since that determines the fluid's current motion.
The idealized model is based on the assumption that the bead's relative velocity is constant;
there is no dynamical change of the fluid motion around the bead.
In reality, when the bead moves through the fluid, the frictional force between the bead and the fluid
depends on the past motion, since that determines the fluid's current motion.
For a stochastic process such as Brownian motion, constant motion is not an accurate assumption.
In addition, the bead and the surrounding fluid have their own mass and inertia, which are also neglected in the idealized model.

Expand All @@ -33,9 +34,13 @@ sizes and higher trap powers the differences can be substantial.
Fast sensor measurement
^^^^^^^^^^^^^^^^^^^^^^^

When fitting a power spectrum, one may ask the question, so why does the fit look good if the model is bad?
.. note::

The following section only applies to instruments which contain fast PSDs.

When fitting a power spectrum, one may ask the question: "Why does the fit look good if the model is bad?"
The answer to this lies in the model that is used to capture the :ref:`parasitic filtering effect<diode_theory>`.
When the parameters of this model are estimated, what can happen is that they "hide" the mis-specification of the model.
When the parameters of this model are estimated, they can "hide" the mis-specification of the model.

Fast detectors have the ability to respond much faster to incoming light resulting in no visible filtering effect in the frequency range we are fitting.
This means that for a fast detector, we do not need to include such a filtering effect in our model and see the power spectrum for what it really is.
Expand Down
4 changes: 2 additions & 2 deletions docs/theory/force_calibration/passive.rst
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ To convert the parameters obtained from this spectral fit to a trap stiffness, t
\kappa = 2 \pi \gamma_0 f_c \tag{$\mathrm{N/m}$}
Here :math:`\kappa` then represents the estimated trap stiffness.
where :math:`\kappa` is the estimated trap stiffness.

We can calibrate the position by considering the diffusion of the bead:

Expand Down Expand Up @@ -76,7 +76,7 @@ It is especially the latter that results in large calibration errors when mis-sp
\begin{align}
\kappa = 2 \pi \gamma(T) f_c &\propto& \eta(T)\\
R_d = \sqrt{\frac{kT}{\gamma(T)D_{volts}}} &\propto& \sqrt{T / \eta(T)}\\
R_d = \sqrt{\frac{kT}{\gamma(T)D_\mathrm{volts}}} &\propto& \sqrt{T / \eta(T)}\\
R_f = R_d \kappa &\propto& \sqrt{T \eta(T)}
\end{align}
Expand Down
Loading

0 comments on commit 814f1e0

Please sign in to comment.