Optical tweezers-based velocimetry: a method to measure microscale unsteady flows

Abstract

In the study of micro-scale biological flows, velocimetry methods based on passive tracers, such as micro-PIV and micro-PTV, are well established to characterize steady flows. However, these methods become inappropriate for measuring unsteady flows of small amplitude, because, on these scales, the motion of passive tracers cannot be distinguished from Brownian motion. In this study, we use optical tweezers (OTs) in combination with Kalman filtering, to measure unsteady microscopic flows with high temporal accuracy. This method is referred to as optical tweezers-based velocimetry (OTV). The OTV method measures the nanometric displacements of a trapped bead, and predicts the instantaneous velocity of the flow by employing a Kalman filter. We discuss the accuracy of OTV in measuring unsteady flows with 1.5–70 $\upmu$m s$^{-1}$ amplitudes and 10–90 Hz frequencies. We quantify how the bead size and the laser power affect the velocimetry accuracy, and specify the optimal choices for the bead size and laser power to measure different unsteady flows. OTV accurately measures unsteady flows with amplitudes as small as 3–6 $\upmu$m s$^{-1}$. We compare the accuracy of OTV and micro-PTV, and characterize the flow regime for which OTV outperforms micro-PTV. We also demonstrate the robustness of OTV by measuring the unsteady flow created by the cilia of green alga Chlamydomonas reinhardtii, and comparing with numerical predictions based on Stokes equations. An open-source implementation of the OTV software in Matlab is available through the 4TU.Centre for Research Data.

Graphic abstract

1 Introduction

Accurate flow velocimetry measurements are crucial in a range of applications from engineering to physics and biology. Established velocimetry methods, developed to measure and visualize flows, include pitot tubes, hot-wire anemometry as well as optical methods relying on passive tracer particles, Laser Doppler Anemometry (LDA), Particle Tracking Velocimetry (PTV) and Particle Image Velocimetry (PIV) (Adrian et al. 2011; Lindken et al. 2009). Determining the appropriate velocimetry method strongly depends on the flow: the characteristic length scale L, the velocity magnitude U, and the timescale of the velocity variations $\tau$. On the micron scale, flow velocity measurements are widely preformed using methods based on passive tracer particles such as micro-PIV (Akbaridoust et al. 2018; Adhikari et al. 2016; Kowalczyk et al. 2007; Kondratieva et al. 2008; Yan et al. 2007; Guasto et al. 2010; Drescher et al. 2010; Williams and Wereley 2009), micro-PTV (Chen et al. 2014; Park and Kihm 2006), and micro-holographic-PTV (Lee et al. 2016; Kim and Lee 2008), see Fig. 1. These methods have been successfully implemented to measure flows in engineering (Kiddy et al. 2000), and biological systems (Poelma et al. 2008, 2010). In all of these velocimetry techniques, the flow velocity is deduced from the displacements of passive beads. On the micron scale, however, particle displacements are also affected by Brownian motion due to thermal fluctuations. The Brownian motion of such tracer particles limits the range of flows that such velocimetry methods can accurately measure.

Consider a flow with a characteristic velocity U, and with a characteristic timescale $\tau$ or frequency of the velocity variation $f \sim 1/\tau$. The bead displacement due to the advection of the particle by the flow scales with $\delta _{\mathrm{adv}} \sim U\tau$, and the bead displacement due to diffusion scales with $\delta _{\mathrm{diff}} \sim \sqrt{D\tau }$, where D is the diffusion coefficient of the tracer particle. The ratio of these timescales corresponds to the Péclet number ${\text {Pe}} \sim U\sqrt{\tau /D}$. The Péclet number is large when the magnitude of the flow is large or when the characteristic timescale for velocity fluctuations $\tau$ is large. In this limit, the displacements of passive tracers are dominated by flow advection, and velocimetry based on passive tracer is justified, see Fig. 1. However, tracer-based methods cannot accurately measure unsteady flows with low Péclet numbers, because the particle displacements due to the flow become comparable to the displacements due to thermal diffusion. This limit is reached when the velocity magnitude scales with $U \sim \sqrt{Df}$, shown as the red line in Fig. 1. This limit is reached for measurements of biological flows, such as those created by ciliated organisms (Wei et al. 2019). Here, we present and fully characterize a velocimetry method using optical tweezers and Kalman filtering, which achieves accurate velocimetry measurements below the diffusion limit.

Fig. 1

Schematic frequency–amplitude diagram for the flow velocimetry studies done at micro-scale. The red curve shows the diffusion limit of a free tracer particle, ${\text {Pe}}\sim$ 1, and the green curve shows the signal-to-noise ratio limit, SNR$^\prime \sim$ 1, for a trapped particle. For the definition of SNR$^\prime$, please see Sect. 4.1

Optical tweezers (OTs) trap beads in a focussed laser beam and are widely used to measure forces on the micrometric and nanometric scales (Marti and Hübner 2010). The high level of accuracy and temporal resolution achieved by OTs in measuring small forces are ideal for velocimetry at small scales. In the presence of a flow, the hydrodynamic force on the trapped bead can be measured with OTs to deduce the flow velocity. In fact, on micron scales, the hydrodynamic force on the bead corresponds to the Stokes drag, which is directly proportional to the flow velocity. OTs can, therefore, be readily used for flow velocimetry purposes. Previous studies have used OTs as a positioning tool to facilitate flow measurements by trapping a particle at a given location, before releasing it to use it as a passive tracer for micro-PTV (Di Leonardo et al. 2006; Leach et al. 2006; Knöner et al. 2005; Nève et al. 2008; Padgett and Di Leonardo 2011). This use of OTs enhances the spatial resolution of the measurements, but the accuracy and temporal resolution remain limited by the diffusion limit. OTs have also been used to measure the flow directly inside microfluidic devices (Eom et al. 2014; Almendarez-Rangel et al. 2018; Nemet and Cronin-Golomb 2002; Mushfique et al. 2008; Nedev et al. 2014). These studies have been limited to measuring the average flow velocity, when the timescale of the flow variations $\tau$ is large and therefore ${\text {Pe}} \gg 1$. Recently, OTs have been used to measure unsteady oscillating flows. Harlepp et al. (2017) have measured the low frequency oscillating flows in a study of the cardiovascular flow of a zebrafish embryo. The pulsing blood flow (roughly 2–3 Hz) was measured by trapping a blood cell. Köhler et al. (2016) used OTs to measure oscillating radial flows at 1–4 Hz, generated by the nonaxisymmetric rotation of a microrotor. In both studies (Harlepp et al. 2017; Köhler et al. 2016), the low frequencies and higher velocities correspond to ${\text {Pe}} \gg 1$, for which other techniques could be used.

However, in many biophysical studies, the timescale $\tau$ associated with the unsteadiness can be very small $\tau \sim$10–100 ms. One example is the quasi-periodic flows generated by cilia, which have frequencies ranging from 15 to 100 Hz (see the shaded area in Fig. 1). Previous velocimetry measurements around ciliated micro-organisms have used micro-PTV and have reported the rate of spatial decay of the average flow velocity (Francis et al. 2009; Olstad et al. 2019; Reiten et al. 2017). Characterizing the unsteady component of the flow with micro-PIV requires to perform time stamped correlation averaging, which is not always possible.

