Optical tweezers (OTs) are renowned for their non-invasiveness in the trapping and manipulation of microscale and nanoscale objects, which find wide applications in biology and nanotechnology [1–3]. In the past few decades, the manipulation of polystyrene (PS) spheres [4], cells [5], DNA molecules [6], and proteins [7] has been realized and several microstructures for optofluidic manipulation have been demonstrated [8–12], by using single-beam (diffraction limited) or near-field optical trapping techniques. Yang et al. proposed a dielectric slot waveguide capable of condensing accessible electromagnetic energy to capture and transport DNA molecules [13]. Cai and Poon utilized a silicon nitride micro-disk resonator with whispering-gallery modes to manipulate PS microparticles [14]. Yoon et al. demonstrated a double-well potential landscape and directly observed the hopping of a single 4-nm particle by using metallic nanoantennas with various gap sizes and different optical pump powers [15]. In addition, several applications for trapping and manipulating particles have been proposed and demonstrated, such as switching [16], sorting [17], storing [18], and conveying [19].

In order to estimate the trapping ability of an optical tweezer system, trapping precision and trapping stability quantification have been proposed and improved continuously. Traditionally, users often adapt an empirical approach for optimizing their experimental parameters. In the previous studies, a stability factor (S), which is proportional to the depth of the potential well (S=ΔU/kBT), is assigned to empirically describe the trapping stiffness [20]. Particles can be considered to be stably trapped when S≫1. Ultimately, S only provides a binary outcome since it only takes care of the energy aspects of the system, which is related to the depth of a potential well. In addition, the particle motion and distribution will be influenced by not only the depth, but also the width and shape of the potential well. Therefore, S alone is incapable of assessing the trapping stability and a new approach is required to evaluate the characteristics of trapping. Then, many researches on trapping or escaping time based on Kramers time have been developed [11,21,22]. Simon and Libchaber have studied the residence time distribution in a double potential well when a 1-μm Brownian particle is trapped [23]. Šiler and Zemanek focused on the transition problem of particles in periodic potential wells and proposed the mean optical trap escape time (MOTET) to quantify the particle motion [24]. Ferrer et al. have researched the influence of fluid viscoelasticity on transitions between a double well [25]. However, the researches have rarely considered and experimentally validated the influence of external disturbances on trapping time.

In the actual applications, the external disturbance cannot be overlooked. As shown in the conceptual image of Fig. 1, Brownian motion [26,27], relative moving of environmental liquid, and the laser-beam-induced thermal convection [28] can affect particle motion and distributions. All the above disturbance forces will diversely deform the potential well in different directions [29] and the escape rate will be increased on the tilting side (lower barrier side of the potential well). As the directional flow velocity continuously increases, the depth of trapping the potential well would reduce to the level of average kinetic energy of the particle. This case can be considered as a weak trapping system resulting in the dynamic equilibrium process of trapping and escaping. Moreover, the residence time of particles in the potential well will satisfy a statistical distribution, which is a relatively strict description for particle trapping in a potential well. Therefore, an accurate prediction of average trapping time or trapping time distribution can be calculated [30,31]. Clearly, based on current researches on static potential wells, there is an undisputed need to develop a comprehensive model for quantifying the trapping stability with an external flow.

All of the above quantifying description is based on obtaining a potential well, but for Gaussian beams, the potential well can be directly derived. In practical applications, a potential well can only be calculated by integrating the optical gradient force. And the force distribution only can be obtained from numerical simulation or experimental measurement, which is not very straightforward. Therefore, we take a Gaussian beam as an example because it is the most commonly used light beam for optical tweezers. The force distribution of a Gaussian beam can be well defined when a particle can be viewed as a dipole and can be uniquely determined by its power and waist width. Therefore, its trapping stability could be quantitatively identified.

