%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1. Introduction
Microswimmers in flow exhibit a wide range of fascinating behaviours, such as directional alignment (Hope et al. Reference Hope, Croze, Poon, Bees and Haw2016), hydrodynamic trapping (Sipos et al. Reference Sipos, Nagy, Leonardo and Galajda2015) and aggregation (Torney & Neufeld Reference Torney and Neufeld2007), all of which are heavily influenced by their surrounding flow environment. The hydrodynamic interactions between microswimmers’ activity and the flow field play a crucial role in shaping these dynamics (Lauga Reference Lauga2016; Mathijssen et al. Reference Mathijssen, Figueroa-Morales, Junot, Clément, Lindner and Zöttl2019; Baker et al. Reference Baker2019). Understanding these hydrodynamic interactions is critical, as they not only govern microswimmer behaviour but also have important applications in areas like targeted drug delivery (Park et al. Reference Park, Jiang, Oncay and Metin2017), bioengineering (Liu et al. Reference Liu2022) and environmental monitoring (Liang et al. Reference Liang, Zhan, Zeng, Xu, Wang, Zhao, Zhang, Guo, Feng and Ma2018). As a result, investigating the physical mechanisms behind microswimmer–flow interactions, and leveraging this knowledge for microswimmer manipulation, has become a key focus in the study of active matter (Bechinger et al. Reference Bechinger, Leonardo, Löwen, Reichhardt, Volpe and Volpe2016; Wheeler et al. Reference Wheeler, Secchi, Rusconi and Stocker2019; Aranson Reference Aranson2022).
Early research on swimming behaviours of microswimmers primarily focused on planar shear flow environments (Zöttl & Stark Reference Zöttl and Stark2012; Kantsler et al. Reference Kantsler, Dunkel, Blayney and Goldstein2014; Lee et al. Reference Lee, Lohrmann, Szuttor, Auradou and Holm2021; Rubio et al. Reference Rubio, Potomkin, Baker, Sen, Berlyand and Aranson2021). When bacteria move in the shear flow, they tend to accumulate in regions of high shear rate (Rusconi, Guasto & Stocker Reference Rusconi, Guasto and Stocker2014), which also is theoretically demonstrated by Vennamneni, Nambiar & Subramanian (Reference Vennamneni, Nambiar and Subramanian2020). They find that the interaction between the shear rate and bacteria is the crucial factor for achieving migratory motion and aggregation behaviour. More recently, the study of microswimmers has been extended to their swimming dynamics in vortex flows (Berman & Mitchell Reference Berman and Mitchell2020). For instance, Marcos & Stocker (Reference Stocker2006) demonstrate that bacteria align with streamlines in strong vortices, a result consistent with our findings for passive particles. For gyrotactic micro-organisms, such as algae, studies by Durham et al. (Reference Durham, Climent and Stocker2011, 2013) have highlighted how gyrotaxis influences their aggregation dynamics in vortical flows. Furthermore, Sokolov & Aranson (Reference Sokolov and Aranson2016) find that elongated bacteria in the vortex flow generated by a rotating microparticle cross the streaminline and form a depletion zone near the vortex core. Tanasijevic & Lauga (Reference Tanasijevic and Lauga2022) have employed a theoretical model to demonstrate that the depletion zone is caused by the interaction between particles and the vortex flow. However, in the Gaussian vortex flow field, the trapped active particles aggregate around the vortex core to form a high-density region (Arguedas-Leiva & Wilczek Reference Arguedas-Leiva and Wilczek2020). Additionally, dumbbell-shaped particles exhibit spiral motion in a two-dimensional steady vortex flow, with the shape of their orbits depending on their initial position and direction (Yerasi, Govindarajan & Vincenzi Reference Yerasi, Govindarajan and Vincenzi2022).
Optical image of the acoustic streaming flow generated by the oscillating microbubble in experiment, reproduced from Ahmed et al. (Reference Ahmed, Ozcelik, Bojanala, Nama, Upadhyay, Chen, Hanna-Rose and Huang2016) under terms of the CC-BY 4.0. licence. Copyright 2016 The Authors, published by Springer Nature.
Given the interaction between the particles and background flow field, the use of microfluidic chips has attracted widespread interest in the biomedical field (Collins et al. Reference Collins, Khoo, Ma, Winkler, Weser, Schmidt, Han and Ai2017; Lu et al. Reference Lu, Martin, Soto, Angsantikul, Li, Chen, Liang, Hu, Zhang and Wang2019; Geng et al. Reference Geng2023; Kutluk, Viefhues & Constantinou Reference Kutluk, Viefhues and Constantinou2024; Liu et al. Reference Liu, Zhang, Zhang, Dong, Jiang, Huang, Li, Xie and Qiu2024). The acoustic streaming flow generated by an ultrasonically driven and trapped microbubble is a primary method used in microfluidic chips for manipulating particles (Marmottant & Hilgenfeldt Reference Marmottant and Hilgenfeldt2004; Chen et al. Reference Chen, Fang, Merritt, Strack, Xu and Lee2016; Zhang et al. Reference Zhang, Bachman, Ozcelik and Huang2020). As shown in figure 1, the trapped cavitated microbubble undergoes radial and tangential oscillations under ultrasonic excitation, producing a stable acoustic streaming flow comprising a pair of vortices (Ahmed et al. Reference Ahmed, Ozcelik, Bojanala, Nama, Upadhyay, Chen, Hanna-Rose and Huang2016). Rogers & Neild (Reference Rogers and Neild2011) find that the behaviour of passive particles within the acoustic streaming flow is influenced by their size and density, making this method widely applicable in the separation of cells (Li et al. Reference Li, Liu, Huang, Ohta and Arai2021; Gao et al. Reference Gao, Wu, Luan, Papautsky and Xu2022). Additionally, experimental evidence demonstrates that Escherichia coli can be collected near the vortex centre of the acoustic streaming flow generated by an oscillating bubble, leading to the formation of a biofilm structure (Yazdi & Ardekani Reference Yazdi and Ardekani2012). Further experiment has demonstrated that multiple oscillating microbubbles of varying sizes can create an acoustic network flow, which facilitates the non-destructive pumping of Escherichia coli (Gao et al. Reference Gao, Wu, Lin, Zhao and Xu2020). Despite these advancements, the effects of acoustic streaming flow on the swimming dynamics and aggregation behaviours of active particles, as well as the underlying physical mechanisms governing these processes, remain poorly understood and require further elucidation.
Inspired by the observed aggregation behaviour of bacteria in an acoustic streaming field during experiments (Yazdi & Ardekani Reference Yazdi and Ardekani2012), we have conducted a detailed investigation into the swimming dynamics of active particles in the acoustic streaming flow generated by an oscillating microbubble. We employ an approximate solution for the acoustic streaming flow (Spelman & Lauga Reference Spelman and Lauga2017), and then we establish a deterministic model to predict the trapped bounded orbits of both passive and active particles. Our results shows that active particles are trapped in the bounded orbits near the vortex centre due to their activity, while passive particles are restricted to moving along the streamlines. The stable motion region of active particle depends on the streaming flow, which further proves the dependence of active micro-organisms on the environment. To further understand the behaviour of particles in a more realistic environment, we introduce rotational noise to analyse the trapping and distribution of both passive and active particles. We find that active particles exhibit a non-uniform distribution within the streaming flow and we identify a flow configuration that maximises the degree of aggregation of the active particles. Our results provide insights for the use of microfluidic chips in the rapid detection and separation of micro-organisms.
2. Set-up and acoustic streaming flows
2.1. Set-up and non-dimensionalisation
The microswimmer is modelled as a spheroidal particle characterised by an aspect ratio
$\lambda$
, defined as the ratio of the long axis to the short axis. We assume that the swimmer maintains a constant intrinsic swimming speed
$v_0\boldsymbol {p}$
, where
$\boldsymbol {p}$
denotes the swimmer’s orientation. For an elongated swimmer, the orientation aligns with the direction of the long axis.
The streaming flow is generated by an acoustically driven oscillating microbubble with radius
$a$
. For simplicity, we assume that the bubble is fixed in an unbounded domain and ignore wall effects. In addition, we assume that the motion of the bubble’s surface is axisymmetric, incorporating both radial and tangential oscillations. The surface of the microbubble is presumed to oscillate at an angular frequency
$\omega$
with a small amplitude
$\epsilon a$
, where
$\epsilon \ll 1$
. We employ spherical coordinates
$(r,\varphi )$
as illustrated in figure 2. Using a Lagrangian framework, the boundary shape of the microbubble is described by its radial position
$R$
and angular position
$\Phi$
, measured at the location of a material point initially at
$(a,\varphi )$
. Due to the axisymmetry, the surface oscillation modes can be expressed as a series of orthogonal basis functions formed by Legendre polynomials, presented as follows:
\begin{align} R(\varphi ,t)&=a-\epsilon a\sum ^{\infty }_{n=0}V_nP_n(\cos \varphi )e^{i(\omega t+\frac {\pi }{2})},\\ \Phi (\varphi ,t)&=\varphi +\epsilon \sum ^{\infty }_{n=1}W_n\left (\frac {\int _{\cos \varphi }^1 P_n(x){\rm d}x}{\sin \varphi }\right )e^{i(\omega t+\frac {\pi }{2})},\nonumber \end{align}
(
$a$
) Illustration of a microswimmer moving in the acoustic streaming flow generated by an oscillating bubble. The microswimmer is depicted as an orange ellipsoid particle located at (
$r$
,
$\varphi$
) with constant velocity
$v_0$
along
$\boldsymbol {p}$
. The magnitude of scaled streaming flow is
$|\psi |/(\epsilon ^2\delta \omega a^3)$
at
$V_0=0.4$
,
$V_1=1$
. (
$b$
),(
$c$
) Time evolution results of the real part of boundary position (b)
$R$
and (c)
$\Phi$
of the oscillating bubble for one period based on the same parameters with representative value
$\epsilon =0.01$
.
where
$V_n$
and
$W_n$
are arbitrary complex constants determined by the surface motion, and
$P_n(x)$
is the Legendre polynomial of degree
$n$
. The flow field is governed by the Navier–Stokes equation:
\begin{align} &\frac {\partial \boldsymbol {u}}{\partial t}+(\boldsymbol {u}\cdot \boldsymbol {\nabla })\boldsymbol {u}=-\frac {1}{\rho }\boldsymbol {\nabla } p+\nu \boldsymbol {\nabla }^2\boldsymbol {u}, \end{align}
where
$\boldsymbol {u}=(u_r, u_{\varphi })$
represents the fluid velocity and
$p$
is the pressure. Since the flow generated by the oscillating bubble is axisymmetric, we introduce a streamfunction
$\psi (r,\varphi )$
. In spherical coordinates, the fluid velocities are represented as follows:
\begin{align} &u_r=\frac {1}{r^2\sin \varphi }\frac {\partial \psi }{\partial \varphi },\quad u_{\varphi }=-\frac {1}{r\sin \varphi }\frac {\partial \psi }{\partial r}. \end{align}
Substituting (2.3) into (2.2), the Navier–Stokes equation becomes
\begin{align} &\frac {\partial \left (D^2\psi \right )}{\partial t}+\frac {1}{r^2}\!\left [\!-\frac {1}{\sin \varphi }\frac {\partial \psi }{\partial r}\frac {\partial \left (D^2\psi \right )}{\partial \varphi }+\frac {1}{\sin \varphi }\frac {\partial \psi }{\partial \varphi }\frac {\partial \left (D^2\psi \right )}{\partial r}+2D^2\psi L\psi \!\right ]=\nu D^2\left (D^2\psi \!\right )\!, \end{align}
where the operators are
\begin{align} &D^2=\frac {\partial ^2}{\partial r^2}+\frac {1}{r^2}\left (\frac {\partial ^2}{\partial {\varphi }^2}-\frac {\cos \varphi }{\sin \varphi }\frac {\partial }{\partial \varphi }\right ),\quad L=\frac {\cos \varphi }{\sin ^2\varphi }\frac {\partial }{\partial r}-\frac {1}{r\sin \varphi }\frac {\partial }{\partial \varphi }\cdot \end{align}
On the surface of the bubble, the fluid velocity satisfies the boundary condition
\begin{align} &u_{r}(R, \Phi )=\frac {\partial R}{\partial t}\left |_{(R, \Phi )}, \quad u_{\varphi }(R, \Phi )=R\frac {\partial \Phi }{\partial t}\right |_{(R, \Phi )}. \end{align}
In addition, the incorporation of an additional boundary condition is crucial to guarantee the absence of any tangential stress exerted on the surface of an oscillating bubble:
\begin{align} &\left. \left [\frac {1}{r}\frac {\partial u_r}{\partial \varphi }+\frac {\partial u_{\varphi }}{\partial r}-\frac {u_{\varphi }}{r}\right ]\right |_{(R, \Phi )}=0. \end{align}
We will sketch the asymptotic solution of (2.4) in § 2.2. We assume that the disturbance of the flow field due to the active particle is negligible. Additionally, the hydrodynamic interactions between the active particle and the bubble surface are neglected in this study, as their effects are significantly weaker compared with the dominant influence of the acoustic streaming flow (Spagnolie et al. Reference Spagnolie, Moreno-Flores, Bartolo and Lauga2015). Consequently, the motion of the particle, as it drifts with the flow, and the rotation of the particle can be determined by Jeffery’s equation (Jeffery Reference Jeffery1922). The position of the particle
$\boldsymbol {r}$
and the orientation
$\boldsymbol {p}$
are governed by the following equations:
\begin{align} \frac {{\rm d} \boldsymbol {r}}{\partial t}&= v_0\boldsymbol {p}+\boldsymbol {u},\nonumber\\ \frac {{\rm d} \boldsymbol {p}}{\partial t}&=\frac {1}{2}(\boldsymbol {\nabla } \times \boldsymbol {u})\times \boldsymbol {p}+\gamma (\boldsymbol {E}\cdot \boldsymbol {p}-(\boldsymbol {p}\cdot \boldsymbol {E}\cdot \boldsymbol {p})\boldsymbol {p}), \end{align}
where
$\boldsymbol {E}= (\boldsymbol {\nabla }\boldsymbol {u}+\boldsymbol {\nabla }\boldsymbol {u^T} )/2$
is the strain rate tensor of the flow field,
$\gamma = (\lambda ^2-1 )/ (\lambda ^2+1 )$
represents the shape parameter of particle,
$\lambda$
represents the ratio of the long axis to the short axis of microswimmer, and
$(-1\lt \gamma \lt 0)$
,
$(\gamma =0)$
and
$(0\lt \gamma \lt 1)$
represent oblate spheroid, sphere and prolate ellipsoid, respectively.
Here, we first focus on obtaining solutions of the acoustic streaming flow field. We choose
$\omega ^{-1}$
as the relevant time scale and the radius
$a$
of the oscillating bubble as characteristic length to non-dimensionalise (2.4). Consequently, the dimension of the streamfunction is
$\omega a^3$
. The dimensionless form of (2.4) is then as follows:
\begin{align} &\frac {\partial \left (D^2\psi \right )}{\partial t}+\frac {1}{r^2}\!\left [\!-\frac {1}{\sin \varphi }\frac {\partial \psi }{\partial r}\frac {\partial \left (D^2\psi \right )}{\partial \varphi }+\frac {1}{\sin \varphi }\frac {\partial \psi }{\partial \varphi }\frac {\partial \left (D^2\psi \right )}{\partial r}+2D^2\psi L\psi \!\right ]=\delta ^2 D^2\!\left (\!D^2\psi \!\right )\!. \end{align}
Here,
$\delta ^2 = \nu / (\omega a^2 )$
, where
$\delta$
represents the ratio of the viscous penetration length scale,
$ (\nu /\omega )^{1/2}$
, to the bubble radius,
$a$
. The kinematic viscosity of a Newtonian fluid is typically
$10^{-6}\,\mathrm {m^2\,\rm s^-{^1}}$
. In experiments, the range of
$\omega a$
is
$O (10^5 ){-} O (10^7 ) \mathrm {\,\unicode{x03BC}m\,\rm s^-{^1}}$
(Yazdi & Ardekani Reference Yazdi and Ardekani2012). The bubble radius is typically set to be
$a=10\,\mathrm {\,\unicode{x03BC}m}$
. Therefore, we can estimate that the range of
$\delta$
is
$O (10^{-1} ) {-} O(1)$
, which means that the viscous penetration length is smaller compared with the radius of the oscillating bubble. Furthermore, the relative oscillating amplitude,
$\epsilon$
, is experimentally measured to be
$O (10^{-2} ){-} O (10^{-1} )$
. Thus, we assume the asymptotic limit
$\epsilon \ll \delta \ll 1$
, which is used to outline the asymptotic analysis in the following section. These estimates of
$\epsilon$
and
$\delta$
are primarily provided to illustrate that the asymptotic conditions can be reasonably satisfied, ensuring that the analytical derivation aligns with plausible experimental parameter ranges. Detailed discussions can be found from Spelman & Lauga (Reference Spelman and Lauga2017).
2.2. Asymptotic solution of the streaming flow
In this section, we outline the main steps to calculate the asymptotic solution of the acoustic streaming flow, reproducing the known results from Spelman & Lauga (Reference Spelman and Lauga2017) to ensure consistency with foundational work and to provide necessary details for the subsequent analysis of particle dynamics. Given that the oscillation of the bubble is of the order of
$\epsilon$
, it is straightforward to expand the flow field in orders of
$\epsilon$
:
\begin{align} &\left \{\psi ,u_r,u_{\varphi }\right \}=\epsilon \left \{\psi ^{(1)}, u^{(1)}_r, {u}^{(1)}_{\varphi }\right \}+\epsilon ^2\left \{\psi ^{(2)}, u^{(2)}_r, {u}^{(2)}_{\varphi }\right \}+o\left (\epsilon ^3\right ). \end{align}
At the order of
$\epsilon$
, (2.9) is expressed as follows:
\begin{align} &\frac {\partial \left (D^2\psi ^{(1)}\right )}{\partial t}=\delta ^2 D^4\psi ^{(1)}. \end{align}
For the boundary conditions, we consider three modes in (2.1):
$(V_0, V_1, W_1)$
. These modes are minimal yet sufficient to generate a streaming flow (Longuet-Higgins Reference Longuet-Higgins1998). By expanding (2.6) and (2.7), we obtain the boundary conditions at the first order:
\begin{align} &\frac {\partial \psi ^{(1)}}{\partial \varphi }\left |_{r=1}=[V_0+V_1\cos \varphi ]\sin \varphi e^{i t}, \nonumber\right. \\ &\frac {\partial \psi ^{(1)}}{\partial r}\left |_{r=1}=-\frac {1}{2} W_1\sin ^2\varphi e^{i t} ,\right.\\ &\left.\left (-\frac {\partial ^2 \psi ^{(1)}}{\partial r^2}+2\frac {\partial \psi ^{(1)}}{\partial r}+\frac {\partial ^2\psi ^{(1)}}{\partial {\varphi }^2}-\frac {\cos \varphi }{\sin \varphi }\frac {\partial \psi ^{(1)}}{\partial \varphi }\right )\right |_{r=1}=0. \nonumber \end{align}
In view of the oscillatory boundary condition (2.12), the first-order solution will be
$\psi ^{(1)} \propto e^{it}$
. By applying separation of variables twice, we can solve for the first-order streamfunction:
\begin{align} \psi ^{(1)}&=e^{i t}\left [V_0(1-\cos \varphi )+\frac {1}{2}\left (\frac {2A\delta ^2\sqrt {r}}{(1+i)^2}K_{\frac {3}{2}}\left (\frac {(1+i)r}{\sqrt {2}\delta }\right )+B\frac {1}{r}\right )\sin ^2\varphi \right ], \nonumber \\ A&=\frac {-1}{K_{\frac {3}{2}}\left (\frac {1+i}{\sqrt {2}\delta }\right )}\left (\frac {1}{\sqrt {2}\delta }(1+i)+1+\sqrt {2}\delta \frac {(i-1)}{4}\right )(W_1+V_1)+O\left (\delta ^2\right ), \\ B&=V_1+\frac {\sqrt {2}\delta (1-i)}{2}(W_1+V_1)-i\delta ^2(W_1+V_1)+O\left (\delta ^3\right ),\nonumber \\ W_1&=-V_1+\sqrt {2}\delta \frac {6}{1+i}V_1+18 i\delta ^2 V_1+O\left (\delta ^3\right ), \nonumber \end{align}
where
$K_{\frac {3}{2}}(r)$
is the modified Bessel function of the second kind of order
$3/2$
. As we can see, the coefficient
$W_1$
is determined by
$V_1$
as a result of the zero tangential stress condition (2.7) imposed on the bubble surface. To obtain a non-zero time-averaged flow, we need to consider the next order. At the second order, the equation becomes
\begin{align} &-\frac {\partial \left (D^2\psi ^{(2)}\right )}{\partial t}+\delta ^2 D^2\left (D^2\psi ^{(2)}\right )\\ &=\frac {1}{r^2}\left [-\frac {1}{\sin \varphi }\frac {\partial \psi ^{(1)}}{\partial r}\frac {\partial \left (D^2\psi ^{(1)}\right )}{\partial \varphi }+\frac {1}{\sin \varphi }\frac {\partial \psi ^{(1)}}{\partial \varphi }\frac {\partial \left (D^2\psi ^{(1)}\right )}{\partial r}+2D^2\psi ^{(1)} L\psi ^{(1)}\right ]. \nonumber \end{align}
This equation cannot be solved explicitly. Instead, we consider its time-averaged solution,
$\langle \psi ^{(2)} \rangle$
, which is independent of time and satisfies
\begin{align} & \delta ^2 D^4\left \langle \psi ^{(2)}\right \rangle\nonumber\\& =\left \langle \frac {1}{r^2}\left [-\frac {1}{\sin \varphi }\frac {\partial \psi ^{(1)}}{\partial r}\frac {\partial \left (D^2\psi ^{(1)}\right )}{\partial \varphi }+\frac {1}{\sin \varphi }\frac {\partial \psi ^{(1)}}{\partial \varphi }\frac {\partial \left (D^2\psi ^{(1)}\right )}{\partial r}+2D^2\psi ^{(1)} L\psi ^{(1)}\right ]\right \rangle . \end{align}
Equation (2.15) can be solved using matched asymptotics. Since the interaction between the bubble and the fluid occurs over a short range, we consider a boundary layer adjacent to the bubble, with a thickness of the order of
$\delta$
. The distance from the bubble is rescaled as
$\eta = (r-1)/\delta$
. When we scale (2.15), the left-hand side is of the order of
$1/\delta ^2$
and the right-hand side can be carefully expanded in orders of
$\delta$
, with the leading term being of the order of
$1/\delta$
. Consequently, the governing equation in the inner region has the following structure:
\begin{align} \frac {\partial ^4\left \langle \psi ^{(2)}\right \rangle }{\partial \eta ^4}&=\left [\delta C_1+\delta ^2 (C_{21}+C_{22}\eta )\right ]e^{-\frac {1+i}{\sqrt {2}}\eta }\sin ^2\varphi +\left [\delta D_1+\delta ^2 (D_{21}+D_{22}\eta )\right ]e^{-\frac {1+i}{\sqrt {2}}\eta }\nonumber \\ &\quad\sin ^2\varphi \cos \varphi +\delta ^2E_2e^{-\frac {1-i}{\sqrt {2}}\eta }\sin ^2\varphi \cos \varphi +O\left (\delta ^3\right ), \end{align}
where
$C_1, C_{21}, C_{22}, D_1, D_{21}, D_{22}, E_2$
are complex constants that are functions of
$V_0$
and
$V_1$
. Correspondingly, we rescale and expand the boundary conditions. The no-slip boundary condition can then be written in the following form:
\begin{align} \frac {\partial \left \langle \psi ^{(2)}\right \rangle }{\partial \varphi }\bigg |_{\eta = 0}&=\left (a_0+\delta a_1+\delta ^2 a_2\right )\sin \varphi \cos \varphi\nonumber\\&\quad +\left (b_0+\delta b_1+\delta ^2 b_2\right )\left (3\cos ^2\varphi -1\right )\sin \varphi +O\left (\delta ^3\right ), \\ \frac {\partial \left \langle \psi ^{(2)}\right \rangle }{\partial \eta }\bigg |_{\eta = 0}&=\left (\delta c_1+\delta ^2 c_2\right )\sin ^2\varphi +\left (\delta d_1+\delta ^2 d_2\right )\sin ^2\varphi \cos \varphi +O\left (\delta ^3\right ), \nonumber \end{align}
where
$a_0, a_1, a_2, b_0,b_1, b_2, c_1, c_2,d_1,d_2$
are complex constants that are functions of
$V_0$
and
$V_1$
. Combining (2.16) and (2.17), the inner flow can be solved order by order. We will match the inner solution to the outer solution. As
$\eta$
increases to infinity, the solution behaves as follows:
\begin{align} \left \langle \psi ^{(2)}\right \rangle &=\left [G_0+G_1\eta +G_2\eta ^2+G_3\eta ^3\right ]\!\sin ^2\varphi +\left (\!H_0+H_1\eta +H_2\eta ^2+H_3\eta ^3\right )\!\sin ^2\varphi \cos \varphi , \end{align}
The constants,
$G_i,H_i,\: i = 0,1,2,3$
, can be expanded as orders of
$\delta$
:
\begin{align} G_i = G_{i0}+\delta G_{i1}+\delta ^2 G_{i2}+O\left (\delta ^3\right ), \: H_i = H_{i0}+\delta H_{i1}+\delta ^2 H_{i2}+O\left (\delta ^3\right ). \end{align}
In the outer region where
$r \gg 1$
, by substituting (2.13) into the right-hand side of (2.15), we find it decays rapidly, of the order of
$r^{-2}e^{-r}$
, due to the appearance of the modified Bessel functions and their derivatives. Consequently, we can neglect these terms and solve the resulting homogeneous equation:
\begin{align} &D^4\left \langle \psi ^{(2)}\right \rangle =0. \end{align}
Considering the form of the inner flow in (2.18) and applying separation of variables to (2.20), the outer flow takes the following form:
\begin{align} &\left \langle \psi ^{(2)}\right \rangle =\left (T_1r^{-1}+T_2r\right )\sin ^2\varphi +\left (T_3r^{-2}+T_4\right )\sin ^2\varphi \cos \varphi , \end{align}
where
$T_i$
(
$i = 1, 2, 3, 4$
) are constants to be determined. Note that we have already applied the boundary condition that the fluid velocity decays to zero at infinity. To match (2.21) to (2.18), we substitute
$r = 1 + \delta \eta$
and expand (2.21) in orders of
$\delta$
. The constants
$T_i$
are also expanded as follows:
\begin{align} T_i = T_{i0}+\delta T_{i1}+\delta ^2 T_{i2}+O(\delta ^3). \end{align}
Comparing the coefficients in (2.19) and (2.22) to the order of
$\delta$
, we have
\begin{align} G_{00} &= T_{10}+T_{20}, G_{01} = T_{11}+T_{21},G_{11} = -T_{10}+T_{20},G_{12} = -T_{11}+T_{21}, \end{align}
\begin{align} &H_{00} = T_{30}+T_{40}, H_{01} = T_{31}+T_{41},H_{11} = -2T_{30},H_{12} = -2T_{31}. \end{align}
In the derivations above, the complex constants,
$C_1, C_2, D_1, D_2, E_2$
,
$a_0, a_1, a_2, b_0, b_1, b_2, c_1, c_2,$
$d_1, d_2$
,
$G_{00}, G_{01}, G_{11}, G_{12}$
,
$H_{00}, H_{01}, H_{11}, H_{12}$
, are given in Appendix A. In the previous derivation, it is readily apparent that the flow field generated by oscillating bubble exhibits the characteristic scale of
$\omega a$
before conducting non-dimensionalisation. Consequently, the dimensional form of the streaming flow can be written as
\begin{align} &u_r=\epsilon ^2\delta \omega a \bar {u}_r,\quad u_{\varphi }=\epsilon ^2\delta \omega a \bar {u}_{\varphi }, \nonumber \\ &\bar {u}_r=2\left (T_{1}r^{-3}+T_{2} r^{-1}\right )\delta ^{-1}\cos \varphi + \left (T_{3}r^{-4}+T_{4} r^{-2}\right )\delta ^{-1}(3\cos ^2\varphi -1), \\ &\bar {u}_{\varphi }=\left (T_{1} r^{-3}-T_{2} r^{-1}\right )\delta ^{-1}\sin \varphi + 2T_{3} r^{-4}\delta ^{-1}\sin \varphi \cos \varphi . \nonumber \end{align}
As previously mentioned, we assume
$\epsilon \ll \delta \ll 1$
. Thus, we set
$V_0=0.4$
and
$V_1=1$
to calculate a representative acoustic streaming flow field. The resulting flow field is illustrated in figure 2(
$a$
), where the colour-coded background represents the magnitude of streamfunction
$\psi$
and the streamlines indicate the flow direction. The steady acoustic streaming flow is characterised by an axisymmetric double vortex structure. Under the same parameter settings, we calculate the time evolution results of the boundary deformations
$R$
and
$\Phi$
over one period, as illustrated in figure 2(
$b$
,
$c$
). The results indicate that the bubble boundary undergoes in-phase oscillations, exhibiting periodic deformation behaviours. In this case, the fluid surrounding the bubble is pushed away along the
$z$
axis and flows in from the
$x$
direction, leading to the occurrence of vortex structure. Additionally, from (2.1), it can be deduced that the bubble undergoes in-phase oscillations only when
$V_0$
and
$V_1$
are taken as real numbers. In this case, the coefficients
$T_{i0}$
(
$i=1,2,3,4$
) of streamfunction
$\psi$
are purely imaginary and the streaming flow is therefore generated at the order of
$\epsilon ^2\delta$
. However, when either
$V_0$
or
$V_1$
is a complex number, the bubble undergoes out-of-phase oscillations, resulting in the leading order of the streaming flow being
$O(\epsilon ^2)$
. Additionally, we select the same
$V_0$
and
$V_1$
to calculate the streaming flow at orders
$\epsilon ^2$
and
$\epsilon ^2\delta$
, which are provided in Appendix A. We find that the streaming flow of order
$\epsilon ^2$
is characterised as a Stokeslet (Longuet-Higgins Reference Longuet-Higgins1998), while vortex structures are generated at order
$\epsilon ^2\delta$
. In our following investigations, we choose
$V_0=0.4$
and
$V_1=1$
to examine the impact of vortex generated by the in-phase oscillation of the bubble on particles in the acoustic streaming flow.
2.3. Motion of an active particle
In this section, we focus on the effect of acoustic streaming flow on the motion of a particle within a two-dimensional framework. We assume that microorganisms respond instantaneously to the time-averaged flow field because of the significant time scale separation between the time-averaging process and the swimming dynamics of the micro-organisms. Specifically, the time scale associated with the bubble oscillations that determine the time-averaged flow field is
$O (10^{-6} ){-}O (10^{-4} )\mathrm {s}$
(Yazdi & Ardekani Reference Yazdi and Ardekani2012), while the time scale of micro-organism motion is
$O (10^{-1} ){-} O (1 )\mathrm {s}$
(Wadhwa & Berg Reference Wadhwa and Berg2022). We introduce a dimensionless quantity
$\alpha =\epsilon ^2\delta \omega a/v_0$
, which represents the ratio of characteristic flow strength
$\epsilon ^2\delta \omega a$
to the intrinsic swimming speed
$v_0$
. We rescale (2.8) and express it in component form as follows:
\begin{align} \frac {{\rm d}r}{{\rm d}t}&=\cos (\theta -\varphi )+\alpha \bar {u}_r, \nonumber \\ \frac {{\rm d}\varphi }{{\rm d}t}&=\frac {1}{r}\sin (\theta -\varphi )+\frac {\alpha \bar {u}_{\varphi }}{r}, \nonumber \\ \frac {{\rm d}\theta }{{\rm d}t}&=\frac {\alpha }{2r}\left (\frac {\partial (r\bar {u}_{\varphi })}{\partial r}-\frac {\partial \bar {u}_r}{\partial \varphi }\right )+\frac {\alpha \gamma }{2}\left [\sin \left (2(\theta -\varphi )\right )\left (-\frac {\partial \bar {u}_r}{\partial r}+\frac {1}{r}\frac {\partial \bar {u}_{\varphi }}{\partial \varphi }+\frac {\bar {u}_{r}}{r}\right )\right .\\ &\quad\left .+\cos (2(\theta -\varphi ))\left (r\frac {\partial }{\partial r}\left (\frac {\bar {u}_{\varphi }}{r}\right )+\frac {1}{r}\frac {\partial \bar {u}_r}{\partial \varphi }\right )\right ]. \nonumber \end{align}
Here, (
$r$
,
$\varphi$
) denotes the position of active particle in the flow field. Hence, we obtained a deterministic model that controls the behaviour of particle motion in the flow field. The terms on the right-hand side of (2.26) represent the deterministic rate of change with respect to
$r$
,
$\varphi$
and
$\theta$
. Based on our estimates, choosing
$\epsilon = O (10^{-2} ){-} O (10^{-1} )$
,
$\delta = O (10^{-1} ){-} O(1)$
and
$\omega a = O (10^5 ) {-} O (10^7 )\mathrm {\,\unicode{x03BC}m\,s^{-1}}$
together with
$v_0=O(10){-} O (100 )\mathrm {\,\unicode{x03BC}m\,s^{-1}}$
leads to a wide range of
$\alpha$
values
$ (O (10^{-1} ){-} O (10^{3} ) )$
. The presence of curved streamlines in the acoustic streaming flow generates additional torques on particles, causing them to intricately align with local streamlines and modify their swimming behaviours. Additionally, the inherent activity of active particles can significantly influence the selection of their swimming direction.
3. Impact of rotational diffusion
3.1. Discrete description
Using a deterministic model (2.26), we have developed an understanding of the swimming dynamics of a single particle in the acoustic streaming flow. However, in real biological environments, many factors, such as random walk behaviour and thermal fluctuations, can influence the movements and distributions of micro-organisms. Randomness plays a crucial role in the spatial heterogeneity of micro-organisms. To gain insights into the effects of random fluctuations, it is necessary to incorporate noise into our deterministic model. Generally, noise primarily affects the motion behaviours of particles by inducing random translational and rotational diffusion. For a swimmer with typical size, rotational diffusion has more significant effects on its movement compared with translational diffusion (Stark Reference Stark2016). Here, we consider the rotational diffusion coefficient
$D_r$
in our model to represent the random noise in the particle’s orientation. To characterise the impact of diffusion on the swimming dynamics of particles, we define the Péclet number as
\begin{align} & Pe=\frac {v_0}{D_r a}. \end{align}
We define
$\mu _r$
,
$\mu _{\varphi }$
and
$\mu _{\theta }$
as drift terms, representing the deterministic rates of change of
$r$
,
$\varphi$
and
$\theta$
in (2.26). The Langevin equation for a particle, incorporating the effect of random rotational diffusion on its orientation, can be written as follows:
\begin{align} \frac {{\rm d}\theta }{{\rm d}t}&=\mu _{\theta }+\sqrt {2 Pe^{-1}}\xi _r(t), \end{align}
where
$\xi _r(t)$
is a random variable that satisfies
$\langle \xi _r(t) \xi _r(t') \rangle = \delta _r(t - t')$
and
$\delta _r$
follows the von Mises distribution. As a reference, for most micro-organisms, the magnitude of the diffusivity is approximately
$O (10^{-2} ){-} O (10^{-1} )\,\mathrm {s^{-1}}$
(Tavaddod et al. Reference Tavaddod, Charsooghi, Abdi, Khalesifard and Golestanian2011; Drescher et al. Reference Drescher, Dunkel, Cisneros, Ganguly and Goldstein2011; Butenko et al. Reference Butenko, Mogilko, Amitai, Pokroy and Sloutskin2012; Koens & Lauga Reference Koens and Lauga2014). As previously mentioned, the velocities of most microorganisms are
$O(10){-} O (100 ) \mathrm {\,\unicode{x03BC}m\,\rm s^-{^1}}$
. By setting
$a \sim 10 \,\mathrm {\,\unicode{x03BC}m}$
, the inverse Péclet number is
$O (10^{-3} )$
–
$O (10^{-1} )$
.
3.2. Continuous description
To describe the statistical characteristics of the particles in the acoustic streaming flow, we introduce the probability distribution function
$P(r,\varphi ,\theta ,t)$
, which represents the probability that the particle is located at position (
$r,\varphi$
) and moving towards
$\theta$
at time
$t$
. The probability distribution function is non-negative everywhere and its total integral is normalised to
$1$
. Here,
$P$
can be obtained by the Langevin equation (3.2). We define an independent Wiener process
$W$
, satisfying
${\rm d}W=\xi _r(t){\rm d}t$
, where the increment
$\mathrm{d}W$
has zero mean and variance
${\rm d}t$
. Therefore, (3.2) can be rewritten as a stochastic differential equation (SDE),
\begin{align} {\rm d}\theta &=\mu _{\theta }{\rm d}t+\sqrt {2 Pe^{-1}}{\rm d}W, \end{align}
where increment
${\rm d}W$
satisfies
\begin{align} &\langle {\rm d}W\rangle =0,\quad ({\rm d}W)^2={\rm d}t. \end{align}
Then, we can obtain the multidimensional It
$\hat {\mathrm{o}}$
formula from the stochastic differential equation (3.3). This formula can be used to transform the stochastic process into an evolution process of probability density. We therefore derive the corresponding Smoluchowski equation as follows:
\begin{align} \frac {\partial P}{\partial t}&=-\frac {\partial }{\partial r}\left (\mu _rP\right )-\frac {\partial }{\partial \varphi }\left (\mu _{\varphi }P\right )-\frac {\partial }{\partial \theta }\left (\mu _{\theta }P\right )+ Pe^{-1}\frac {\partial ^2P}{\partial \theta ^2}, \end{align}
where
\begin{align} \mu _r&=\cos (\theta -\varphi )+\alpha \bar {u}_r,\nonumber \\ \mu _{\varphi }&=\frac {1}{r}\sin (\theta -\varphi )+\frac {\alpha \bar {u}_{\varphi }}{r}, \nonumber \\ \mu _{\theta }&=\frac {\alpha }{2r}\left (\frac {\partial (r\bar {u}_{\varphi })}{\partial r}-\frac {\partial \bar {u}_r}{\partial \varphi }\right )+\frac {\alpha \gamma }{2}\left [\sin \left (2(\theta -\varphi )\right )\left (-\frac {\partial \bar {u}_r}{\partial r}+\frac {1}{r}\frac {\partial \bar {u}_{\varphi }}{\partial \varphi }+\frac {\bar {u}_{r}}{r}\right )\right . \\ &\quad\left .+\cos (2(\theta -\varphi ))\left (r\frac {\partial }{\partial r}\left (\frac {\bar {u}_{\varphi }}{r}\right )+\frac {1}{r}\frac {\partial \bar {u}_r}{\partial \varphi }\right )\right ]. \nonumber \end{align}
4. Results
4.1. Motion of a single particle
We investigate the swimming dynamics of both passive and active particles under identical initial conditions in the acoustic streaming flow. For the numerical results, we solve (2.26) using a fourth-order Runge–Kutta method, as employed in the study by Arguedas-Leiva & Wilczek (Reference Arguedas-Leiva and Wilczek2020), ensuring a sufficiently small time interval to maintain stability and accuracy. The particles are initially positioned at
$r_0=3$
,
$\varphi _0=\pi /5$
, with an orientation angle of
$\theta _0=5\pi /4$
. Upon exposure to the streaming flow, these particles tend to align with the streamlines. Over time, their positions and orientations exhibit periodic variations, indicating that they are trapped in stable bounded orbits. Figure 3
$(a)$
demonstrates the stable bounded orbits for a spherical particle (
$\gamma =0$
) at different flow strengths
$\alpha$
. We observe that the passive particle moves along the streamline passing through its initial position, whereas the stable bounded orbit of the active particle is closer to the vortex centre and has a smaller range compared with that of a passive particle. In this context, the range of a stable bounded orbit refers to the extent of variation in both the radial and polar directions of the orbit, describing its spatial coverage. This difference of orbit indicates that the inherent activity of active particles enables them to cross streamlines to enter the region near the vortex centre, leading to a reduction in the range of their stable bounded orbits.
Stable bound orbits of passive and active particles: (
$a$
) spherical (
$\gamma =0$
); (
$b$
) elongated (
$\gamma =1$
). We consider different flow strengths,
$\alpha =10$
,
$20$
,
$100$
in different colours. The initial position of each particle is
$(r_0, \varphi _0)=(3, \pi /5)$
(black circular and elliptical point) and the initial orientation angle is
$\theta _0=5\pi /4$
(black arrows). The magnitude of the scaled streamfunction is
$|\psi |/\alpha a^2$
.
Additionally, we observe that the stable bounded orbit of an active particle depends on the flow strength
$\alpha$
. When
$\alpha =10$
, the active particle crosses the streamlines and is trapped in the region near the vortex centre. As
$\alpha$
increases to
$20$
, the stable bounded orbit becomes closer to the vortex centre compared with the case of
$\alpha =10$
. However, for
$\alpha =100$
, the range of stable bounded orbits for the active particles increases compared with the case of
$\alpha =20$
. In this context, the acoustic streaming flow dominates the behaviour of the active particle. As the particle moves into regions far from the oscillating bubble, the cross-flow velocity of the streaming flow aligns it with the streamlines, reducing its ability to enter the interior of the vortex, and ultimately trapping it in the stable bounded orbit farther from the vortex centre.
In figure 3
$(b)$
, we demonstrate the stable bounded orbits of elongated particles (
$\gamma =1$
) under the same initial conditions. We find that in both spherical and elongated cases, passive particles move along the streamlines. Elongated active particles are also trapped in the stable bounded orbits near the vortex centre similar to the spherical cases. However, at the given flow strength
$\alpha$
and initial conditions, the range of the stable bounded orbit of an elongated active particle is systematically smaller than that of a spherical particle. This could be due to the fact that elongated particles are more easily able to cross streamlines. According to Jeffery’s orbit (Jeffery Reference Jeffery1922), elongated particles exhibit more significant and complex rotational behaviours compared with spherical particles. Specifically, when elongated particles are in motion, their rotations tend to follow Jeffery’s orbit, which results in greater alignment with the flow direction and facilitates their ability to move across streamlines. Therefore, it is reasonable to expect that elongated particles, when they have velocity, will have a greater tendency to cross streamlines, leading to a smaller overall motion range. However, we hypothesise that elongated particles may migrate towards regions with higher shear rates, particularly near the vortex centre. This is similar to the findings of Nitsche & Hinch (Reference Nitsche and Hinch1997), who reported that elongated particles in shear flows migrate towards regions with higher shear rates, such as near the wall. However, due to the complexity of the flow field in our study, the exact cause of this migration cannot be fully confirmed.
Dependence of the stable motion range of active particle on the flow strength
$\alpha$
. Motion ranges are measured by (
$a$
)
$\Delta r$
along the radial direction and (
$b$
)
$\Delta \varphi$
along the polar direction. These ranges are normalised by the
$\Delta r_0$
and
$\Delta \varphi _0$
, which represents the motion range of the passive particle along the radial and polar direction, respectively. We consider the spherical (
$\gamma =0$
) and elongated (
$\gamma =1$
) active particles under two different initial condition: IC1 (
$r_0=3$
,
$\varphi _0=\pi /5$
,
$\theta _0=5\pi /4$
) and IC2 (
$r_0=4$
,
$\varphi _0=\pi /4$
,
$\theta _0=\pi /4$
).
The stable bounded orbit of an active particle depends on the flow strength
$\alpha$
and the shape parameter
$\gamma$
, indicating that the interaction between the background flow field and the particle can significantly influence the particle’s range of motion. We further investigate the dependence of the stable motion ranges of both spherical and elongated active particles on flow strength
$\alpha$
. We use
$\Delta r$
and
$\Delta \varphi$
to represent the motion range of particles along the radial and polar directions, respectively, when the dynamics reach a steady state. Specifically, once the motion reaches a steady periodic state, we extract the maximum and minimum values of
$r$
and
$\varphi$
during the periodic motion and calculate their differences to obtain
$\Delta r$
and
$\Delta \varphi$
. These ranges are then normalised by the corresponding values for passive particles, denoted by
$\Delta r_0$
and
$\Delta \varphi _0$
. To ensure consistency across all cases, we use identical initial conditions, which eliminates the influence of initial conditions on the results. In figure 4, we present the stable motion ranges
$\Delta r$
and
$\Delta \varphi$
of spherical (
$\gamma =0$
) and elongated (
$\gamma =1$
) active particles under different flow strengths
$\alpha$
. We consider two initial conditions: the first case
$(r_0,\varphi _0,\theta _0)=(3,\pi /5,5\pi /4)$
(IC1) and the second case
$(r_0,\varphi _0,\theta _0)=(4,\pi /4,\pi /4)$
(IC2). We note that although we present results for two different initial conditions (IC1 and IC2), the ratio of motion ranges between active and passive particles remains qualitatively similar in both cases, indicating that the observed trends are not sensitive to the initial conditions. We observe that as the flow strength
$\alpha$
increases, the motion ranges
$\Delta r$
and
$\Delta \varphi$
of the particles initially decrease for all cases. When the flow strength
$\alpha$
reaches a critical value of
$30$
, the motion ranges
$\Delta r$
and
$\Delta \varphi$
of active particles are minimised. As the flow strength
$\alpha$
continues to increase, the stable motion range of the active particle increases until it converges to
$1$
.
This non-monotonic behaviour arises from the interaction between the intrinsic activity of the active particles and the external acoustic streaming flow. At low flow strengths, the particle’s intrinsic velocity dominates its motion and its orientation is random, leading to a relatively broad motion range. At very high flow strengths, the particle’s intrinsic velocity becomes much smaller than the flow velocity, causing its behaviour to closely resemble that of a passive particle. In this regime, the particle’s motion is restricted to follow the streamlines, with its behaviour determined by its initial conditions. However, when the velocity scales of the particle and the flow become comparable, the particle exhibits more intricate dynamics, such as crossing streamlines and remaining confined nearer to the vortex centre. Consequently, the motion range is relatively large at both low and high flow strengths, resulting in the observed non-monotonic dependence of motion range on flow strength.
Additionally, we find that the stable motion range
$\Delta r$
and
$\Delta \varphi$
for active particles is nearly identical under different initial conditions. This indicates that the initial position and orientation of the particle have little effect on its motion range. Moreover, under the given flow strength
$\alpha$
and initial conditions, the stable motion range of an elongated particle is smaller than that of spherical particle. This further indicates that the elongated geometry of the particle enhances its ability to cross streamlines and enter the region near the vortex centre. Therefore, we demonstrate that the interaction between the flow field and active particles has a crucial impact on the stable bounded orbits and motion ranges of active particles.
4.2. Dependence of trapping on the flow strength
We investigate the trapping of passive and active particles by the acoustic streaming flow in the presence of rotational diffusion. In the deterministic model, the stable bounded orbits of spherical and elongated particles are qualitatively similar. Therefore, we primarily focus on the case of an elongated particle in the following analysis. We initialise 1500 particles uniformly distributed within the region
$1\leqslant r\leqslant 10$
and
$0\leqslant \varphi \leqslant \pi /2$
, with their initial orientation angle set randomly within the range (
$0,2\pi$
). Under the influence of the acoustic streaming flow, the particles exhibit two possible behaviours: escaping from the flow field or being trapped within the streaming flow field. Particles are considered to be trapped when their motion stabilises into a bounded orbit within the computational domain, indicating confinement by the streaming flow. However, particles are considered to escape when they cross the boundaries
$r=1$
,
$r=10$
,
$\varphi =0$
and
$\varphi =\pi /2$
. To account for this, we apply the absorbing boundary conditions to remove particles that reach the boundaries
$r=1$
,
$r=10$
,
$\varphi =0$
and
$\varphi =\pi /2$
. We numerically solve (3.2) by using the Euler–Maruyama method to calculate the trajectories until the particles either escape or their motion ranges stabilise. Thus, we can obtain the stable distribution positions of the trapped particles under different flow strengths.
Stable distribution of trapped passive and active particles represented by orange points. We consider different flow strength
$\alpha =10$
,
$\alpha =20$
,
$\alpha =100$
. The corresponding parameters are
$ Pe^{-1}=0.01$
and
$\gamma =1$
. The magnitude of the scaled streamfunction is
$|\psi |/\alpha a^2$
.
Figure 5 illustrates the stable distribution of the trapped passive and active particles in the acoustic streaming flow. The orange points shown in figure 5 correspond to the instantaneous positions of individual particles after the system has reached a statistically steady state, providing a discrete visualisation that complements the continuum probability density functions. We observe that passive particles are distributed dispersedly within the region of closed streamlines. In contrast, active particles exhibit a tendency to aggregate in specific regions of the flow field, with their stable distribution area being significantly smaller than that of passive particles. When
$\alpha =10$
, the distribution region of the particles is relatively confined, with active particles aggregating in regions far from the oscillating bubble. When
$\alpha =20$
, active particles are primarily distributed in the region near the vortex centre. As the flow strength
$\alpha$
increases further, reaching
$\alpha =100$
, the distribution region expands and becomes more dispersed compared with the case at
$\alpha =20$
. These observations indicate that the stable distribution region of active particles depends on the flow strength.
Dependence of trapping rate
$N_{ {trap}}/N$
of active elongated particles (
$\gamma =1$
) on the flow strength.
$N_{ {trap}}$
is the number of trapped particles, while
$N$
is the total number of particles initially located in the streaming flow.
Then, we define the trapping rate as the proportion of active particles that remain confined within the computational domain after their motion stabilises in a steady state. To calculate this, the positions of all particles are tracked throughout the simulation, and those that remain within the domain are identified. The trapping rate is then determined as the ratio of the number of confined particles,
$N_{ {trap}}$
, to the total number of particles initially introduced into the computational domain,
$N$
. In figure 6, we present the trapping rate of active particles under different flow strengths. We observe that the trapping rate of active particles initially increases with the flow strength
$\alpha$
, as the particle’s intrinsic activity allows for stable confinement within the vortex region. When
$\alpha$
reaches a critical value of
$30$
, the number of trapped active particles reaches its maximum, with the majority of active particles being confined within the acoustic streaming flow field. When the flow strength
$\alpha$
exceeds the critical value and continues to increase, the number of active particles trapped in the acoustic streaming flow gradually decreases. As previously mentioned in figure 4, at low flow strength
$\alpha$
, the intrinsic activity of particle dominates its motion, leading to a larger motion range, which increases the likelihood of particles escaping under the influence of rotational diffusion, resulting in a lower trapping rate. As
$\alpha$
increases to intermediate values, the flow strength becomes comparable to particle activity, reducing their motion range and confining them closer to the vortex centre. This confinement minimises escape probabilities and maximises the trapping rate. However, at high flow strengths
$\alpha$
, the flow field dominates the particle’s behaviour and the motion range increases again, particularly near the vortex boundary. The increased motion range, combined with rotational diffusion, results in a higher probability of escape, causing the trapping rate to decline.
4.3. Distribution of the trapped particles
We further investigate the distribution of trapped passive and active particles under different flow strengths. On the one hand, following the same approach as previously described, we initialise 1500 particles uniformly distributed within the region
$1\leqslant r\leqslant 10$
and
$0\leqslant \varphi \leqslant \pi /2$
, with their initial orientation angle
$\theta$
set randomly within the range (
$0,2\pi$
). The boundaries of this calculation region, defined by
$r=1$
,
$r=10$
,
$\varphi =0$
and
$\varphi =\pi /2$
, are set as absorbing to remove particles that reach the boundaries. Then, we calculate the distribution positions of all trapped particles by solving the Langevin equation until the particle distributions converge over time, indicating that a steady state has been reached. Subsequently, the calculation domain is divided into a uniform grid of
$200$
cells in the
$r$
direction and
$100$
cells in the
$\varphi$
direction. We determine the number of particles in each grid cell and normalise these counts by the total number of particles. By averaging the normalised counts over the angular direction for each radial segment and over the radial direction for each angular segment, we can obtain stable radial distribution function (RDF) and angular distribution function (ADF) of particles.
On the other hand, we numerically solve the Smoluchowski equation (3.5) by employing the Crank–Nicolson method. This computation is performed within the region defined by
$1\leqslant r\leqslant 10$
,
$0\leqslant \varphi \leqslant \pi /2$
and
$0\leqslant \theta \leqslant 2\pi$
. The initial probability distribution function
$P$
within the calculation region has a uniform distribution that satisfies the normalisation condition. To maintain consistency with the boundary conditions employed in the Langevin equation, we impose absorbing boundary conditions with
$P=0$
at the boundaries
$r=1$
,
$r=10$
,
$\varphi =0$
and
$\varphi =\pi /2$
. Additionally, a periodic condition is imposed to
$\theta$
. The numerical solution is carried out on a mesh grid of dimensions
$200\times 100\times 400$
, providing sufficient resolution to capture the distribution’s dynamics. The solution is iteratively computed until a steady-state probability distribution function
$P_s$
is obtained, with convergence monitored throughout the process to ensure accuracy. Therefore, the stable RDF and ADF of particles can be calculated by
\begin{align} &\textrm {RDF}(r)=\iint P_s\, {\rm d}\varphi {\rm d}\theta ,\quad \textrm {ADF}(\varphi )=\iint P_s r\, {\rm d}r{\rm d}\theta . \end{align}
Stable distribution of (
$a$
) passive and (
$d$
) active particles obtained from the Langevin model. The distribution of particles are determined by using RDF and ADF: (
$b,c$
) passive particles and (
$e,f$
) active particles, which are compared with results obtained from Smoluchowski equation. The corresponding parameters are
$\alpha =60$
,
$ Pe^{-1}=0.01$
and
$\gamma =1$
.
Figure 7 demonstrates the stable distribution of the trapped passive and active particles at
$\alpha =60$
,
$ Pe^{-1}=0.01$
and
$\gamma =1$
. Additionally, in figures 7(
$b$
,
$c$
) and 7(
$e$
,
$f$
), we compare the stable RDF and ADF of these particles obtained from the Langevin model (3.2) and the Smoluchowski equation (3.5). The consistency between the results from both models confirms the accuracy and effectiveness of our approach. Figure 7(
$a$
–
$c$
) shows that the distribution of passive particles is relatively uniform in the radial direction, with particles primarily distributed in regions of closed streamlines in the angular direction. Passive particles initially located on unclosed streamlines escape from the flow field due to the incoming flow, whereas only particles located on the closed streamlines can remain confined within the flow field and continue moving along streamlines because of their lack of activity. Thus, passive particles are uniformly distributed in the region composed of closed streamlines. However, in figure 7(
$d$
–
$f$
), active particles exhibit a significant concentration in the central region of the vortex, leading to a non-uniform distribution in both the radial and angular direction. Unlike passive particles, active particles can cross the streamlines laterally and enter the region near the vortex centre, regardless of whether their initial positions are on the enclosed streamlines. Therefore, active particles located within the acoustic streaming flow region can be trapped and collected in the smaller region near the vortex centre compared with that of passive particles, resulting in a non-uniform distribution.
Stable (
$a$
) RDF and (
$b$
) ADF of active particles. We consider different flow strength
$\alpha =5,10,30,50,100$
in different colours. Other corresponding parameters are
$\gamma =1$
,
$ Pe^{-1}=0.01$
.
To further elucidate the dependence of the stable distribution of active particles on the flow strength
$\alpha$
, we calculate the radial distribution function and angular distribution function of active particles under different flow strengths within the range
$(1\lt \alpha \lt 100)$
. In figure 8, we demonstrate the typical results for
$\alpha =5$
,
$10$
,
$30$
,
$50$
and
$100$
. We observe that, at
$\alpha =5$
, active particles concentrate in the region further away from the vortex centre in the radial direction. However, as the flow strength
$\alpha$
continues to increase, the concentration position of active particles in the radial direction is closer to the vortex centre and the particles become more concentrated compared with the case of
$\alpha =5$
. When the flow strength reaches a critical value of
$\alpha =30$
, almost all trapped active particles are collected in the vortex centre of the acoustic streaming flow field, where the degree of accumulation reaches its maximum in both radial and angular directions. Thus, for
$\alpha \leqslant 30$
, increasing the flow strength
$\alpha$
enhances the accumulation degree of active particles. However, when
$\alpha$
continues to increase, as in the cases of
$\alpha =50, 100$
, the accumulation degree of active particles decreases in both radial and angular directions. The accumulation behaviour of active particles strongly depends on individual particle dynamics. At very low or high flow strengths, the particle’s motion range is large, leading to weak accumulation. When the flow is comparable to the particle’s velocity, the motion range of particle becomes smaller, causing a more pronounced accumulation. Consequently, accumulation behaviour of active particles shows a non-monotonic dependence on flow strength. Our results show that the interaction between the acoustic streaming flow and particle activity dominates the distribution position and accumulation degree of active particles.
4.4. Effect of noise on the distribution of particles
Generally, rotational diffusion significantly affects the distribution of active particles. Thus, we calculate the radial distribution function and angular distribution function of active particles under
$ Pe^{-1}=0$
,
$0.01$
,
$0.02$
and
$0.05$
. Figure 9 demonstrates the obtained results. We can observe that as
$ Pe^{-1}$
increases, the aggregation degree of active particles in both radial and angular directions decreases. In other words, an increase in rotational diffusion leads to a more dispersed distribution of active particles under a given flow strength
$\alpha$
. This is because the enhanced rotational diffusion increases the randomness of their orientation, reducing the influence of the flow field on their directional movement. As a result, active particles are more likely to diffuse throughout the streaming flow field rather than being concentrated in specific regions, leading to a more dispersed overall distribution.
Stable (
$a$
) RDF and (
$b$
) ADF of active particles at
$\gamma =1$
and
$\alpha =25$
. We consider different Péclet number,
$ Pe^{-1}=0, 0.01, 0.02, 0.05$
.
5. Concluding remarks
We investigate the swimming dynamics of both passive and active particles within the acoustic streaming flow generated by an oscillating bubble. Through the establishment of a deterministic model, we demonstrate the existence of stable bounded orbits for both passive and active particles in the absence of noise. Our findings indicate that a passive particle can only move along the streamline passing through its initial position, while an active particle has the ability to cross streamlines and become trapped in stable bounded orbits closer the vortex centre compared with that of a passive particle. Then, we define the flow strength
$\alpha$
, which is the radio of characteristic flow strength of acoustic streaming flow to the intrinsic swimming speed of particles. We find that the flow strength dominates the stable bounded orbit of the active particle. Therefore, we further investigate the effect of flow strength
$\alpha$
on the stable motion range of active particles. We observe that as the flow strength increases, the stable motion range of an active particle initially decreases. However, after exceeding a critical flow strength, the range of particle motion gradually increases with
$\alpha$
until it becomes identical to that of a passive particle. This is attributed to the interaction between particle activity and the flow field. Furthermore, we demonstrate that the features of both the stable bounded orbits and the trends in their stable motion range variation with flow strength remain consistent for the spherical and slender shape particles.
We further introduce rotational noise to examine how the trapping and distribution of passive and active particles are influenced by flow strength, using the Langevin and Smoluchowski equations. We find that the distribution range and trapping rate of the active particles initially increases and then decreases with increasing flow strength
$\alpha$
, which is consistent with the observations of a single particle. However, active particles exhibit a significant concentration in the central region of the vortex, leading to a non-uniform distribution in both the radial and angular direction. This phenomenon provides insight into the observed accumulation of Escherichia coli in the acoustic streaming flow field (Yazdi & Ardekani Reference Yazdi and Ardekani2012). Additionally, as the flow strength increases, the concentration region of the active particle is closer to the vortex centre, while the accumulation degree of active particles initially increase and then decreases. Additionally, we identify a flow strength that maximises the accumulation degree of active particles within the acoustic streaming flow. The particle accumulation strongly depends on the motion range of individual particle, resulting in the non-monotonic dependence on the flow strength, which represents a key insight from our study. This result has practical implications for microfluidic systems, suggesting that flow strength can serve as a control parameter to modulate particle behaviour for targeted aggregation or dispersion.
We remark several limitations of the current study and explore potential avenues for future research. First, in our study, we assume both the motion of particles and the streaming flow field to be two-dimensional. This simplification facilitates a more straightforward investigation of particle dynamics in the acoustic streaming flow. However, allowing particles to move in the third dimension can potentially reveal new behaviours not observed in the two-dimensional case. The emergence of complex vortex structures in three-dimensional acoustic flows may lead to altered particle aggregation characteristics and the appearance of scattering behaviour (Ahmed et al. Reference Ahmed, Ozcelik, Bojanala, Nama, Upadhyay, Chen, Hanna-Rose and Huang2016). Second, our study considers the flow field and particle motion to be unbounded, whereas real-world microfluidic environments are typically bounded. In typical microfluidic conditions, an oscillating microbubble is ususlly trapped within microfluidic chips, introducing additional boundary conditions such as wall-induced hydrodynamic interactions. These effects may alter the flow field and particle dynamics (Dey et al. Reference Dey, Buness, Hokmabad, Jin and Maass2022), which are not accounted for in our current model. Incorporating reflective or other boundary conditions with bounded flow fields will provide a more comprehensive understanding of particle dynamics in acoustic streaming flows, which is essential for translating theoretical findings into practical applications in microfluidic devices. Third, our analysis is based on an asymptotic solution of the acoustic streaming flow, which requires
$\epsilon \ll \delta \ll 1$
. While this approximation provides a tractable way to capture the essential vortex structure and define a meaningful flow strength scale, it may not perfectly represent real experimental conditions, especially as
$\epsilon$
and
$\delta$
approach
$O (10^{-1} )$
and
$O(1)$
, respectively. In such regimes, higher-order effects, boundary conditions or other complexities may become significant, and more detailed numerical simulations or direct experimental measurements would be required to accurately quantify the flow field and particle dynamics. Finally, we have focused on the motion of non-gyrotactic particles but neglecting inertial effects. Future research can explore the role of inertia and gyrotaxis, especially in micro-organisms like algae, to better understand particle behaviour in acoustic streaming flows and its implications for biological and fluid systems.
Our study provides valuable insights into the hydrodynamic mechanisms underlying the swimming and aggregation behaviours of active particles within acoustic streaming flow fields. These findings extend beyond basic research and have broader implications across various fields. For instance, the focusing mechanisms of particles under varying flow strengths, as demonstrated in our study, offer indirect evidence supporting the existence of microbial heterogeneity within the human body’s environment (Weiss et al. Reference Weiss, Niedermeier, von Strempel, Burrichter, Ring, Meng, Kleigrewe, Lincetto, Hübner and Stecher2023; She et al. Reference She2024), highlighting the intricate interdependence between micro-organisms and their non-biological surroundings, particularly in biofilm formation. Moreover, the aggregation behaviour of active particles in acoustic streaming flows can be harnessed for practical applications such as developing disposable and portable microfluidic devices for rapid microbial detection (Rossi et al. Reference Rossi, Coulon, Ma, Zhang and Yang2023). By optimising flow strength, these devices can collect micro-organisms more quickly to detectable levels, significantly advancing diagnostic technology. In environmental monitoring, similar mechanisms can be applied to efficiently detect pollutant particles in water and air, improving the speed and accuracy of pollution detection. Furthermore, the control of cell aggregation in tissue engineering and the development of high-precision laboratory equipment in microfluidics and acoustic manipulation can benefit from our insights, paving the way for innovations in areas like single-cell analysis and particle sorting.
Acknowledgements
We would like to thank Professor Eva Kanso and Dr Zhiwei Peng for their invaluable discussions and the anonymous reviewers for their valuable comments in improving this manuscript.
Funding
We acknowledge support from the National Key Research and Development Program of China No. 2021YFA1000200, 2021YFA1000201 (to W.T.) and the National Natural Science Foundation of China under grant No. 12372258 (to Y.M.).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Coefficients and flow characteristics
In the section of sketching the asymptotic solution, the relevant complex constants,
$C_1, C_2$
,
$D_1, D_2, E_2$
,
$a_0, a_1, a_2, b_0, b_1, b_2, c_1, c_2,d_1,d_2$
,
$G_{00}, G_{01}, G_{11}, G_{12}$
,
$H_{00}, H_{01}, H_{11}, H_{12}$
, are given as follows:
\begin{align} &C_1=\frac {3\sqrt {2}}{4}(1+i)\bar {V}_0V_1,\quad C_{21}=\left (\frac {39}{2}-24\sqrt {2}\right )\bar {V}_0V_1,\quad C_{22}=-6(1+i)\bar {V}_0V_1,\nonumber \\ &D_1=\frac {3\sqrt {2}}{4}(1+i)\bar {V}_1V_1,\quad D_{21}=\left (-30+36\sqrt {2}\right )\bar {V}_1V_1,\quad D_{22}=-9(1+i)\bar {V}_1V_1,\nonumber\\ & E_2=3i\bar {V}_1V_1,\quad a_0=V_0\bar {V}_1i+\frac {3i}{2}\bar {V}_0V_1,\quad a_1=-\frac {3\sqrt {2}}{2}(1+i)\bar {V}_0V_1, \end{align}
\begin{align} &a_2=39\bar {V}_0V_1, \quad b_0=\frac {27i}{40}\bar {V}_1V_1,\quad b_1=-\frac {3\sqrt {2}}{4}(1+i)\bar {V}_1V_1,\quad b_2=(30-3i)\bar {V}_1V_1,\nonumber \\ &c_1=\frac {11i}{4}\bar {V}_0V_1+\frac {i}{2}V_0\bar {V}_1,\quad c_2=\left [-\frac {1}{4}+\left (6-\frac {39\sqrt {2}}{4}\right )(1+i)\right ]\bar {V}_0V_1,\nonumber \\ &d_1=\frac {11i}{4}\bar {V}_1V_1,\quad d_2=\left [\left (9-\frac {27\sqrt {2}}{2}\right )(1+i)-\frac {27}{20}\right ]\bar {V}_1V_1,\nonumber \\ &G_{00}=\frac {3i}{4}\bar {V}_0V_1+\frac {i}{2}V_0\bar {V}_1,\quad G_{01}=0, \quad G_{11}=\frac {5i}{4}\bar {V}_0V_1+\frac {i}{2} V_0\bar {V}_1,\quad G_{12}=-\frac {1}{4}\bar {V}_0V_1, \nonumber\\ &H_{00}=\frac {27i}{40}\bar {V}_1V_1,\quad H_{01}=0,\quad H_{11}=\frac {5i}{4}\bar {V}_1V_1, \quad H_{12}=-\frac {27}{20}\bar {V}_1V_1,\nonumber \\ &T_{10}=-\frac {i}{4}\bar {V}_0V_1,\quad T_{11}=\frac {1}{8}\bar {V}_0V_1,\quad T_{20}=\bar {V}_0V_1i+\frac {i}{2}V_0\bar {V}_1,\quad T_{21}=-\frac {1}{8}\bar {V}_0V_1,\nonumber \\ &T_{30}=-\frac {5i}{8}\bar {V}_1V_1,\quad T_{31}=\frac {27}{40}\bar {V}_1V_1,\quad T_{40}=\frac {13i}{10}\bar {V}_1V_1,\quad T_{41}=-\frac {27}{40}\bar {V}_1V_1. \end{align}
In addition, we present the steady acoustic streaming flow with
$V_0=0.4$
and
$V_1=1$
, which are used in our study for particle dynamics in figure 10. The separatrix (
$\psi =0$
) is shown as the boundary between the closed and open streamline regions.
Steady acoustic streaming flow generated by the oscillating bubble for
$V_0=0.4$
,
$V_1=1$
. The magnitude of the scaled stream function is
$|\psi |/(\epsilon ^2\delta \omega a^3$
). The red line
$\psi =0$
represents the separatrix, which separates the closed-streamline regions from the open-streamline regions.
Appendix B. Dependency of steady orbits on initial conditions
In this appendix, we include the full time evolution of trajectories for various initial conditions in the
$(r, \varphi )$
space, including conditions IC1 and IC2 explicitly mentioned in the text, at
$\alpha =20$
. As shown in figure 11, the results illustrate the transient dynamics of the active particles, showing how their positions evolve over time before reaching a steady state. We observe that under different initial conditions, the trajectories of active particles undergo varying evolutions. However, the trajectories of captured particles ultimately converge to similar stable bound orbits, which demonstrates the stability of the particle’s bounded trajectories. For certain initial conditions in figure 11(
$d$
), particles may either escape the flow field without being captured or collide with the bubble, which also reveals the complexity of particle dynamics within the acoustic streaming flow field.
Time evolution results of the radial
$r$
and angular position
$\varphi$
of active particles with different initial conditions: (
$a$
) IC1 (
$r_0=3$
,
$\varphi _0=\pi /5$
,
$\theta _0=5\pi /4$
) and IC2 (
$r_0=4$
,
$\varphi _0=\pi /4$
,
$\theta _0=\pi /4$
); (
$b$
)
$r_0=3,3.1,3.5, \varphi _0=\pi /6, \theta _0=0$
; (
$c$
)
$r_0=3.5, \varphi _0=\pi /8,\pi /6,\pi /3,\theta _0=0$
; (
$d$
)
$r_0=3.5, \varphi _0=\pi /8, \theta _0=0,\pi /2,5\pi /4,7\pi /4$
. The dots represent the initial positions of active particles. Arrows represent their initial motion directions, given by
$(\dot {r}(r_0,\varphi _0,\theta _0),\dot {\varphi }(r_0,\varphi _0,\theta _0))$
. The transparent lines represent the detailed time evolution results, while the solid lines represent their stable orbits. The corresponding parameters are
$\alpha =20$
,
$\gamma =1$
.