Here, we discuss the performance of Optical Tweezers-based Velocimetry (OTV) in measuring oscillatory flows of different magnitudes U and frequencies f. Our method uses an OT to trap a bead in the flow and uses Kalman filtering to deduce the instantaneous flow velocity from the bead position. The accuracy of the velocity measurements is assessed by imposing a range of pre-calibrated flow fields and comparing our velocimetry measurements with the pre-calibrated values. To do this, we use periodic oscillatory uniform flows of a wide range of U and f as pre-calibrated flow fields. The performance of OTV depends on the stiffness of the trap and the size of the trapped bead. We characterize both and discuss the optimal choices for the trap stiffness and the bead size depending on the flow characteristics (U, f). We quantify the accuracy of the OTV to measure unsteady flows with ${\text {Pe}}>1$ and ${\text {Pe}} \leqslant 1$. Our OTV method can accurately measure unsteady flows as small as 3–10 $\upmu$m s$^{-1}$ at frequencies of 10–100 Hz. It bears emphasis, that the OTV is versatile and robust in measuring any unsteady flows and is not limited to measuring strictly periodic flows. We demonstrate this through direct measurements of ciliary flows, which are pseudo-periodic flows with a broad frequency spectrum. This measurement also confirms the compatibility of OTV with biological flow measurements. This method has been used recently by Wei et al. (2019) to measure the unsteady flow around C. reinhardtii. The Matlab implementation of the Kalman filter by EKF/UKF toolbox (Hartikainen et al. 2011) is available online through the 4TU. Centre for Research Data (Ghoddoosi Dehnavi and Tam 2020).

Fig. 2

The simplified schematic of our experimental setup. The flow chamber is connected to a piezo stage, which is used to generate a known background flow. The inset shows two opposing forces acting on the trapped bead: the hydrodynamic force exerted by the flow, and the trapping force exerted by the OT on the bead

2 Experimental setup

2.1 Optical tweezers setup

In optical tweezers, a particle with a refractive index different from the surrounding fluid can be trapped at the focal point of a highly focused laser beam (Keen et al. 2007; Neuman and Block 2004). Near the focal point, the trapping force increases linearly with the particle’s displacement from the focal point, following Hooke’s law, $F_{t} = k x(t)$, where $F_{t}$ represents the force applied by the trap, k is the trapping stiffness, and x(t) is the time-varying displacement of the particle from the focal point (see Fig. 2). An accurate measurement of $F_{t}$ relies on an accurate measurement of the position x(t). This requires the addition of an instrumentation to precisely detect the position of the bead. Here, we use back focal plane interferometry for positional detection (Gittes and Schmidt 1998; Lang et al. 2002; Farré et al. 2012). In this technique, a separate detection laser, in combination with a position-sensitive detector (PSD), are used to measure the position x(t), see Fig. 2.

A schematic of our setup is presented in Fig. 2. A Nd:YAG ($\lambda$ = 1064 nm) laser is focused by a water immersion objective (NA=1.20, 60$\times$) mounted on an inverted bright field microscope (Nikon Eclipse Ti-U). Spherical polystyrene beads of radii $a=$ 0.5–2.5 $\upmu$m are used for trapping. Back focal plane interferometry is performed with a detection laser ($\lambda$ = 880 nm). The detection laser is aligned to the trapping laser and focused by the same objective. The laser beam goes through the trapped bead, is focused again through the condenser on the back focal plane, and is detected by a position-sensitive detector (PSD, First Sensor DL100-7). The experimentally acquired electrical signal is converted back into bead position following the same methodology as (Lang et al. 2002). In parallel, videos of the microscopy are synchronously recorded by an sCMOS camera (LaVision PCO.edge) at 100–1000 Hz.

2.2 Generating an unsteady flow: piezo stage calibration

The flow chamber is fabricated by sticking a cover-slip of 175 $\upmu$m thickness to a glass slide with double-sided tape, sealed with silicone grease after filling the sample. The height of the flow chamber is 180 $\upmu$m and the flow is measured in the middle of the flow chamber. For simplicity, we only impose flows along the x-axis. In this study, we characterize oscillatory flows $u(t)= U\sin (2\pi ft)$, where U is the amplitude of the flow and f the frequency. It bears emphasis that OTV is not limited to measuring purely periodic oscillatory flows, but these flows provide the simplest characterization of the operational range of OTV. These uniform flows are generated through the controlled translation of the piezo stage on which the sealed flow chamber is fixed. The bulk motion of the flow chamber results in a uniform flow around the trapped bead. In this study, we generate flows with amplitudes ranging from $U=$ 1.5 to 66 $\upmu$m s$^{-1}$, and frequencies ranging from $f=$ 10 to 90 Hz, corresponding to Péclet numbers ranging from ${\text {Pe}}=$ 0.24 to 63. These flow velocities are measured independently with OTV, and compared with the piezo displacements to characterize the accuracy of the method. This requires a high degree of precision for the calibration of the piezo displacements. This calibration is challenging for two reasons. First, in this range of (U, f), the motion of the piezo stage does not precisely follow the position input, and is non-linear. Therefore, the displacements of the piezo stage cannot be assumed to correspond to the input displacements and require a separate calibration. Second, the oscillatory flows investigated here ($U=$ 1.5–66 $\upmu$m s$^{-1}$ and $f =$ 10–90 Hz) imply very small displacements of the stage. For example, the smallest flow velocity $U=$ 1.5 $\upmu$m s$^{-1}$ corresponds to a maximum displacement of the piezo stage of only $\delta _{{\text {adv}}} \approx$ 50 nm at 10 Hz and is even smaller, $\delta _{{\text {adv}}} \approx$ 5 nm, at 90 Hz. Such small displacements are below the camera resolution of our setup, for which the pixel size is of $\sim$ 100 nm.

To tackle this limitation, we developed a reliable calibration method, which measures the piezo displacement with back focal plane interferometry. We manufactured sample surfaces with polystyrene beads imbedded inside a hard layer of NOA 81. For this purpose, we spin-coated a mixture of liquid polymer NOA 81 and micro-beads on the surface of a glass slide, and cured the mixture under UV light. The calibration surface was then placed on the piezo stage. The hard bonds between the beads and the substrate guarantees that the bead motion precisely follows the piezo stage motion. We then focussed a laser beam through the center of a bead imbedded in the hard NOA 81 layer and activated the piezo stage. The motion of the spherical bead causes a deflection of the laser beam, and back focal plane interferometry was used to deduce the position of the bead as a function of time. The accuracy of this technique is $\sim 1.2$ nm, and this allowed us to calibrated piezo-stage motions as small as 5 nm. This calibration was performed separately for $N = 9$ times for each velocity magnitude U and frequency f to verify the consistency of the piezo stage motion. In this paper, we denote the calibrated velocity of the piezo stage as $u_0(t)$, and refer to it as the ground truth velocity.

Fig. 3

Trap stiffness calibration with the drag force method. Calibration for x- and y-direction are shown in blue and red, respectively. The insets show two examples of the recorded bead positions, x(t), inside the trap. Amplitudes of the bead displacements are obtained by fitting sine waves to the measured x(t)

2.3 Stiffness calibration

Measuring the force exerted on the bead requires the accurate calibration of the trap stiffness k. The stiffness k depends on the refractive index of the medium, the wavelength and power of the trapping laser, and the properties of the trapped bead such as refractive index, size, and shape.

In this study, we employ the drag force method for calibrating the stiffness, as it is the most accurate for measuring hydrodynamic forces exerted on the bead (Capitanio et al. 2002; Jones et al. 2015). In this method, the trapped beads are subjected to known flows u(t), and hence known Stokes hydrodynamic drag forces $F_{h} = 6\pi \mu a u(t)$. We impose uniform flows of amplitudes $U =$ 16–330 $\upmu$m s$^{-1}$ at a low frequency of $f=$ 4 Hz. For these flows, the Reynolds number is negligible ($\sim 10^{-6}$), and the recovery force from the OT, $F_{\mathrm{t}} = k x(t)$, is simply equal to $F_{\mathrm{h}}=6\pi \mu a u(t)$. The amplitude of the bead displacement is extracted by fitting the PSD signal with a sinusoid. The stiffness calibration is demonstrated in Fig. 3, where the measured amplitudes of the bead displacements are plotted against the imposed hydrodynamic drag force, or equivalently the flow velocity. The linear trend (solid lines through the marker) indicates that Hooke’s law for the trap force holds true within this range of flow velocities. The stiffness k of the OT is obtained as the inverse of the slope, see Fig. 3. Within the focal plane, the trap stiffness in the x- and y-direction differ slightly $\sim 20\%$, and they are, therefore, calibrated individually, see Fig. 3.