In this work, we have studied trapping stability when Brownian motion, relative movements, and external flow are considered. The kinetics of an optically trapped object have been formulated firstly in the static symmetric potential well by the Fokker–Planck (FP) equation. Based on the revised equation, the statistics of particle position in a flow-skewed potential well, as well as a trapping time that provides a simple but accurate descriptor of the trapping stability, have been calculated. To verify the model, we have set up an optical tweezer system with a back-focal-plane interferometer that can monitor the Brownian motion trajectory of the trapped particles. The experimentally verified model can serve as a look-up function to determine trapping stability for a given set of operation conditions, implying strong potential for many branches of physical and biological applications.

We first consider the case of an ideal static symmetric potential well. When only gradient force (Fop), viscous resistance (Fdrag), and Brownian motion (FBr) are considered, the generated potential well will be symmetric, which is termed as a static potential well. Fop can be obtained based on the method of Maxwell’s stress tensor (MST) when the electric field distribution is defined. In this case, the position and motion trajectory of a single spherical dielectric particle, which will be influenced by the randomness of Brownian motion, can be numerically simulated with a given time interval using a Langevin equation [32,33]. Note that the Langevin equation used in the simulation is an overdamped Langevin equation, which is valid for most of the liquid background environment. When the environment material becomes vacuum or gas, inertia and an underdamped Langevin equation should be considered [34]. In addition, the effect of hydrodynamic memory has also been neglected [35]. Furthermore, the distribution of particles can be simulated when multiple particles are calculated simultaneously and individually (termed as simulation method in the following). Moreover, for a large number of particles, in the perspective of statistics, the distributions of Brownian particles and potential well should satisfy the relationship based on thermodynamic theory and governed by the Fokker–Planck (F-P) equation as [36] dP(x,t)dt=ddx[U′(x)γP(x,t)]+kBTγd2P(x,t)dx2,(1)where kB is the Boltzmann constant, T is ambient temperature, γ is the Einstein-Stokes coefficient, U(x) represents the potential well, and P(x,t) is the probability density distribution of particles. x and t represent spatial coordinate and time. Obviously, P(x,t) in Eq. (1) can describe both the temporal and spatial evolution characteristics of particle distributions (termed as theoretical method in the following). When particle distributions tend to equilibrate, the diffusion current also tends to vanish and steady distribution will be established. In this case, supposing a Gaussian-shaped potential well and escape model schemed in Fig. 2(a), a 3D potential well and steady distributions by both simulation and theoretical methods have been shown in Fig. 2(b). Most of the particles are concentrated at the bottom of the potential well. Therefore, for a steady distribution obtained from the F-P equation, the population of trapped particles is positively correlated to the depth of potential well U(x).

According to Eq. (1), besides the temporal and spatial evolution characteristics, diffusion current can also be obtained and further an escape rate can be calculated under evolutionary transients [24]. Intuitively, we have also proposed the mean first trapping time (MFTT) τ, defined as the mean duration that a particle can stay in a certain potential well, reciprocal of escape rate. By the simulation method, MFTT of a single particle can be calculated as the time moving from the depth to the barrier of the potential well. When the theoretical method is considered, taking the potential distribution into the F-P equation, we can have the probability distribution of trapped particles, the escape rate, as well as the MFTT of particles: τ=12·γkBT|∫−∞+∞e−S(x)dx∫baeS(x)dx|,(2)where γ represents the coefficient of viscous resistance. a and b represent the barrier and bottom of the potential well. After a Taylor expansion, MFTT can be simplified by ignoring higher-order terms above the third order [24]: τ≅14·2πγ|U′′(a)U′′(b)|eS,(3)where S is the stability factor mentioned in Section 1 and represents the depth of the potential well. Obviously, An Arrhenius relation exists between the trapping time and the depth of the potential well.

