docs: document the coupling factor

changelog.md

@@ -1,5 +1,11 @@
 # Changelog
 
+## v1.4.0 | t.b.d.
+
+#### New features
+
+* Added model to correct for bead-bead coupling when using active calibration deep in bulk with two beads. See [`tutorial`](https://lumicks-pylake.readthedocs.io/en/latest/tutorial/force_calibration.html#active-calibration-with-two-beads-far-away-from-the-surface) and [`theory`](https://lumicks-pylake.readthedocs.io/en/latest/theory/force_calibration/active.html#bead-bead-coupling) for more information.
+
 ## v1.3.2 | t.b.d.
 
 #### Bug fixes

docs/refs.bib

@@ -456,3 +456,36 @@ @article{Rabiner1989HMM
   year={1989},
   publisher={IEEE}
 }
+
+@article{stimson1926motion,
+  title={The motion of two spheres in a viscous fluid},
+  author={Stimson, Margaret and Jeffery, George Barker},
+  journal={Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character},
+  volume={111},
+  number={757},
+  pages={110--116},
+  year={1926},
+  publisher={The Royal Society London}
+}
+
+@article{goldman1966slow,
+  title={The slow motion of two identical arbitrarily oriented spheres through a viscous fluid},
+  author={Goldman, AJ and Cox, RG and Brenner, H},
+  journal={Chemical Engineering Science},
+  volume={21},
+  number={12},
+  pages={1151--1170},
+  year={1966},
+  publisher={Elsevier}
+}
+
+@misc{the_fenics_project_developers_2023_10432590,
+  author={The FEniCS Project Developers},
+  title={FEniCS/dolfinx: v0.7.3},
+  month={dec},
+  year={2023},
+  publisher={Zenodo},
+  version={v0.7.3},
+  doi={10.5281/zenodo.10432590},
+  url={https://doi.org/10.5281/zenodo.10432590}
+}

docs/theory/force_calibration/active.rst

@@ -91,3 +91,83 @@ hydrodynamically correct model for the peak:
 We can also include a distance to the surface like before. This results in an expression for the drag
 coefficient :math:`\gamma` that depends on the distance to the surface which is given by the same
 equations as listed in the section on the :doc:`hydrodynamically correct model<hyco>`.
+
+Bead-bead coupling
+------------------
+
+The active calibration method presented in the previous sections relies on oscillating the nanostage with a known amplitude and frequency.
+The fluid in the flow-cell follows the stage motion.
+This in turn exerts a drag on the bead that leads to a sinusoidal displacement of the bead from the trap center.
+The amplitude of the detected displacement (measured in Volts) and the stage amplitude are then quantified.
+From the stage amplitude (measured in microns, since the stage position is calibrated) an expected bead displacement is calculated.
+
+When using two beads, the flow field around the beads is reduced (because the presence of the additional bead slows down the fluid).
+The magnitude of this effect depends on the bead diameter, distance between the beads and their orientation with respect to the fluid flow.
+Streamlines for some bead configurations are shown below (simulated using FEniCSx :cite:`the_fenics_project_developers_2023_10432590`).
+
+.. image:: figures/streamlines.png
+  :nbattach:
+
+As a result, the bead moves less than expected for a given stage motion.
+Since the displacement sensitivity (microns/V) is given by the ratio of the expected bead displacement (in microns) to detected displacement (in Volts) and we detected less displacement than expected (lower voltage amplitude), we obtain an artificially higher displacement sensitivity than expected.
+To correctly take this into account, we need to take into account what happens to the fluid around the beads.
+
+Considering the fluid velocity and viscosity, we can conclude that we typically operate in the regime where viscous effects are dominant (creeping flow).
+This can be validated by calculating the Reynolds number for the flow.
+Filling in the maximal velocity we expect during the oscillation, we find the following expression.
+
+.. math::
+
+    Re = \frac{\rho u L}{\mu} = 2 \pi f A d \frac{\rho}{\mu}
+
+Here :math:`\rho` refers to the fluid density, :math:`u` the characteristic velocity, :math:`L` the characteristic length scale and :math:`\mu` the viscosity.
+For microfluidic flow, this value is typically smaller than `1`.
+
+In this limit, the Navier-Stokes equation reduces to the following expressions:
+
+.. math::
+
+    \begin{align}
+    -\nabla \cdot ( \nabla u + pI ) & = 0\\
+    \nabla \cdot v & = 0
+    \end{align}
+
+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.
+
+For two beads aligned horizontally in creeping flow, 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.
+
+For beads aligned perpendicular to the flow direction, we use a model from :cite:`goldman1966slow`.
+
+Considering the linearity of the equations that describe creeping flow :cite:`goldman1966slow`, we can combine these two solutions by decomposing the incoming unit velocity :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}_{radial}`.
+
+.. math::
+
+    \begin{align}
+    v_{aligned} & = (\vec{e}_{radial} \cdot\vec{e}_{osc}) c_{aligned}\\
+    v_{perp} & = (\vec{e}_{perp} \cdot \vec{e}_{osc})  c_{perp}
+    \end{align}
+
+This provides us with a fractional velocity along those axes, but we still need to project this velocity back to the oscillation axis (since this is where we measure our amplitude):
+
+.. math::
+
+    c_{total} = v_{aligned} (\vec{e}_{radial} \cdot \vec{e}_{osc}) + v_{perp} (\vec{e}_{perp} \cdot \vec{e}_{osc})
+
+The response of this combined model for equally sized beads can be calculated as follows::
+
+    diameter = 1.0
+    l_d = np.arange(1.01, 8, 0.1) * diameter
+    zeros = np.zeros(l_d.shape)
+    plt.plot(l_d, lk.coupling_correction_2d(l_d, zeros, diameter, is_y_oscillation=False), label="horizontal alignment [Stimson et al]")
+    plt.plot(l_d, lk.coupling_correction_2d(zeros, l_d, diameter, is_y_oscillation=False), label="vertical alignment [Goldman et al]")
+    plt.plot(l_d, lk.coupling_correction_2d(l_d / np.sqrt(2), l_d / np.sqrt(2), diameter, is_y_oscillation=False), label="diagonal alignment")
+    plt.ylabel('Correction factor [-]')
+    plt.xlabel("l/D [-]")
+    plt.legend()
+    plt.savefig('Correction factor.png', dpi=100)
+
+.. 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`.

docs/tutorial/force_calibration.rst

@@ -261,8 +261,8 @@ Axial force calibration can be performed by specifying `axial=True`::
 
     force_model = lk.PassiveCalibrationModel(bead_diameter, distance_to_surface=5, axial=True)
 
-Active calibration
-------------------
+Active calibration with a single bead
+-------------------------------------
 
 Active calibration has a few benefits.
 When performing passive calibration, we base our calculations on a theoretical drag coefficient which depends on parameters that are only known with limited precision:
@@ -394,6 +394,42 @@ However, if we do not provide the height above the surface, we can see that the
     Consequently, when this model is selected, this parameter affects both passive and active calibration.
     For more information on this see the :doc:`theory section on force calibration</theory/force_calibration/force_calibration>` section.
 
+Active calibration with two beads far away from the surface
+-----------------------------------------------------------
+
+.. warning::
+
+    The implementation of the coupling correction models is still alpha functionality.
+    While usable, this has not yet been tested in a large number of different scenarios.
+    The API can still be subject to change *without any prior deprecation notice*!
+    If you use this functionality keep a close eye on the changelog for any changes that may affect your analysis.
+
+When performing active calibration, we get a smaller fluid velocity around the beads than expected when calibrating with two beads in a dual trap configuration.
+This leads to a smaller voltage readout than expected if there's no coupling and therefore a higher displacement sensitivity (microns per volt).
+Failing to take this into account results in a bias.
+Pylake offers a function to calculate a correction factor to account for the lower velocity around the bead.
+Appropriate coupling correction factors for oscillation in x can be calculated as follows::
+
+    factor = lk.coupling_correction_2d(dx=5.0, dy=0, bead_diameter=bead_diameter, is_y_oscillation=False)
+
+Here `dx` and `dy` represent the horizontal and vertical distance between the beads, while the parameter `bead_diameter` refers to the bead diameter.
+Note that all three parameters have to be specified in the same spatial unit (meters or micron).
+The final parameter `is_y_oscillation` indicates whether the stage was oscillated in the y-direction.
+
+The obtained correction factor can be used to correct the calibration factors::
+
+    Rd_corrected = factor * calibration["Rd"].value
+    Rf_corrected = calibration["Rf"].value / factor
+    stiffness_corrected = calibration["kappa"].value / factor**2
+
+To correct a force trace, simply divide it by the correction factor::
+
+    corrected_force1x = f.force1x / factor
+
+.. note::
+
+    This coupling model is only valid far away from the surface.
+
 Fast Sensors
 ------------
 

lumicks/pylake/__init__.py

@@ -31,6 +31,7 @@
     ActiveCalibrationModel,
     PassiveCalibrationModel,
     viscosity_of_water,
+    coupling_correction_2d,
 )
 from .force_calibration.power_spectrum_calibration import (
     fit_power_spectrum,