3 Methodology: optical tweezers-based velocimetry

3.1 Velocimetry measurements: working principle

The flow chamber is first filled with a buffer solution of suspended beads of radii $a=$ 0.5–2.5 $\upmu$m at a low concentration. The flow chamber is subsequently sealed and placed on the piezo stage. A suspended bead is then trapped with the OT in the middle of the flow chamber, 90 $\upmu$m away from the top and bottom walls, and the stiffness of the trap is calibrated for this particular bead, see Sect. 2.3. We generate a uniform flow u(t) around the trapped bead, which we want to measure. This flow is created by activating the oscillatory motion of the piezo stage (U, f). The buffer solution in the sealed chamber follows the motion of the stage, thereby creating a uniform flow $u(t) = U \sin (2 \pi f t)$ around the trapped bead. The bead experiences two external forces: the force due to the optical tweezers $F_{t} = k \varDelta x$, and the hydrodynamic force due to the relative motion between the bead and the external flow $F_{h}(t)$. Hence, the dynamics of the bead is governed by the Boussinesq–Basset–Oseen equations:

$$\begin{aligned} m \ddot{x} = -k x + F_{h}(t), \end{aligned}$$

(1)

where the hydrodynamic force $F_{h}(t)$ includes the Stokes drag, the Basset history force, and the added mass force (Kim and Karrila 1991; Grimm et al. 2012). For a spherical bead in an oscillatory flow of frequency f, such as in our experiments, the expression for $F_{h}(t)$ can be written with complex notations as, $F_{h}(t) = -\gamma \left( 1+\lambda + \frac{1}{9}\lambda ^2 \right) ({\dot{x}} - u(t))$, where $\lambda ^2 = - 2\pi i f a^2/\nu$, i represents the imaginary unit, and $\gamma =6\pi \mu a$. The second term in $F_{h}(t)$ represents the Basset history force and scales with $\Vert \lambda \Vert \sim {\mathcal {O}}$ ($10^{-2}$–$10^{-3}$), while the last term represents the added mass effect and scales with $\Vert \lambda ^2\Vert \sim {\mathcal {O}}$ ($10^{-3}$–$10^{-6}$). Both are small and negligible compared with the first term, which corresponds to the Stokes drag. Furthermore, the polystyrene beads have almost the same density as the surrounding water, and hence the acceleration term on the left side of Eq. (1) has the same order of magnitude as the added mass term and is also negligible in our experiments. This reduces Eq. (1) to a first-order equation for the position x of the trapped bead:

$$\begin{aligned} {\dot{x}} + 2\pi f_c~x = u(t), \end{aligned}$$

(2)

where $f_c$ is the corner frequency of the OT $f_c = k/(2\pi \gamma )$, which characterizes its time of response. Our optical tweezers setup provides time resolved measurements of the bead position x, from which we can deduce the flow velocity u(t) using Eq. (2). Equation (2) is a first-order linear ordinary differential equation, characterized by the time constant $1/f_c$. In the limit, when the timescale of the velocity fluctuations is larger than $1/f_c$, ${\dot{x}}$ in Eq. (2) becomes negligible, and the flow velocity u(t) is directly proportional to x(t). This is the case for a sinusoidal flow u(t) at a frequency lower than $f_c$ for example. In this case, the flow velocity u(t) can be directly deduced from the in-phase sinusoidal position of the bead x(t), such that $u(t) = k/\gamma \, x(t)$, see the bode diagram in Fig. 5. This highlights a fundamental aspect of OTV, namely that measuring the bead position x(t) provides a direct measurement of the flow velocity u(t). This is different for velocimetry methods based on passive tracer particles such as PIV and PTV, for which measurements of particle displacements (i.e. the difference between successive particle positions) are required to infer the velocity. It bears emphasis that, when the timescale of the fluctuations is smaller than $1/f_c$, such as for a sinusoidal flow of frequency larger than $f_c$, the bead position will be phase-delayed and the amplitude of the oscillations in the bead position will decrease, see Fig. 5.

Fig. 4

Three illustrative examples of measuring unsteady flows by OTV, as follow: b $f=$ 10 Hz and $U=$ 10 $\upmu$m s$^{-1}$, d $f=$ 50 Hz and $U=$ 20 $\upmu$m s$^{-1}$, and f $f=$ 90 Hz and $U=$ 36 $\upmu$m s$^{-1}$. The gray curves show the bead displacements, x(t), the blue curves present the flow velocities obtained from the OTV measurements and using the Kalman filter, and the red curves correspond to the ground truth velocities obtained from the piezo stage calibration. In all the examples, the corner frequency of the trap is $f_c=$ 44 Hz. a, c, and e represent the Fourier spectrum of x(t) in (b), (d), (f), respectively. The green dash-lines are Lorentzian fits, which correspond to the Fourier spectrum of x(t) due to a pure Brownian motion of a trapped bead

3.2 Kalman filter

Using Eq. (2), we can in principle deduce the flow velocity u(t) from the measurements of x(t), but this operation involves computing the time derivative ${\dot{x}}(t)$. Our measurements of x(t) from OTV are affected by thermal fluctuations. The noise due to the Brownian motion of the bead can significantly affect the measurements of x(t) when the flow magnitude u is small. To accurately deduce the flow velocity u(t) from our noisy measurements of x(t), we implement a Kalman filter. We consider the discrete-time dynamics of the system at uniform sampling intervals $\varDelta t$. The discrete-time state space model corresponding to Eq. (2) can be written as:

$$\begin{aligned} \begin{aligned} \textit{x}_{n}&=(1-2\pi f_c\varDelta t) \textit{x}_{n-1}+ \varDelta t\ \textit{u}_{n-1} + \textit{w}_{n-1}, \\ \textit{y}_n&= \textit{x}_n + \textit{v}_n, \end{aligned} \end{aligned}$$

(3)

where $x_n$ denotes the discrete-time bead position in the trap at time $t = n \varDelta t$, and $y_n$ denotes the measurement of the bead position obtained from our instrumentation. $\textit{w}_n$ and $\textit{v}_n$ represent, respectively, the process noise and the measurement noise (Pei et al. 2017), and $\textit{u}_n$ represents the flow velocity. It bears emphasis that, in our application of Kalman filtering, we want to estimate the input, i.e. the velocity $\textit{u}_n$. We, therefore, define a model that describes the velocity $u_n$ as a dynamic phasor (Monti et al. 2016; de la Serna and Rodriguez-Maldonado 2011), see Eq. 5 in appendix A. Dynamic phasors, with the form of $a(t)\sin (2\pi f_0 t+\phi (t))$, allow us to approximate any unsteady flow with time varying amplitude, phase, and frequency. We model the input velocity u(t) as a summation of three different dynamic phasors with central frequencies $f_0 = [f-1.2, f+0.25, f+1.2]$, where f is the intrinsic frequency of the flow. We choose on purpose not to include the flow frequency f as one of the dynamic phasors, as f is not known a priori in practical applications. The full detail of the velocity model can be found in appendix A. A forward–backward Kalman filter is implemented to approximate the flow velocities based on state space model in Eq. 3. For more details, please see the appendix B, an open source implementation of the OTV software is available through the 4TU.center for research data (Ghoddoosi Dehnavi and Tam 2020).

3.3 Velocimetry measurements