When a Gaussian-shaped potential well is supposed, the second derivative of U(x) can be described in the equation precisely. Notably, in order to simplify the calculation, ±3σ has been defined as the barrier of the potential well because the trapping time will show no evident difference when the barrier extends beyond 3σ. σ is the standard deviation of the Gaussian potential well. Then, the MFTT for spherical dielectric particles in a Gaussian symmetric potential well can be expressed as τ=τ0eSeS,(4)where τ0 represents the trapping time when the well depth is 1kBT and is called the characteristic value of trapping time for a unit potential well. τ0 satisfies the equation τ0=CDησ2, where C can be represented by exp[3.25+ln3π2/(42kBT)], and D is particle diameter. Therefore, the escape rate of particles will decrease and trapping time becomes longer with the increase of D. η represents the liquid viscosity, which is related to temperature. In addition, η used in our simulation has neglected Faxén’s correction because the diameter of a capillary is much larger than particle size [37]. σ is the standard deviation of the Gaussian potential well (representing the width). It is worth noting that the actual parameter that affects the trapping time is the curvature at the bottom and edge of the potential well as in Eq. (3). However, in Gaussian potential wells, the curvature is directly related to the width, so we will take width as one of the important influencing factors in our later discussions. According to Eq. (4), trapping time in a potential well will be related to not only the depth, but also the width of the potential well. Therefore, traditional descriptions such as S on trapping stability will be not accurate enough because of the neglect of width. In addition, MFTT is also proportional to the size of trapped particles because greater optical gradient force will be exerted on particles with larger diameters.

Moreover, to reveal the exponential relationship between τ and S more intuitively, the logarithm of MFTT lnτ=lnτ0+S−lnS−1 is introduced. Therefore, the exponential relationship can be transformed into a quasi-linear relationship. For the trapping of particles with diameters of 200, 300 nm, and in a symmetric Gaussian potential well with standard deviations of 25, 100, 150 nm, we can directly calculate their MFTT shown as the straight lines in Fig. 2(c) by Eq. (2), which coincides well with the numerical simulation method represented as the scattered points. It also shows that there is a good linear relationship between lnτ and Gaussian well depth when the depth is much larger than 1kBT. In addition, we have also calculated MFTT with Eq. (4) shown as the dashed lines in Fig. 2(c) and found that the results of trapping time will be decreased a little because of the simplification of equations during the derivations. Here, we have discussed the cases when S≤10 because with the increase of S, MFTT will be much larger and the trapping tends to be stable. It is worth noting that the intersections on the y axis give the characteristic value τ0, and these curves for different traps are parallelly transformed. An animated visualization of particle motion in a static potential well is also presented in a supporting video (Visualization 1).

For the specific case when a Gaussian beam is incident and the particle diameter satisfies D≪λ/10, the electric field distribution is definite on the basis of Gaussian parameters because particles can be approximated as dipoles, which will not influence the field. Therefore, the width and depth of the potential well under a Gaussian beam can be calculated directly. We apply the form of a Gaussian beam in the formula of gradient force F=14α∇|E|2, where α represents electric polarizability. The field intensity of Gaussian beam E can be formulated by E=E0exp(−r2ω02), where ω0 represents its waist radius and r represents the radial coordinate for the Gaussian beam. Then, its corresponding potential well can be calculated by integrating F along the propagation direction as U=α2|E0|2exp(−2r2ω02). We can find that there is a linear relationship between the waist of the beam and the waist of the potential well, i.e., σ=k1ω0, and k1 should be 0.5. Since the waist of the potential well σ can be preliminarily predicted by parameters of the Gaussian beam, the trapping time for the unit potential well τ0 can be represented by CgDηω02, where the constant Cg=k12exp[3.25+ln3π2/(42kBT)]. In addition, we can find that the depth of the potential well can also be represented by the beam power P and beam waist ω0 as S=k2Pω02, where k2=α2πkBT.

As an example, assuming a Gaussian beam with waist radius ω0 of 1 μm and power of 200 mW, we compared the result of potential wells obtained by simulation and theoretical methods as shown in Fig. 2(d). Notably, the width for comparison here is full width at half-maximum (FWHM) because of simple markings in the figure. For the simulation case, we have simulated the force exerted on a 50-nm particle. Then, potential wells have been calculated by integrating the force along propagation direction and fitted to Gaussian curves shown as the purple dashed line in Fig. 2(d). FWHM and S are calculated as 1.14 μm and 1.2 (σ can be fitted as 0.48 μm) so that the trapping time can be derived as 317.6 ms. For the theoretical case, σ is calculated as 0.5 μm according to the waist of the Gaussian beam; S can be calculated as 1.0, which is close to the simulation calculations. In addition, the trapping time can be derived as 350.7 ms. The two cases match very well with similar depths, widths, and further trapping times, which shows the availability of our theoretical model.

Now we consider a more realistic scenario with external disturbance. In an actual optical tweezer system, uniform external flow, capillary force, and laser-beam-induced thermal convection are unavoidable because the focused laser in a capillary will generate a photothermal effect and further turbulence in different directions will be produced. When the photothermal effect caused by the laser-beam-induced thermal convection is considered, both the particle dispersion and the trapping stability should be corrected [28]. Here we apply this correction using COMSOL Multiphysics simulations. Figure 3(a) shows the relative motion of the aqueous medium due to thermal convection generated by the focused trapping beam. A self-excited fluid flow field (with velocity up to 12 μm/s) is presented along both the vertical and lateral directions around the focal spot so that potential wells with different inclinations in orthogonal directions will be generated. In such case of laser-beam-induced thermal convection, we have neglected the change of temperature-dependent viscosity due to the photothermal effect. In addition, even for the case of laser-beam-induced dynamic flow with non-uniform distribution as shown in Fig. 3(a), the mean trapping time in different locations can also be calculated, since the flow distribution around the tiny particle could be assumed uniform.

Considering a dynamic potential well, a more general scenario, all the external effects or perturbation will be reflected in the modification of the potential well. The relative uniform motions between the fluid and the beam (with either a static beam in a flow, or a moving/delivering beam in a static medium) can be regarded as a component of the viscous resistance force [38] calculated by 6πRηu shown as Fdrag in Fig. 3(b), where R represents the radius of particles. u denotes the velocity of uniform external flow. The force curves will move down so that trapping potential well integrated from optical force will tilt on one side shown in Fig. 3(c). Assuming the external flow is uniform, the so-called effective dynamic potential well satisfies Ueff(x)=U(x)+γux, where x is the displacement of the particle relative to the potential well, and U(x) represents the originally static symmetric potential well. When relative velocity and Brownian motion are considered, we have simulated and visualized the 2D and 3D potential wells and particle distribution with the relative velocity of 200 μm/s shown in Figs. 3(d) and 3(e). The presence of relative-motion-induced resistance force will tilt the originally balanced potential well through modifying the barrier heights in opposite manners, i.e., of unbalanced barriers ΔU1 and ΔU2 at either side as shown in Fig. 3(d). With the effective potential well distribution, as shown in Fig. 3(e), we can also directly calculate the particle distribution by using the F-P equation, similar to that in the symmetric case. In addition, we notice that the lower barrier side always has a higher particle density, hence a higher chance of escape from the well. Since the potential well becomes asymmetric, the escape rate is unbalanced. The net escape rate is the sum of the rates on either side of the barrier, and the MFTT could be expressed as τ=2τ1τ2/(τ1+τ2), where τ1 and τ2 can be calculated using Eq. (3) with ΔU1 and ΔU2. In a weak trapping system, the depth of the potential well on the tilting side will be much smaller than the other, i.e., ΔU1≪ΔU2, so we can approximately consider ΔUeff=ΔU1 and τ2→∞, while MFTT can be simplified as τ=2τ1. We can define the effective stability factor Seff=ΔUeff/kBT. Therefore, according to Eq. (4), trapping time in a dynamic potential well can be calculated directly by τ=2τ0eSeffeSeff.(5)