The working principle of our OTV method is demonstrated by measuring the sinusoidal flow generated by the motion of the piezo stage described in Sect. 2.3. Figure 4 represents three different flow measurements corresponding to flow magnitudes of $U=10,\ 20,\ 36$$\upmu$m s$^{-1}$ with respective frequencies of $f=10,\ 50,\ 90$ Hz. In Fig. 4b, d and f, the raw measurements of x(t) related to each flow are represented in gray. The periodic motion of the bead inside the trap can be clearly seen from the periodicity in the raw signal for x(t). On the left side of Fig. 4, the frequency spectrum of x(t) is shown. Each frequency spectrum of x has a peak that exactly corresponds to the characteristic frequency of the flow (see Fig. 4a, c, e). For frequencies higher than the corner frequency of the trap $f_c$, the magnitude of the spectrum decays with a slope of one. This behaviour is related to the response time of the trap $1/f_c$. Theoretically, the Brownian motion of a trapped bead has a lorentzian frequency spectrum, expressed as ${\mathcal {F}}(x(t)) \varpropto 1 /\sqrt{\left( f_c^2+f^2 \right) }$ (Gittes and Schmidt 1997), where ${\mathcal {F}}$ denotes the Fourier transform (see the green-dashed curves in Fig. 4a–c). The relative values of the characteristic timescale of the flow and the trap ($\tau$ and $1/f_c$) have significant implications in the accuracy of the flow measurement, which will be discussed in the next sections.

The blue lines in Fig. 4b, d, f correspond to our flow measurements using the Kalman filter. Since Eq. 2 is a first-order system, the position of the bead $\frac{k}{\gamma }x(t)$ will be modulated by a phase difference and an amplitude gain compared with the external flow u(t) (see the blue and gray lines in Fig. 4b, d, f). The ground truth unsteady flow velocity (calibrated piezo velocity) is represented in red in Fig. 4b, d, f. Our estimation of the unsteady velocity is in excellent agreement with the actual flow velocity with a maximum error of 4.2% in amplitude of the flow, and a maximum error of 0.7% in the phase of the velocity (see Sect. 4.2 for more details on the accuracy of OTV). Here, we want to emphasize again, that in all these examples, the exact flow frequency is not included as one of the three dynamic phasors in our model for the state representation of the input $u_k$, see Sect. 3.2.

To further validate the performance of our approach and the implementation of the Kalman filter, we extract the phase and the gain from the OTV experiments for a range of magnitudes $U=$ 6–70 $\upmu$m s$^{-1}$ and frequencies $f =$ 10–90 Hz, and compare it to the phase and gain predicted by Eq. 2. The results are presented in Fig. 5. The phase and the gain for each experiment are represented by the red circles in Fig. 5, and are obtained by comparing $u_n$ and $\frac{k}{\gamma }x_n$, computed by the Kalman filter. The phase and gain from the Kalman filter accurately match the theoretical phase and gain from the first-order system in Eq. (2). This indicates that the Kalman filter effectively reconstruct the flow velocities with different intrinsic frequencies from the noisy measurements of x(t).

Fig. 5

Theoretical a gain, and b phase of a first-order system related to Eq. (2). The red circles are the phase and gain deduced from OTV measurements, by comparing $u_n$, the predicted velocity from the Kalman filter, and $2\pi f_c x_n$ the measured displacements from the PSD. The corner frequency of the trap is $f_c=$ 44 Hz

4 Assessment of the accuracy and characteristics of OTV

The performance of our velocimetry method depends on the characteristics of the optical tweezers and in particular on the trap stiffness and the properties of the trapped bead. In this section, we characterize the influence of these OT parameters to formulate guidelines to optimize the performance of the OTV, for any given characteristics of the flow to measure. In addition, we determine the limits and the accuracy of the OTV.

Fig. 6

Signal-to-noise ratio for measuring 50 Hz unsteady flows with two different bead radii: a 2.5 $\upmu$m, and b 0.5 $\upmu$m. The solid lines represent the SNR calculated based on the experimental data. The dash lines represent SNR$^\prime \sim$ 1, which theoretically approximate the SNR. Error bars represent the standard deviation of SNR for different repetitions of the same experiments. All measurements in (a) were performed using the same bead and likewise for (b)

Fig. 7

a, d Examples of measuring an unsteady flow ($U = 6.17$$\upmu$m s$^{-1}$ and $f = 50$Hz) with beads of radius 2.5 and 0.5 $\upmu$m, respectively. The corner frequencies of the trap, corresponding to 2.5 and 0.5 $\upmu$m beads, are 48 Hz and 169 Hz, respectively. We define two measurement errors. First, the error in amplitude, $\varepsilon _{u}$, is deduced from the amplitude of the measured velocity $A_u(t)$, see (b) and (e). Second, the error in phase, $\varepsilon _{\phi }$, is deduced from the unwrapped phase of the measurements $\phi _u(t)$, see (c), (f). We calculate $\phi _u(t)$ and $A_u$ with the Hilbert transform of u(t). The shaded areas in (b, c), and (e, f), represent how much the measured $\phi _u(t)$ and $A_u(t)$ differ from the ground truth

4.1 Influence of the trap stiffness and the bead size on the signal-to-noise ratio

The performance and accuracy of OTV can be characterized by the signal-to-noise ratio (SNR) of our measurements of the bead displacements, x(t), inside of the trap. We estimate the SNR as $\langle A_{x} \rangle /\sigma (v_n)$, where $\langle A_{x} \rangle$ is the average amplitude of $x_n$ obtained from the Kalman filtering, and $\sigma (v_n) = \sigma (y_n - x_n)$ is the standard deviation of noise in our measurements, as explained in the Eq. 3. As discussed in Sect. 3.3, the response time of the trap depends on the corner frequency $f_c = k/(2\pi \gamma )$. To determine the optimal OT parameters, we measured the same oscillatory flow for OTs with different corner frequencies, and report the SNR as a function of $f_c$. The corner frequency is a function of both the stiffness of the trap, k, and the hydrodynamic drag coefficient of the trapped bead, $\gamma$. In our experiments, k was varied by changing the power of the trapping laser, and $\gamma$ was varied using beads of different radii between 0.5 $\upmu$m and 2.5 $\upmu$m.

Figure 6 presents the SNR for three different sinusoidal flows, all at the same frequency $f=$ 50 Hz, and with velocity magnitudes $U =$ 6, 10, 19 $\upmu$m s$^{-1}$. We first consider measurements performed using a bead of 2.5 $\upmu$m radius, see Fig. 6a. As expected, we find that the SNR increases with the magnitude of the external flow measured, regardless of $f_c$. In addition, it is noteworthy that, for all velocity magnitudes, the SNR shows a clear maximum when the corner frequency of the trap is $f_c = f =$ 50 Hz. This indicates that, to maximize the SNR, the trap stiffness k should always be chosen, such that the corner frequency of the trap precisely coincides with the frequency of the sinusoidal flow.

The maximum SNR value can be rationalized by comparing the displacements of the bead due to the external flow with the Brownian displacements of the bead. First, the bead displacements x(t), measured with the PSD, scale with $U/2\pi \sqrt{f^2+f_c^2}$, see Fig. 5a. Second, our measurements are affected by the displacements due to thermal fluctuations, which scale with $\sqrt{k_B T/2 \pi f_c \gamma }$. Therefore, the theoretical prediction for the SNR can be defined as SNR$^\prime = U\sqrt{\frac{\gamma }{2 \pi k_B T}} \sqrt{\frac{f_c}{f^2+f_c^2}}$. SNR$^\prime$ is represented in Fig. 6 (see dashed lines), and agrees with the SNR values that we estimated from experiments. When the corner frequency of the trap is low ($f_c < f$), corresponding to a low trap stiffness, the bead displacements are significantly affected by thermal fluctuations. In addition, the response time of the trap is then high, leading to a low amplitude of the bead displacements, see Fig. 5a. Both effects lead to a low SNR for $f_c<f$. For $f_c>f$, on the other hand, the OT is stiffer, which decreases the noise amplitude due to thermal fluctuations. However, the amplitude of the signal (bead displacements due to the external flow) decreases as well, leading to an overall decrease in SNR for large $f_c$. In between these two limits, the SNR reaches a theoretical maximum when the corner frequency of the trap coincides with the frequency of the external flow.

Next, to investigate the effect of the bead size, we repeat the same experiments with a 0.5 $\upmu$m bead radius and different laser powers. Using smaller beads is advantageous when a high spatial resolution is required, e.g., measuring the flow very close to a boundary such as a solid wall. The use of a smaller bead, implies that the corner frequency, $f_c = k/2\pi \gamma$, is much larger. In this case, the corner frequency cannot coincide with the frequency of the flow ($f \sim$ 50 Hz), because it requires lowering the trap stiffness to very low values, for which it is no longer possible to trap the bead inside the OT. The minimum corner frequency which we are able to reach with a 0.5 $\upmu$m bead is around 169 Hz. It is important to consider that at a given laser power, SNR$^\prime$ for 0.5 $\upmu$m beads are smaller than for 2.5 $\mu$m beads, see Fig. 6b.

In Fig. 7a and d, we illustrate an example of OTV measurement of the same flow ($U =$ 6.1 $\upmu$m s$^{-1}$ and $f =$ 50 Hz), with beads of 0.5 and 2.5 $\upmu$m radii. The corner frequency of the trap for a 0.5 and a 2.5 $\upmu$m bead radius are $f_c =$ 169, 50 Hz, respectively. The corresponding SNR$^\prime$ numbers for 0.5 and 2.5 $\upmu$m beads are 0.29 and 0.86, respectively. Unlike the measurements done with the 2.5 $\upmu$m bead, the 0.5 $\upmu$m bead measurements fail to accurately capture the flow. It is evident that the low SNR of the measurements with 0.5 $\upmu$m beads leads to poor accuracy of the measurement in Fig. 7d. This example illustrates how the choice of the bead size can affect the accuracy of the OTV. Smaller beads can be beneficial to measure flows with large characteristic frequencies f; for instance, beads of 0.5 $\upmu$m radius can be used to measure flows with $f>169$ Hz optimally.

Fig. 8

OTV accuracy in measuring unsteady flows with various amplitudes and frequencies: a the error in the amplitude, $\varepsilon _u$, and b the error in the phase, $\varepsilon _{\phi }$. The red curve represents the diffusion limit of a passive tracer particle, and the white curve corresponds to SNR$^\prime \sim$ 1 for a trapped bead

4.2 Accuracy of the measurements

The accuracy of the velocimetry method is quantified by comparing the experimental OTV measurements u(t) with the known flow velocities from the calibration experiments $u_0(t)$. Unsteady, oscillatory flows are characterized by the amplitude of the velocity and the phase of the oscillations. Hence, we define two errors to quantify the accuracy of velocity measurements: the error in amplitude $\varepsilon _u$ and the error in phase $\varepsilon _{\phi }$. To do this, we compute the time varying amplitude $A_u(t)$ and the phase $\phi _u(t)$ of the measured velocities u(t) by computing the Hilbert transform of H(u(t)) (Feldman 2011). The analytic signal of u(t), defined by ${\hat{u}}(t) = u(t)+i~H(u(t))$, is a complex-valued time series, which has a modulus and a phase. The modulus of ${\hat{u}}(t)$ corresponds to the envelope of the amplitude of u(t), and represents the flow velocity magnitude at each time $A_u(t)$.

Fig. 9

Comparing the accuracy of OTV for measuring unsteady flows, $f =$ 50 Hz and $U=$ 1.5–70 $\upmu$m s$^{-1}$, using beads of 0.5 and 2.5 $\upmu$m radii: a the error in the amplitude, $\varepsilon _u$, and b the error in the phase, $\varepsilon _{\phi }$. The blue and red dashed lines correspond to SNR$^\prime \sim$ 1 for 0.5 $\upmu$m and 2.5 $\upmu$m bead radii, respectively

The argument of ${\hat{u}}(t)$ represents the phase of the oscillatory flow $\phi _u(t)$. $A_u(t)$ and the unwrapped phase $\phi _u(t)$ are illustrated, respectively, in Fig. 7b, c, e and f. Figure 7b, c, e and f represents respectively the $A_u(t)$ and $\phi _u(t)$ of the flow measurement of the sample flow $U =$ 6.17 $\upmu$m s$^{-1}$ and $f =$ 50 Hz performed with the 2.5 $\upmu$m and 0.5 $\upmu$m bead. We define the measurement errors, $\varepsilon _u$ for the amplitude, and $\varepsilon _{\phi }$ for the phase, by comparing the amplitude and phase of the measured velocity ($A_u(t)$, $\phi _u(t)$) with those of the calibrated ground truth (${A_u}_0(t)$, ${\phi _u}_0(t)$) such that:

$$\begin{aligned} \begin{aligned} \varepsilon _{\phi } = 1/2\pi \sqrt{\frac{1}{\tau }\int _0^\tau \left( \phi _u(t) - {\phi _u}_0(t) \right) ^2 {\text {d}}t}, \\ \varepsilon _{u}=1/U \sqrt{\frac{1}{\tau } \int _0^\tau \left( A_u(t) - {A_u}_0(t) \right) ^2 {\text {d}}t}. \end{aligned} \end{aligned}$$

(4)

Fig. 10

An example of measuring the same unsteady flow ($U =$ 6.17 $\upmu$m s$^{-1}$ and $f =$ 50 Hz) with the OTV and the micro-PTV methods. A bead of 2.5 $\upmu$m radius is used in both of the methods. a The Fourier spectrum of the micro-PTV’s tracer particle displacements, x(t). b The gray line shows the displacements of the free tracer particle, and the blue line is the predicted velocity using the Kalman filter. c The Fourier spectrum of the bead displacements, x(t), for a bead trapped in the optical tweezers. d The gray line is the scaled bead displacements, $\frac{k}{\gamma }x(t)$, and the blue line is the velocity predicted by the Kalman filter. The corner frequency of the trap is $f_c =$ 48 Hz

Figure 7b, c, e, and f illustrates $A_u(t)$, $\phi _u(t)$, ${A_u}_0(t)$, and ${\phi _u}_0(t)$ for a particular measurement. We estimate $\varepsilon _{\phi }$ and $\varepsilon _{u}$ for $U =$ 1.5–70 $\upmu$m s$^{-1}$ and $f = 10-90$ Hz, and present results in Fig. 8. For each flow condition, we repeat the experiment nine times, each recording for 3 s. For all of the measurements in Fig. 8, we used a 2.5 $\upmu$m bead and a fixed trap stiffness of 13.9 pN $\upmu$m$^{-1}$ corresponding to $f_c =$ 44 Hz. The insets in Fig. 8 are close-up views of the error maps for $U<$ 25 $\upmu$m s$^{-1}$. The red curve in Fig. 8 represents the diffusion limit, also discussed in Fig. 1. Both errors in phase and amplitude of velocities are very small far from the diffusion limit, when ${\text {Pe}} \ge 1$, and remain small below the diffusion limit, when ${\text {Pe}} \le 1$, see Fig. 8. At very low flow amplitudes and higher frequencies, the errors increase. This corresponds to measurements in the regime, when the SNR of the OT signal decreases and SNR $<1$. SNR $\sim 1$ corresponds to the diffusion limit of the trapped particle. This limit is reached when $U\sim \sqrt{\frac{2 \pi k_B T}{\gamma }} \sqrt{\frac{f^2+f_c^2}{f_c}}$, and is represented by the white curve in Fig. 8. Below this limit, neither the phase nor the amplitude of the measured velocities can be measured accurately.