We demonstrate the effective depth and MFTT in the dynamic potential well with a range of flow velocity, as shown in Fig. 3(f). With the increasing relative velocity and particle size, the viscous resistance force will be enhanced, and the effective depth shown as dotted lines will be reduced, corresponding to the left axis. On the right side, the MFTT exponentially decreases with respect to an increase of relative velocity. The theoretical results (solid curves) coincide well with simulated results (scatter points). An animated visualization of particle motion in a flow-skewed potential well with different flow velocities is presented in a supporting video (Visualization 2). Interestingly, when the relative velocity increases to a critical value (uc) to tilt the well to an extent, the trapping time will decrease dramatically, which is corresponding to the condition that the barrier on this side disappears and the particle will immediately escape. The critical velocities uc for 500 nm, 300 nm, and 200 nm particles are 239, 399, and 598 μm/s, respectively. More specifically, the critical velocity is determined by a balance between the gradient trapping force, which can be influenced by the size of particles and power of the Gaussian beam, and the external directional viscous resistance force.

A Gaussian beam is also considered in such a dynamic condition. Moreover, the parameters of the dynamic potential well σeff (effective standard deviation) and Seff (effective stability factor) cannot be calculated by Gaussian parameters directly because the potential well has deformed due to directional flow velocity. Therefore, we have defined a and b as the solution of Feff=0, which represent the inflection point and the depth of the dynamic potential well as shown in Figs. 3(b) and 3(c) so that S1 can be represented as ΔUeff/kBT, where ΔUeff=Ueff(a)−Ueff(b), which is formulable based on Gaussian beams.

The width of the potential well will also influence the characteristics of trapping. The potential well is integrated by the force exerted on the particle, so a wider potential well will correspond to a smaller force when the depth of potential wells is consistent. Therefore, a very small Fdrag caused by external flow will make the dynamic potential well reach the critical velocity. In a static potential well, dispersion time of particles will become longer with the wider static potential well. Consequently, the trapping characteristics will be different when the width of the potential well is variant even with the same depth as shown in Eqs. (4) and (5). Then, we have analyzed a three-dimensional ellipsoidal potential well formed by a perpendicular focus point with different widths in the major axis (longitudinal z direction) and the minor axis (transverse x direction). Its trapping stability and critical velocity will differ. Even though the depth of wells is the same for both directions, the potential well along the major axis has a larger width with a smaller gradient force, which is more likely to be affected by flow-velocity-induced viscous force. It is worth noting that scattering force Fsc is also an important factor that can influence the trapping characteristics in the longitudinal z direction. When only Fsc is considered without external flow, we can calculate the equilibrium position or MFTT by adding an additional item Fsc in the Langevin equation. In the theoretical model, Fsc can be viewed as a contributing factor of external disturbances, which will generate external perturbance and tilt the potential well. Here, we mainly want to study the influence of uniform external flow on potential wells with different widths. Figures 4(a) and 4(c) show the biasing of the trapping force and the tilting of the potential well, respectively, due to the flow velocity in the transverse direction, while Figs. 4(b) and 4(d) are for the case of flow in longitudinal direction. Since the trapping force has a lower peak in longitudinal direction, a slight biasing of the viscous force with 64 μm/s flow will make the force peak become nearly zero as shown in Fig. 4(b). As a result, the potential barrier on the left side becomes nearly zero as shown in Fig. 4(d). Definitely, the particle will escape from the longitudinal side. The dynamic motions of particles in this kind of ellipsoidal potential wells are shown in Visualization 3.