We use the errors, $\varepsilon _{\phi }$ and $\varepsilon _u$, to compare the accuracy of OTV with a 2.5 $\upmu$m bead radius to that with a 0.5 $\upmu$m bead radius. We measure unsteady flows at a fixed frequency of $f =$ 50 Hz and decrease the velocity magnitudes from $U=$ 70 $\upmu$m s$^{-1}$ to $U =$ 1.5 $\upmu$m s$^{-1}$, using beads of both sizes. For each bead size, we use the laser power which gives us a corner frequency leading to the highest SNR: $f_c =$ 169 Hz for 0.5 $\upmu$m beads and $f_c =$ 50 Hz for 2.5 $\upmu$m beads. For both bead sizes, $\varepsilon _{\phi }$ and $\varepsilon _{u}$ increase with decreasing amplitude U of the measured velocity, see Fig. 9a, b.

The theoretical limit SNR$^\prime \sim 1$ corresponds to $U =$ 26.8 $\upmu$m s$^{-1}$ and $U =$ 7.1 $\upmu$m s$^{-1}$, for bead radii $a = 0.5~\upmu$m and 2.5 $\upmu$m, respectively (see the red and blue dash lines in Fig. 6). For $U \gtrsim$ 26.8 $\upmu$m s$^{-1}$, the theoretical signal-to-noise ratio is such that SNR$^\prime > 1$ for both bead sizes. For these flows, the errors are small for both bead sizes, with the amplitude error $\varepsilon _{u}<$14%, 7.5%, and the phase error $\varepsilon _{\phi }<$2%, 1.5%, respectively, for 0.5 and 2.5 $\upmu$m beads. For $U <$ 26.8 $\upmu$m s$^{-1}$, the accuracy of the OTV with a 0.5 $\upmu$m bead substantially decreases, especially for $\varepsilon _{\phi }$, see Fig. 6. For the 2.5 $\upmu$m bead, however, when 7.1 $<U<$ 26.8 $\upmu$m s$^{-1}$, the OTV measurements remain accurate, $\varepsilon _u<$20% and $\varepsilon _{\phi }<$4%. For $U<$ 7.1 $\upmu$m s$^{-1}$, the signal-to-noise ratio is such that SNR$^\prime <1$ for the 2.5 $\upmu$m bead as well, and the accuracies of measurements start to decrease. As Fig. 9 suggests, measurements with the 2.5 $\upmu$m bead show much better accuracies than the ones with the 0.5 $\upmu$m bead. For $U<$ 3 $\upmu$m s$^{-1}$, the errors in amplitude and phase are large for both cases, which means the measured velocity is not representative of the actual flow. For 50 Hz flows, OTV with 2.5 $\upmu$m beads allows us to measure flows as small as 3–5 $\upmu$m s$^{-1}$, while with 0.5 $\upmu$m beads, OTV is limited to flows larger than 20 $\upmu$m s$^{-1}$.

4.3 Comparison between OTV and micro-PTV

Fig. 11

Comparing the velocimetry errors of the OTV and the micro-PTV methods with beads of 2.5 $\upmu$m radius: a the error in the amplitude, $\varepsilon _u$, and b and the error in the phase, $\varepsilon _{\phi }$. Error bars represent the standard deviation of $\varepsilon _{\phi }$ and $\varepsilon _u$ for different repetitions of the same experiments. The corner frequency of the trap is $f_c =$ 48 Hz

Fig. 12

A case study of measuring the unsteady flow velocity by OTV close to the cilia of C. reinhardtii. a A 2.5 $\upmu$m trapped bead, 35 $\upmu$m away from the C. reinhardtii, is used to measure the flow velocity. The corner frequency of the trap is $f_c =$ 32.3 Hz. b and d present the Fourier spectra of the bead displacements in x and y-direction, respectively. c and d are comparison of the measured velocity by OTV experiments with BEM simulations. $u_x(t)$ and $u_y(t)$ are the velocity of the flow in the x and y-direction, respectively

In this section, we compare the accuracy of OTV and micro-PTV measurements. The micro-PTV measurements are preformed by first trapping a bead in the flow cell and releasing it just before the onset of piezo stage motion. We do so to keep the bead in the focal plane of the camera and in the middle of the flow cell. The bead motion is recorded at a frame rate of 990.85 fps by the camera. The bead position is extracted from the recordings, using the TrackMate extension in Fiji (Tinevez et al. 2017). A Laplacian of Gaussian segmentation with quadratic sub-pixel localization method is used in TrackMate to detect the bead position. The micro-PTV data is filtered following the same approach as the one described for the OTV in Sect. 3.2. We use the same Kalman filter, described in Eq. (3), for a trap stiffness of $k=0$. The OTV measurements are performed with a 2.5 $\mu$m bead inside a trap of stiffness $k=15$ pN $\upmu$m$^{-1}$, corresponding to $f_c =$ 48 Hz. Figure 10b and d present measurements of the same unsteady flow, $U =$ 6.1 $\upmu$m s$^{-1}$ and $f=$ 50 Hz using, respectively, micro-PTV and OTV. Brownian motions clearly dominate the micro-PTV measurements, and despite the Kalman filter, the velocity cannot be measured, see Fig. 10b. This can also be seen from the spectrum of the bead displacements in Fig. 10a. In comparison, OTV measurements of the same flow are far less noisy, and the Kalman filter can reconstruct the actual flow velocity accurately, see Fig. 10c, d. The example of Fig. 10 corresponds to ${\text {Pe}}\sim$ 0.35 < 1, for which methods based on passive tracers will be inaccurate, see Sect. 1. For OTV, on the other hand, the SNR$^\prime \sim$ 1, and the measurement errors remain low, see Fig. 8.

We quantify the differences in accuracy with the phase and amplitude errors. Figure 11 compares $\varepsilon _u$ and $\varepsilon _{\phi }$ for micro-PTV and OTV. Both $\varepsilon _u$ and $\varepsilon _{\phi }$ are lower for OTV than for micro-PTV. The flow amplitudes corresponding to ${\text {Pe}} \sim$ 1 and SNR$^\prime \sim$ 1 are respectively $U= 17.9$ and 7.1 $\upmu$m s$^{-1}$, and are indicated with dashed lines in Fig. 11. Both methods are appropriate for high-amplitude flows, however, for lower amplitudes, when ${\text {Pe}}<1$, micro-PTV is no longer accurate, see the increase of $\varepsilon _u$ and $\varepsilon _{\phi }$ in Fig. 11. For OTV, on the other hand, the errors remain low even for ${\text {Pe}}<1$, and slowly increase for very low amplitudes, when SNR$^\prime <1$, see Fig. 11. For unsteady flows at 50 Hz, the use of micro-PTV is limited to flow amplitudes larger than $\sim$ 20 $\upmu$m s$^{-1}$, while OTV can measure 50 Hz flows with amplitudes as low as 3–5 $\upmu$m s$^{-1}$.

5 Case study: unsteady flow around C. reinhardtii

Finally, we demonstrate the efficiency and robustness of OTV by measuring unsteady biological flows, in the $f =$10–100 Hz frequency range. We use OTV to measure the 2D flow field around the flagella of the micro-organism C. reinhardtii. The cilia, or flagella, of this micro-organism are approximately 10–12 $\upmu$m long, and beat at about 40–70 Hz. We capture one cell with a glass pipette at 120 $\upmu$m above the bottom surface of a flow cell, to limit the influence of wall effects on our velocity measurements. It bears emphasis that the beating of the flagella is not strictly periodic and the beating frequency fluctuates over the course of time. In addition, the beating frequency is not known a priori. We measure the velocity with OTV, using a bead trapped 35 $\upmu$m away from the center of the cell body. We use a bead with radius $a=$ 2.7 $\upmu$m and the trapping laser has a stiffness of $k=$ 9.85 pN $\upmu$m$^{-1}$ and a corner frequency of $f_c=$ 32.3 Hz, see Fig. 12a. Here, we define the bead displacement in x- and y-direction as x(t) and y(t), with the corresponding velocities $u_x(t)$ and $u_y(t)$. Figure 12b, d represents the spectrum of x(t) and y(t). In this case study, the main frequency peak is broad, see Fig. 12b, compared to the flows generated with the piezo stage, see Fig. 4c, and the signal consists of several harmonics, here at 50, 100 and 150 Hz (Fig. 12b). Each harmonic in the x(t) or y(t) signal correspond to a different harmonic in $u_x(t)$ or $u_y(t)$ with a different phase and gain. Each harmonic of $u_x(t)$ or $u_y(t)$ will affect the bead displacements x(t) and y(t) differently according to the first-order dynamics of the bead described in Eq. (2), and characterized by the Bode diagrams on Fig. 5. Figure 12c, e represents the measured bead displacements and the reconstruction of the velocity $u_x(t)$ and $u_y(t)$ from the Kalman filter. The presence of several harmonics in the x(t) and y(t) signals leads to reconstructions of the input velocities $u_x(t)$ and $u_y(t)$, which significantly differ from x(t) and y(t) in amplitude, shape and phase. To verify the accuracy of the velocimetry method, we compare the OTV measurements with the numerical boundary element method (BEM) simulations (Wei et al. 2019). We model the cell body and the pipette as a single axisymmetric object. Briefly, our BEM simulations use a combination of stresslet singularities on the solid surfaces, and stokeslet and rotlet singularities along the cell-pipette centreline to impose the non-slip boundary conditions on the solid boundaries. Furthermore, we use slender body theory to represent the flagella. The shapes and deformations of the flagella are directly extracted from the experiment recordings and used in BEM simulations to compute the flow fields (see Wei et al. (2019) for more details). In Fig. 12, the predicted velocities by the BEM simulations confirm the accuracy of OTV to measure the amplitude and phase of the unsteady flow generated by the C. reinhardtii. Here, we want to emphasize that although the flow frequency in this example is fluctuating with time, and is not known a priori, the Kalman filter can accurately reconstruct the unsteady flow velocity.

6 Conclusion

In this paper, we established a velocimetry method, called OTV, to accurately measure micro-scale unsteady flows. We implement a Kalman filter model to predict instantaneous flow velocities from the optical tweezers measurements. We characterized the accuracy of the OTV measurements for different bead sizes and corner frequencies of the OT. The velocimetry method is the most accurate for the highest value of the SNR. This is reached when the trap’s corner frequency coincides with the characteristic frequency of the unsteady flow. Therefore, for an optimal accuracy of the OTV measurements, it is recommended to always tune the corner frequency of the OT to the characteristic frequency or time scale of the flow. The corner frequency depends on several parameters including the bead radius and the trap stiffness, which also depends on the bead size, its index of refraction, the design of the trap and the laser power. In general, larger beads correspond to smaller corner frequencies and will measure more accurately flows with lower characteristic frequencies. In our experiments, 2.5 $\upmu$m beads could be used to optimally measure unsteady flows with characteristic frequencies lower than $\sim$300 Hz, while 0.5 $\upmu$m beads are advantageous to measure flows with larger characteristic frequencies of up to 3000 Hz. We quantify the OTV accuracy to measure unsteady flows with amplitudes 1.5–70 $\upmu$m s$^{-1}$ and frequencies 10–90 Hz. Our results demonstrate that the OTV method is capable of measuring unsteady flows below the diffusion limit of a free tracer particle, ${\text {Pe}}<1$, which are not measurable with micro-PIV or micro-PTV methods. This higher accuracy lies in the fact that a trapped particle can offer a higher SNR than a free particle. For micro-PTV, the error in the phase is very large, when measuring flows below the diffusion limit ${\text {Pe}} < 1$. On the other hand, the OTV method can accurately measure flows below the diffusion limits. For 50 Hz flows below the diffusion limit, which are typical to the ciliary flows, the maximum OTV errors in the amplitude and phase are $\sim$16.1% and $\sim$4.3%, respectively. Due to this high temporal accuracy, OTV can be effectively used in studies related to measuring unsteady flows in microfluidic and biofluidic studies, especially studies related to microorganisms.

Appendices

Appendix

A velocity model

We model the flow velocity, u(t), as a summation of different harmonics with the central frequencies $f_{0}^j$, as in (de la Serna and Rodriguez-Maldonado 2011):

$$\begin{aligned} \begin{aligned} u(t)&= \sum _{j} a^j(t)\cos (2\pi f_{0}^j t+\phi ^j(t)) \\&= \mathfrak {R}\left( \sum _{j} p^j(t)e^{i2\pi f_{0}^j t}\right) =\sum _{j}\mathfrak {R}({\mathcal {U}}^j(t)), \end{aligned} \end{aligned}$$

(5)

where $a^j(t)$ and $\phi ^j(t)$ are the time-varying amplitude and phase of u(t) corresponding to the $j_{th}$ central frequency, and $\mathfrak {R}$ represents the real part. $p^j(t) = a^j(t) e^{i\phi ^j(t)}$ is referred to as the dynamic phasor and ${\mathcal {U}}^j(t)$ is the rotating phasor. Equation (5) implicitly suggests that the instantaneous frequency related to jth harmonic in Eq. (5) can dynamically change as $f^j(t) = f_0^j + {\dot{\phi }}^j(t)/2\pi$ (Monti et al. 2016). The state transition of the dynamic phasor $p^j(t)$ at time $t=n \varDelta t$ can be estimated by a K order Taylor expansion (de la Serna and Rodriguez-Maldonado 2011). In this study we use a Taylor expansion up to the second order, such that:

$$\begin{aligned} \varvec{{\mathcal {P}}}^j_{n} = \begin{pmatrix} 1 &{} \varDelta t &{}\frac{\varDelta t^2}{2!}\\ 0 &{} 1 &{} \varDelta t \\ 0 &{} 0 &{} 1 \end{pmatrix} \varvec{{\mathcal {P}}}^j_{n-1} = \tilde{\varvec{\varPhi }} \varvec{{\mathcal {P}}}^j_{n-1}, \end{aligned}$$

(6)

where $\varvec{{\mathcal {P}}}^j_n=[p^j_n,\, {\dot{p}}^j_n,\, \ddot{p}^j_n]^T$ is the state vector, ${\dot{\square }}$ and $\ddot{\square }$ stand for the first- and second-order time derivatives, and $p^j_n$ is the truncated approximation of $p^j(t)$ at time $t = n \varDelta t$. Therefore, the velocity $u_n$ can be approximated as:

$$\begin{aligned} \begin{aligned} u_n&= \sum _{j} \mathfrak {R}\left( {{\varvec{h}}}^T \varvec{{\mathcal {P}}}^j_n e^{i2\pi f_0^j t}\right) \\&=\frac{1}{2}\sum _{j} \begin{pmatrix} {{\varvec{h}}} \\ {{\varvec{h}}} \end{pmatrix}^T \begin{pmatrix} \varvec{{\mathcal {U}}}^j_n \\ \bar{\varvec{{\mathcal {U}}}}^j_n \end{pmatrix}, \end{aligned} \end{aligned}$$

(7)

where ${{\varvec{h}}} = [1,\, 0,\, 0]^T$, $\varvec{{\mathcal {U}}}^j_n=\varvec{{\mathcal {P}}}^j_n e^{i2\pi f_0^j t}$, and the bar represents the complex conjugate. The state transition of $\varvec{{\mathcal {U}}}^j_n$ can be derived using Eq. (6) (de la Serna and Rodriguez-Maldonado 2011):