Further, we have also explored the relationship of MFTT calculated by the simulation method and Eq. (5) on the basis of Fig. 4. As shown in Table 1, we can find that MFTT in a dynamic potential well calculated by different methods keeps consistency. MFTT calculated by the theoretical method is based on Eq. (4), so the trapping time is a little smaller than that by the simulation method as shown in Fig. 2(c). Furthermore, even a small flow velocity in longitudinal direction generated by laser-induced thermal convection, as small as 12 μm/s in Fig. 3(a), is sufficient to introduce a significant tilting of the potential well shown in Fig. 4, and dramatical reduction of trapping time. In addition, when velocity is up to 64 μm/s, the effective depth of the dynamic potential well in longitudinal direction is near zero and the trapping time is tiny. However, even with the same depth of the static potential well, in transverse direction, ΔUeff can still maintain 5.2kBT and MFTT will be more than 100 ms. Therefore, MFTT will also be seriously affected by the width of the potential well. We have also shown errors between the two methods. Errors may be increased when trapping time is very short because a very tiny difference will lead to a large ratio error. Interestingly, when a dynamic potential well is applied, the error becomes smaller, because the error induced by a one-side infinite potential well assumption partially compensates for the error of the theoretical model shown in Fig. 2(c). Overall, when the trapping characteristic value of an optical tweezer system is determined, one can accurately predict particle capture events and trapping time for a given set of system parameters. It largely extends the applications of the stability analysis based on trapping time.Table 1.

Comparison of MFTT Calculated by Simulation and Theoretical Method

In the above analysis, all the calculations and theories are used to obtain mean trapping time. On the other hand, the experimental values are randomly distributed. Given the energy-dependent nature of the system, it should satisfy a certain statistical distribution [39]. Since the escape process is a homogeneous Poisson process and each escape event is a low probability random event, from the perspective of statistics the probability density distribution of trapping times p(t) can be expressed as a Poisson distribution p(t)=C1e−C2t [23,40]. Since the sum of probability density distribution is one and its expectation is the average value of trapping time, it can be found that the coefficient C1 is equal to C2, which characterizes the escape rate of particles. Therefore, the probability density distribution of trapping times satisfies p(t)=1τ¯(u)e−tτ¯(u),(6)where τ¯(u) represents MFTT when the velocity of uniform external flow is u. Obviously, the intercept of the probability density distribution of trapping time on the y axis represents the escape rate, that is, the reciprocal of MFTT, which is a function of flow velocity in a dynamic potential well. With the increasing flow velocity, MFTT will be reduced and the theoretical distribution of trapping time will decrease faster.

In order to study the characteristics of trapping time distribution experimentally, we have set up an optical tweezers system operating simultaneously in a laminar-flow-supporting capillary, as shown in Fig. 5(a). The system is designed to capture the flowing particles with different relative velocities, particle sizes, and beam power levels. A 976-nm laser is tightly focused by a high numerical aperture (NA) objective to form a focus spot, and a 40× objective is used for collecting light travelling in the forward-scattering direction. A back-focal-plane interferometer (BFPI) with a quadrant photodiode device (QPD) is used for monitoring the trap by recording the Brownian motion trajectory of particles [41]. The trapping duration is measured with millisecond precision.

Figure 5(b) shows a Brownian motion trajectory (one of the dimensions) recorded by the BFPI, with clear partitions of an initially empty trap (signal stays at the center), particle entering the trap with presence of a laminar flow (signal jumps quickly from the upstream side to the downstream side of the trapping beam), random walking at the bottom position of the skewed potential well for a period (signal fluctuation at the bottom position), and final escaping (signal returns to the beam center). Based on this technique (an experimental video of the signals is presented in Visualization 4 and Visualization 5), we are able to observe a weak trapping process and measure the trapping duration. Figure 5(c) shows the actual experimental system.

Shown in Fig. 5(d) is the QPD electrical signal of trapping a single 500-nm PS with 57-μm/s flow rate and 35-mW optical power. Two specified enlarged figures of particle trapping signal are also shown in Figs. 5(e) and 5(f). According to the figures, the electrical signal on the QPD shows a great change in the transmitted light signal when a single particle is trapped. The sudden downward step is that the particle just enters the potential well and then is trapped stably (at this time, the QPD signal is displayed as rising). During the stable trapping time, the signal remains at a height and floats with the Brownian motion, and then is released after stable trapping for a period of time in the flow rate environment.