$$\begin{aligned} \begin{pmatrix} \varvec{{\mathcal {U}}}^j_{n} \\ \bar{\varvec{{\mathcal {U}}}}^j_{n} \end{pmatrix} = \begin{pmatrix} \psi ^j \tilde{\varvec{\varPhi }} &{} {{\varvec{0}}}\\ {{\varvec{0}}} &{} {\bar{\psi }}^j \tilde{\varvec{\varPhi }} \end{pmatrix} \begin{pmatrix} \varvec{{\mathcal {U}}}^j_{n-1} \\ \bar{\varvec{{\mathcal {U}}}}^j_{n-1} \end{pmatrix}, \end{aligned}$$

(8)

where $\psi ^j = e^{i2\pi f_0^j \varDelta t}$, and $\tilde{\varvec{\varPhi }}$ is defined in Eq. (6). Equations (7)–(8) provide us with a model to approximate the velocities of the unsteady flows.

B Kalman filter

We rearrange the state space model from Eq. (3) to an extended version in order to include Eqs. (7)–(8) in our model. The extended state vector is defined as $\varvec{{\mathcal {X}}}_n = [x_n,\, \varvec{{\mathcal {U}}}^{1T}_n,\,\bar{\varvec{{\mathcal {U}}}}^{1T}_n,\ldots ,\, \varvec{{\mathcal {U}}}^{jT}_n,\,\bar{\varvec{{\mathcal {U}}}}^{jT}_n]^T$, and therefore, equation 3 becomes:

$$\begin{aligned} \begin{aligned} \varvec{{\mathcal {X}}}_n&= \mathbf{A } \varvec{{\mathcal {X}}}_{n-1} +{{\varvec{w}}}_{n-1},\\ y_n&= \mathbf{H } \varvec{{\mathcal {X}}}_n + v_n,~~~~~~~~~~ 1\le n\le N, \end{aligned} \end{aligned}$$

(9)

where $\mathbf{H } = [1,\,0,\ldots ,\,0]$, $y_n$ is the measurements of the bead position, $\varvec{w}_n$ is the process noise with the covariance matrix $\mathbf{Q }$, and $v_n$ is the measurement noise with the covariance r. The transition matrix $\mathbf{A }$ is given by:

$$\begin{aligned} \begin{aligned} \mathbf{A } =\left( \begin{array}{c|ccc} 1-\frac{k}{\gamma } \varDelta t&{}\frac{\varDelta t}{2}{{\varvec{h}}}_2^T&{}\cdots &{}\frac{\varDelta t}{2}{{\varvec{h}}}_2^T\\ \hline \\ 0 &{} \varvec{\varPhi }_1 &{} &{} \\ \vdots &{} &{} \ddots &{} \\ 0 &{} &{} &{} \varvec{\varPhi }_j\\ \end{array} \right) , \end{aligned} \end{aligned}$$

(10)

where ${{\varvec{h}}}_2=[{{\varvec{h}}}^T,\, {{\varvec{h}}}^T]^T$, and $\varvec{\varPhi }_j$ is the transition matrix in Eq. (8). The Kalman filter utilizes the state space model in Eq. (9) to estimate the state vector $\varvec{{\mathcal {X}}}_{n}$ from the noisy measurements $y_n$. In this study we implement a forward–backward Kalman filter rather than a single Kalman filter. Using a single Kalman filter can cause unwanted time lags in estimating the flow velocity u(t). A forward-backward Kalman filter not only increases the accuracy of filtering, but also prevents any possible time lag in estimating the velocity u(t). The forward part of the filter recursively calculates the best estimate of $\varvec{{\mathcal {X}}}_{n}$ based on the measurements before time instance n, as follow:

$$\begin{aligned} \varvec{{\mathcal {X}}}_{n|n-1}= & {} \mathbf{A } \varvec{{\mathcal {X}}}_{n-1|n-1}, \end{aligned}$$

(11)

$$\begin{aligned} \varvec{\varSigma }_{n|n-1}= & {} \mathbf{A }\varvec{\varSigma }_{n-1|n-1}\mathbf{A }^T+\mathbf{Q }, \end{aligned}$$

(12)

$$\begin{aligned} \mathbf{K }_{n}= & {} \varvec{\varSigma }_{n|n-1}\mathbf{H }^T(\mathbf{H }\varvec{\varSigma }_{n|n-1}\mathbf{H }^T+r)^{-1}, \end{aligned}$$

(13)

$$\begin{aligned} \varvec{{\mathcal {X}}}_{n|n}= & {} \varvec{{\mathcal {X}}}_{n|n-1} +\mathbf{K }_{n}(y_n-\mathbf{H }\varvec{{\mathcal {X}}}_{n|n-1}), \end{aligned}$$

(14)

$$\begin{aligned} \varvec{\varSigma }_{n|n}= & {} (\mathbf{I }-\mathbf{K }_{n}\mathbf{H })\varvec{\varSigma }_{n|n-1}. \end{aligned}$$

(15)

The forward filtering, which is presented in Eqs. (11)–(15), includes two main steps. Equations (11) and (12) represent the first step, which gives an a priori estimation of the state vector $\varvec{{\mathcal {X}}}_{n|n-1}$, and state-error covariance matrix $\varvec{\varSigma }_{n|n-1}$. This estimation is purely based on the state space model in Eq. (9). In the second step, the measurements $y_n$ are combined with the state vector estimation by the Kalman gain $\mathbf{K }_{n}$. Therefore, the posterior estimation of the state vector $\varvec{{\mathcal {X}}}_{n|n}$ and state-error covariance matrix $\varvec{\varSigma }_{n|n}$ can be calculated through Eqs. (13)– (15).

The backward filter is a Rauch $\acute{\mathrm{T}}$ung $\acute{\mathrm{S}}$triebel smoother (Maybeck 1982) that estimates the state vector, conditional on all the measurements $y_n$, as follow:

$$\begin{aligned} \mathbf{C }_{n}= & {} \varvec{\varSigma }_{n|n} \mathbf{A }^T \varvec{\varSigma }_{n+1|n}^{-1}, \end{aligned}$$

(16)

$$\begin{aligned} \varvec{{\mathcal {X}}}_{n|N}^s= & {} \varvec{{\mathcal {X}}}_{n|n} + \mathbf{C }_{n}(\varvec{{\mathcal {X}}}_{n+1|N}^s-\varvec{{\mathcal {X}}}_{n+1|n}), \end{aligned}$$

(17)

$$\begin{aligned} \varvec{\varSigma }_{n|N}^s= & {} \varvec{\varSigma }_{n|n}+ \mathbf{C }_{n}(\varvec{\varSigma }_{n+1|N}^s-\varvec{\varSigma }_{n+1|n}) \mathbf{C }_{n}^T, \end{aligned}$$

(18)

where $\mathbf{C }_{n}$ is defined as the smoothing estimator gain matrix. $\varvec{\varSigma }_{n+1|n}$ and $\varvec{{\mathcal {X}}}_{n+1|n}$ in Eqs. (17) and (18) can be calculated similar to Eqs. (11) and (12). Equations (16)–(18) estimate the smoothed state vector $\varvec{{\mathcal {X}}}_{n|N}^s$ and the state-error covariance matrix $\varvec{\varSigma }_{n|N}^s$ backward in time, from time instance $n=N$ to $n=1$, with the boundary conditions $\varvec{{\mathcal {X}}}_{N|N}^s = \varvec{{\mathcal {X}}}_{N|N}$ and $\varvec{\varSigma }_{N|N}^s = \varvec{\varSigma }_{N|N}$.