The experimental results and theoretical distribution of trapping time for particles with a diameter of 300 nm in the aforementioned Gaussian potential well are shown in Fig. 6(a). To compare the theoretical and experimental results more clearly, the distribution curves of trapping time are converted to logarithmic scale and shown in Fig. 6(b). The slope of the curves in Fig. 6(b) corresponds to the intercept of the curves shown in Fig. 6(a). For three relative velocities 57, 113, and 227 μm/s, MFTTs τ calculated by theory are 9.3 ms, 6.9 ms, and 5.4 ms, respectively, and MFTTs τexpr measured in our experiment are 10.6±0.2ms, 7.2±0.11ms, and 5.3±0.07ms. When the relative velocity is 113 μm/s and the beam power is 15 mW, MFTTs obtained by theory, simulation, and experiment of 300 nm PS are 6.9 ms, 6.8 ms, and 7.2 ms. Furthermore, in order to analyze the trapping time more comprehensively, the experiments with different particle diameters (200, 300, 500 nm) and trapping beam powers (15, 25, 35 mW) have also been carried out. As seen from Figs. 6(c) and 6(d), MFTT is proportional to particle size and trapping beam power. In conclusion, the experimental results coincide well with the theoretical predictions of trapping time. Our results demonstrated that the distribution of trapping time satisfies Poisson distribution from the perspective of statistics, and it is significantly affected by relative velocity, particle size, and beam power.

In conclusion, we have proposed a comprehensive method from a statistical point of view and quantified the trapping stability by trapping time and particle distribution. Specifically, we have implemented the MFTT for describing particle dispersion and transient trapping process in a potential well, especially for weak trapping systems. For a practical optical tweezers system excited by Gaussian beams, we have demonstrated that when the physical parameters (beam power, spot size, and liquid properties) are determined, the width and depth of potential wells can be calculated directly when the particle diameter is less than λ/10. Besides the analysis of the above static symmetric condition, we have also presented the dynamic potential well model when relative motion is introduced by regarding the effect of the relative velocity as a directional viscous resistance force. Then, MFTT can be calculated directly by considering the parameter of flow velocity. The differences between the method of theory and simulation have also been analyzed. Furthermore, we have set up an optical tweezers system in a microfluidic channel to capture the flowing particles and measure the trapping time at different relative velocities, particle sizes, and beam powers to verify our model and the Poisson distributions of MFTT. In total, our model quantitatively reveals the effects of fluidic flow velocity, particle size, and optical power on trapping time, which will find applications where optical trapping stability is challenged by external perturbations in weak trapping conditions.

The sample used in the trapping time measurement experiment is polystyrene spheres. The trapping time of polystyrene spheres with diameters of 200, 300, 500 nm in a water environment was measured.

The velocity field generated by thermal convection is simulated by using a finite-element mode solver (COMSOL Multiphysics) together with optimized mesh size and perfectly matched layer. The parameters of the microfluidics channel and the source are consistent with the experiment.

The electromagnetic field is solved by using a finite difference time domain solver (Lumerical FDTD). In the simulation, we simulated the Gaussian beam and set the microfluidics channel at its waist position. On the other hand, the force of PS particles is simulated and calculated by the Maxwell stress tensor (MST) method in FDTD software. The parameters of source and PS particles are consistent with the experiment.

The motion of the potential well and particles is solved by MATLAB. Through the iterative algorithm based on the Euler method, we solved the motion equations and calculated the trapping time of several groups of PS particles, each group of 10,000 particles, and the time interval is one-tenth of the characteristic time of particles because it can strike a balance between the accuracy and the efficiency of the simulation. MFTT is obtained by the average value of each group, and the distribution of trapping time is represented by the histogram of simulated results.