- 1Department of Physics, University of California, Merced, CA, United States
- 2Department of Physics and Astronomy, Bucknell University, Lewisburg, PA, United States
We investigate a model for the dynamics of ellipsoidal microswimmers in an externally imposed, laminar Kolmogorov flow. Through a phase-space analysis of the dynamics without noise, we find that swimmers favor either cross-stream or rotational drift, depending on their swimming speed and aspect ratio. When including noise, i.e., rotational diffusion, we find that swimmers are driven into certain parts of phase space, leading to a nonuniform steady-state distribution. This distribution exhibits a transition from swimmer aggregation in low-shear regions of the flow to aggregation in high-shear regions as the swimmer’s speed, aspect ratio, and rotational diffusivity are varied. To explain the nonuniform phase-space distribution of swimmers, we apply a weak-noise averaging principle that produces a reduced description of the stochastic swimmer dynamics. Using this technique, we find that certain swimmer trajectories are more favorable than others in the presence of weak rotational diffusion. By combining this information with the phase-space speed of swimmers along each trajectory, we predict the regions of phase space where swimmers tend to accumulate. The results of the averaging technique are in good agreement with direct calculations of the steady-state distributions of swimmers. In particular, our analysis explains the transition from low-shear to high-shear aggregation.
1 Introduction
The interaction between self-propelled particles and fluid flows is central to many active matter systems, including swimming bacteria [1], Janus particles [2, 3], and microtubule-based active nematics [4]. A key issue for these systems is understanding how the combined effects of fluid advection and self-propulsion determine the macroscopic properties of the active suspension. This is a nontrivial question even for dilute suspensions of active particles subjected to externally imposed, time-independent fluid flows. In this case, connecting the behavior of individual particles to the macroscopic properties of the fluid-particle suspension may be difficult because active particles can exhibit chaotic dynamics in such flows, e.g., vortex flows [5–9].
For sufficiently simple flows—such as time-independent planar shear flows—swimmers do not exhibit chaos, but the connection between the swimmer trajectories and the macroscopic properties of the suspension has remained elusive. Experiments on swimming microorganisms in microfluidic channel flows reveal that the swimmer concentration profile is nonuniform across the channel and dependent on the type of microorganism and the ratio of the swimming speed relative to the flow speed [1, 10]. Intriguingly, for certain species, the shape of the concentration profile exhibits a transition from being peaked on either side of the channel center to being peaked at the center of the channel as the flow speed increases [10]. Because the magnitude of the fluid shear is peaked at the edges of the channel and goes to zero at the center of the channel, this is characterized as a transition from low-shear depletion to high-shear depletion. In Ref. [1], it is shown that in order to capture this key feature of the experimental data, the model for swimmer motion must take into account fluctuations in the swimming direction caused by rotational diffusion or run-and-tumble swimming. In particular, a Fokker-Planck model [11] predicts a transition from low-shear to high-shear depletion in channel flows as a function of swimmer shape and rotational diffusivity, consistent with the experiments [10].
The Fokker-Planck model provides an accurate quantitative description of the swimmer density in the channel flow, but it does not provide a clear mechanism for the transition from low-shear to high-shear depletion. In particular, it obscures the link between the concentration profile and the swimmer trajectories one actually observes in the presence of both fluid flows and rotational noise. It also leaves open the question of how nonuniform distributions in the swimmer’s phase space arise in the first place, given that typical models of swimmer motion in planar shear flows are conservative dynamical systems perturbed by noise [12–14]. Because conservative systems cannot possess attractors, which typically account for the high-density regions of phase space in noisy dynamical systems, the cause of density variations in the swimmer phase space is more subtle. For elongated swimmers in planar shear flows, the dynamics is conservative but non-Hamiltonian, similar to the case of elongated swimmers in axisymmetric vortex flows [15]. Therefore, the dynamics does not locally preserve phase-space volume, and consequently density variations occur even in the absence of noise. When noise is included, however, the swimmer distributions clearly favor certain trajectories over others, an effect that cannot be accounted for by the non-Hamiltonian character of the dynamics.
In this paper, we study a model of an elongated swimmer in a planar Kolmogorov flow and elucidate the connection between the swimmer trajectories with noise and the swimmer density in phase space. We choose to study the planar Kolmogorov flow—that is, a spatially periodic, alternating shear flow—rather than the channel flow, because this allows us to ignore the effect of boundary conditions on the swimmer dynamics, which can be quite complex and system dependent [16]. We calculate the steady-state distributions of swimmers in the flow numerically, which exhibit nonuniform concentration profiles in the cross-stream direction that are similar to those observed in channel flows. In particular, we map out the transition from low-to high-shear depletion as a function of swimmer speed, shape, and rotational diffusivity. Then, taking advantage of a conserved quantity possessed by our model, we derive a reduced drift-diffusion model for the swimmer dynamics in the limit of weak noise. The reduced model allows us to calculate the likelihood of observing particular swimmer trajectories in the presence of noise. Combining this with the slow-down of trajectories in phase-space, we are able to fully explain the nonuniformity of the swimmer steady-state distributions. The predictions of the reduced model are quantitatively accurate in the small-diffusion limit and explain the transition from low-to high-shear depletion in terms of the relative weighting of different swimmer trajectories.
This paper is organized as follows. In Sec. 2, we describe the phase-space structure of the swimmer in the Kolmogorov flow without noise. In Sec. 3, we add rotational diffusion to our model and calculate the probability density and depletion of swimmers in the Kolmogorov flow as a function of swimmer speed, shape, and rotational diffusivity. In Sec. 4, we introduce the reduced model for swimmer dynamics with weak diffusion, and we use the model to explain our observations from the previous section. Concluding remarks are in Sec. 5. Supplementary Appendix A1 contains details of our numerical method for solving the swimmer Fokker-Planck equation.
2 Deterministic Dynamics of a Swimmer in the Kolmogorov Flow
The 2D, laminar Kolmogorov flow is given by
Here the constant U is the maximum fluid speed at y = 0 and y = ±πw, where 2πw is the spatial period of the flow. We refer to y = 0 and y = ±πw as the centerlines of the flow, in analogy with the Poiseuille flow, which has maximum speed on the centerline of the channel. In the Kolmogorov flow, the equations of motion of an ellipsoidal swimmer with position r = (x, y) and swimming direction θ relative to the x axis, as shown in Figure 1, are [5, 9].
where the first equality in Eq. 2c is Jeffery’s equation for an ellipsoidal particle in an incompressible flow [17]. We have used the notation (⋅),x ≡ ∂(⋅)/∂x. The parameter V is the swimming speed, and α = (γ2 − 1)/(γ2 + 1) is the swimmer shape parameter, where γ is the ellipsoid aspect ratio. The parameter α takes values between −1 and 1, with α = 1 corresponding to a rod swimming parallel to its major axis, α = 0 corresponding to a circular swimmer, and α = −1 corresponding to a rod swimming perpendicular to its major axis. A schematic of the model is shown in Figure 1.
FIGURE 1. Schematic of an ellipsoidal swimmer with shape parameter 0 < α < 1 in the Kolmogorov flow (not to scale).
In this paper, we focus on the dynamics in the (y, θ) plane, since the system is translationally invariant in the x direction. Furthermore, we work with the non-dimensionalized version of Eq. 1, which we obtain by normalizing by the velocity scale U and the length scale w. This leads to
with the rescaled swimming speed v0 = V/U. Because Eq. 3 is periodic in both y and θ, we interpret the phase space topologically as a two torus, i.e., a donut. Remarkably, system (3) is conservative, as is the dynamics of ellipsoidal swimmers in 2D and 3D Poiseuille flows [12, 13], spherical gyrotactic swimmers in Kolmogorov flows [14], and spherical swimmers in 2D vortex flows [15]. The conserved quantity may be derived in the standard way. Dividing Eq. 3a by Eq. 3b, we obtain
Separating variables in Eq. 4, we obtain
Eq. 5 may be integrated, which introduces an integration constant Ψ. This implies that Ψ (y, θ) is a constant of motion, and it is given by
When α = 0, i.e., for circular swimmers, Eq. 3 has a Hamiltonian structure
As a result, each trajectory of Eq. 3 lies on a constant-Ψ contour and the phase space is foliated by periodic trajectories almost everywhere. Illustrative phase portraits are shown in Figures 2A–C for α = 0.9 and Figures 2D–F for α = −0.7. The phase portraits are organized around the four fixed points of Eq. 3. For all v0 and all α < 1, the fixed points are (y1, θ1) = (0, 0), (y2, θ2) = (0, π), (y3, θ3) = (π, π) and (y4, θ4) = (π, 0). The fixed points (y1, θ1) and (y3, θ3) are saddles, while (y2, θ2) and (y4, θ4) are centers, and the two fixed points within each pair map into each other by the shift-flip symmetry
FIGURE 2. Phase portraits of Eq. 3 for (A–C) α = 0.9 and (D–F) α = − 0.7. (A) v0 = 0.2. (B) v0 = v*(α) = 0.430. (C) v0 = 0.8. (D) v0 = 0.02. (E) v0 = v*(α) = 0.285. (F) v0 = 0.6. Red dots are the fixed points of Eq. 3. Red curves are the separatrices emanating out of the saddle fixed points at (y, θ) = (0, 0) and (y, θ) = (±π, π). The color of each remaining curve corresponds to a distinct value of Ψ.
of Eq. 3. The centers are surrounded by periodic orbits, as shown in Figure 2. This can either be viewed as a consequence of the conservative nature of the system, or a consequence of the time-reversal symmetries
Indeed, all of the fixed points are invariant under both of the symmetries (8), in particular the centers, which guarantees that they are surrounded by periodic orbits [18].
The time-reversal symmetries (8) also force the invariant manifolds of the saddles to be homoclinic orbits. These homoclinic orbits (red curves in Figure 2) act as separatrices between two topologically distinct types of periodic trajectories (Figures 2A,C,D,F). Either 1) a trajectory is enclosed by the separatrices of a saddle, and hence wraps around the center inside the separatrices, or 2) the trajectory is in between the separatrices of the two saddles, and hence does not wrap around a center. In other words, the homoclinc orbits separate orbits which are 1) contractible on the torus from those which are 2) non-contractible. We refer to the regions enclosed by a saddle’s separatrices as the islands, and the remaining regions as jets.
For sufficiently small v0 (Figures 2A,D), the spatially varying shear suppresses sustained cross-stream drift and promotes rotational drift. The islands contain swimmer trajectories that oscillate around the centerlines, while the swimmer’s orientation oscillates back and forth, never completing a 2π rotation. Meanwhile, the jets contain swimmer trajectories that oscillate in the y direction, never crossing the centerlines, while the swimmer orientation continuously drifts either clockwise or counterclockwise. On the other hand, for sufficiently large v0, (Figures 2C,F), a swimmer’s self-propulsion can prevail and give rise to sustained motion in the y-direction. While the islands are qualitatively similar to the small v0 case, the jets now either transport swimmers in the +y or −y direction, depending on θ. In these jets, a swimmer’s orientation never completes a full rotation, but oscillates around θ = π/2 or θ = 3π/2, respectively.
The transition from rotational drift to cross-stream drift in the jets occurs at a shape-dependent swimming speed v∗(α), at which system (3) exhibits a global bifurcation. The bifurcation takes place when the separatrices of each of the saddle points coincide exactly, as shown in Figures 2B,E. This occurs when the level sets containing the saddles (0, 0) and (π, π) are identical, meaning Ψ (0, 0) = Ψ (π, π). Using Eq. 6, this relation leads to v∗(α), given by
The bifurcation curve is plotted in Figure 3. This transition from rotational to cross-stream drift is an interesting aspect of the swimmer dynamics in the Kolmogorov flow, which has no direct counterpart in channel flows. A similar phenomenon has been observed numerically for gyrotactic swimmers in the Kolmogorov flow [14].
FIGURE 3. Swimmer phase space structure as a function of the swimming speed v0 and shape parameter α, where α = 0 corresponds to circular particles. The red curve is v*(α) (Eq. 9), the set of parameters where the saddles are heteroclinically connected and phase space bifurcates from containing rotational drift trajectories to cross-stream drift trajectories.
3 Stochastic Dynamics and Noise-Driven Aggregation
In both laboratory experiments and nature, real swimmers are subject to fluctuations in their motion which cause their trajectories to deviate from Eq. 3. We consider the situation where the swimming direction Eq. 3b is perturbed by white noise, while the translational motion continues to follow Eq. 3a. This description is a simple model for biological microswimmers, like swimming bacteria. For such swimmers, translational fluctuations are negligible compared to their self-propulsion, but rotational fluctuations are significant due to the fluctuating active forces that propel the swimmer [19, 20] and potentially run-and-tumble behavior as well [1]. Eq. 3 is replaced by the stochastic differential equation
where σ is the dimensionless strength of rotational diffusion. In Eq. 10b, dW is the infinitesimal increment of a Wiener process, meaning for every instant of time t, dW is a zero-mean normally-distributed random variable with variance dt. We denote a particular realization of the Wiener process as
where PeR = U/DRw is the rotational Péclet number.
3.1 Steady-State Probability Density
Of particular interest is the steady-state probability density P (y, θ) of Eq. 10, which has been measured in experiments on the channel flow [1, 10]. The density P satisfies the steady-state Fokker-Planck equation associated with Eq. 10 [21, 22],
For circular swimmers, α = 0 and Eq. 12 is linear in the partial derivatives of P. Hence, the uniform distribution P = const is a solution, and for σ ≠ 0, we suspect that this solution is unique. Therefore, perfectly circular swimmers with rotational diffusion attain a uniform distribution in phase space. This is to be expected, given that the dynamics in the absence of noise is Hamiltonian and thus area-preserving, and the added noise is isotropic. For α ≠ 0 on the other hand, we expect a non-uniform steady-state distribution, due in part to the non-Hamiltonian nature of the deterministic dynamics. The periodicity of Eq. 12 in y and θ means it is straightforward to solve numerically using Fourier transforms. This aspect also makes the Kolmogorov flow more convenient to work with than a channel flow. We describe our numerical method in Supplementary Appendix A1, and our code is available in Ref. [23]1. Example solutions are plotted in Figure 4. The plots are centered on one of the stable equilibria, (y, θ) = (0, π), and the separatrices enclosing each island are shown as the dotted red curves. The density inside the central y ∈ [ −π/2, π/2] region corresponds to rightward fluid flow (in the x direction), and it is related to the density in the leftward-flow region by the symmetry (7). When α ≠ 0, Figure 4 shows that P (y, θ) is highly nonuniform, featuring sharp peaks in certain regions.
FIGURE 4. Steady-state probability densities P (y, θ) for various combinations of σ, v0, and α indicated on each panel. Red dotted curves are the separatrices emanating from the saddles (y, θ) = (0, 0) and (y, θ) = (±π, π). Panels (A–D) are for parameters before the bifurcation from rotational to cross-stream drift, and panels (E) and (F) are after this bifurcation.
For α > 0 (i.e., swimming parallel to the particle’s long axis) and swimming speeds below the bifurcation [v0 < v∗(α), Figures 4A,B,D], the density is mostly concentrated near θ = 0 and π, with the y position of each peak lying close to the separatrices between the islands and the jets. At these angles, the deterministic part of the angular velocity [Eq. 10b] is minimized, meaning the swimmer spends more time near θ = 0 and π. Hence, this explains the increased density around those angles. When the noise is small (Figure 4A), there are three peaks around each island: one at the saddle point at θ = 0 (which appears as two peaks in Figure 4 due to our cutting the torus at θ = 0), and one on each side of the separatrix. Meanwhile, the density in the center of the islands, around the center equilibrium points, is small and at a local minimum. This suggests that the centers are unstable with respect to perturbations from rotational diffusion, which is not a priori obvious since centers are linearly stable. The small-noise distribution (Figure 4A) approximately respects the time-reversal symmetries (8). As rotational diffusion increases (Figure 4D), the peak around the saddle point splits into two peaks, and the distributions increasingly break the time-reversal symmetries. This is evidenced by the asymmetry of Figure 4D with respect to reflections about y = 0 and reflections about θ = π. The increased noise also causes the peaks near the saddle at (y, θ) = (0, 0) = (0, 2π) to begin merging with the peaks along the middle of the separatrix surrounding the island centered on y = ±π. This merging becomes more pronounced for higher v0, because the separatrices get closer together (Figure 4B).
For α > 0 and swimming speeds above the bifurcation (Figure 4E), the density is peaked on both of the saddles, with most of the remaining density concentrated in the cross-stream jet region. The occurrence of peaks on the saddles may be surprising, because these fixed points are linearly unstable and hence one might not expect density to accumulate nearby. Again, the density inside the islands is small. Hence, we can sum up the behavior of α > 0 swimmers as follows. The effect of small noise is to eject swimmers from the centers of the islands and make them aggregate near the separatrices. Increasing the intensity of noise tends to preferentially push the swimmer density into the jet region, as seen when comparing Figures 4A,D, but the density is still strongly peaked at specific angles near θ = 0 and π. The aggregation of swimmer density near the separatrices is consistent with similar behavior seen in Monte Carlo simulations of swimmers in Poiseuille flow [1]. These effects are more pronounced the more elongated a swimmer is, i.e., the closer α is to 1. As α → 0 and the swimmer shape becomes more circular, the distribution P (y, θ) tends towards the uniform distribution, so the aggregation effects gradually diminish.
On the other hand, for α < 0 (i.e., swimming perpendicular to the particle’s long axis), the density is peaked both at the center fixed points and away from the islands at the angles θ = π/2 and 3π/2 (Figures 4C,F). When v0 is small and below the bifurcation, the peaks around the centers have a narrow width in the y direction, as seen in Figure 4C. Meanwhile, there are prominent peaks for θ = π/2 and 3π/2 in the jet region. These angles correspond to the minimal angular speeds in Eq. 10b for the α < 0 case, where the swimmers spend more time. The peaks around the centers suggest that for α < 0, the centers are metastable with respect to perturbations from rotational diffusion. This is a stark difference from the α > 0 case, for which P (y, θ) is at a local minimum at the centers, and it is all the more striking because the linear stability type of the center fixed points does not depend on α. As v0 increases, the peaks around the centers intensify while the peaks at θ = π/2 and 3π/2 gradually fade. For a sufficiently large v0 (Figure 4F), the peaks at the centers completely dominate the probability distribution. Though we have selected a v0 above the bifurcation in Figure 4F, the peaks at the centers become dominant compared to those at θ = π/2 and 3π/2 for values of v0 well below the bifurcation. This behavior is revealed by examining the swimmer cross-stream concentration profiles.
3.2 Cross-Stream Concentration Profile and Depletion
Using the calculated densities P (y, θ), we investigate the depletion of swimmers from low- or high-shear regions of the flow, which is known to exhibit a complex dependence on the swimmer parameters (v0, α, σ). Following prior work [1, 10, 11], we define a depletion index to quantify the nonuniformity of the spatial distribution of swimmers. The spatial distribution is given by
The depletion index ID is defined by computing the ratio of the concentration in the low-shear regions of the flow to the uniform concentration, and subtracting this from 1. For Poiseuille flows, the cutoff for the low-shear region is usually taken to be the central half-width of the channel. In the case of the Kolmogorov flow, which is like a periodically alternating Poiseuille flow, we take a similar definition as the central half-width of a region in which the flow points in one direction. Hence, we define ID as
The first term in Eq. 14 is the concentration in the rightward low-shear region, and the second two terms comprise the concentration in the leftward low-shear region (see Figure 1). For a uniform density, ID = 0; for low-shear depletion of swimmers, ID > 0; for high-shear depletion, ID < 0.
Example results of these calculations are shown in Figure 5 for the parameters used in Figure 4. We observe a diverse set of concentration profiles that reflect the strong dependence of the density P (y, θ) on the system parameters. For α > 0 and small v0 (Figures 5A,D), C(y) is peaked on either side of each centerline as a result of the large peaks along the separatrices seen in Figures 4A,D. Because the islands are very narrow in the y direction when v0 is small, the concentration peaks on either side of the centerline are close together. Hence, we obtain a negative depletion index, indicating high-shear depletion. For higher values of v0, the separatrices [and the peaks of P (y, θ)] move away from the centerlines (Figure 4B), leading to a positive depletion index (Figure 5B). However, when v0 is sufficiently large and past the bifurcation, the accumulation of density on the saddles at y = 0 and ±π (Figure 4E) leads again to a negative depletion index (Figure 5E). For α < 0 and small v0 (Figure 5C), the concentration profile has peaks both at the centerlines and in the high-shear regions, owing to the P (y, θ) peaks at the centers and the jet regions, respectively (Figure 4C). Whether one observes a positive or negative ID depends on the relative prominence of these two sets of peaks in phase space. As v0 increases and the peaks at the centers dominate (Figure 4F), the swimmer concentration only features peaks at the centerlines (Figure 5F), leading to a negative depletion index.
FIGURE 5. Concentration profiles C(y) for various combinations of σ, v0, and α corresponding to those used in Figure 4. Panels A–D are for parameters before the bifurcation from rotational to cross-stream drift, and panels E and F are after this bifurcation. The dashed lines show the uniform concentration profile, to highlight the nonuniformity of the swimmer concentrations C(y). The depletion index ID for each parameter is noted on each panel.
We perform a systematic exploration of the (v0, α, σ) parameter space. The results are plotted in Figure 6 as isosurfaces of ID. The transitions from high-shear to low-shear depletion (ID < 0 to ID > 0) are indicated by the grey surfaces. Clearly, these transitions are highly dependent on v0 and α, though they only depend weakly on σ for the range of σ values considered here. Surprisingly, the transitions appear to be independent of the bifurcation in the swimmer phase space discussed in Sec. 2, indicated by the black dotted curve in Figure 6.
FIGURE 6. Isosurfaces of the depletion index ID in the (v0, α, σ) parameter space. Blue surfaces: ID = −0.15. Grey surfaces: ID = 0. Red surfaces: ID = 0.1. The black dotted curve in the (v0, α) plane is the bifurcation curve shown in Figure 3.
The speed- and shape-dependent depletion behavior of swimmers in the Kolmogorov flow is similar to the experimental observations reported in Refs. [1, 10] for a channel flow. These experiments measured C(y) and ID as a function of maximum flow speed for several swimming microorganisms, with differing shapes and stroke patterns. We qualitatively mimic this type of experiment by taking cuts of Figure 6 at fixed σ and α, leading to ID (v0) as plotted in Figure 7. High v0 corresponds to low flow speed, while small v0 corresponds to high flow speed. Figure 7 shows that in all cases where α > 0, the depletion index is maximized at an intermediate v0. Furthermore, for extremely elongated swimmers such that α is close to 1, ID > 0 for a wide range of v0 (Figure 7C). These observations are consistent with experiments on slender swimming bacteria in channel flows, which found that ID > 0 and is maximized for an intermediate flow speed [1]. On the other hand, for more rounded swimmers for which α is closer to 0, the depletion index goes from positive to negative as v0 decreases, staying negative over a fairly wide range of v0 (Figures 7A,B). This is indicative of a transition from low-to high-shear depletion with increasing flow speed, which has been observed in experiments on motile phytoplankton in channel flows [10]. The most extreme high-shear depletion with increasing flow speed depends on σ and α. For the value of σ used in Figure 7, ID (v0) becomes most negative with decreasing v0 for intermediate values of α (i.e., intermediate aspect ratios, Figure 7B), while the effect is much more modest when α is close to 0 (Figure 7A, aspect ratio close to 1).
FIGURE 7. Depletion index ID dependence on v0 for selected values of α > 0 and σ = 0.1625. (A) α = 0.2. (B) α = 0.5. (C) α = 0.9.
4 Weak-Noise Behavior via the Averaging Principle
In Sec. 3, we showed that rotational diffusion leads to a pronounced, nonuniform probability density of swimmers in phase space. Nonuniform steady states are typical for dissipative—i.e., non-conservative—dynamical systems perturbed by noise. In that case, the probability density is peaked around the dissipative system’s attractors (e.g., stable fixed points or limit cycles) and is shaped by the balance between diffusion and phase space contraction around the attractors. In our case, Eq. 3 is conservative, meaning the phase space does not possess attractors. Hence, the dynamics of Eq. 3 alone cannot explain the regions of phase space where the swimmer density accumulates. This effect is specifically caused by the interplay between rotational diffusion and conservative dynamics in Eq. 10 [1]. In this section, we derive a model that captures this interplay by applying an averaging principle for conservative systems that is valid in the weak-noise limit. The result is a reduced drift-diffusion model that describes how a swimmer randomly moves across the deterministic orbits of its phase space. This model allows us to decompose the swimmer’s steady-state distribution P (y, θ) into the product of a probability density on the space of swimmer orbits, parametrized by the deterministic constant of motion Ψ, and a kinematic factor accounting for the orientation-dependent angular velocity of elongated particles.
4.1 Reduced Drift-Diffusion Model
To gain insight into the interplay between noise and the conservative dynamics, we calculate the evolution equation for the function Ψ (y(t), θ(t)) under Eq. 10. Because Ψ is a function of the stochastic process (y(t), θ(t)), its evolution equation must be derived using Itô’s Lemma [24], which yields the stochastic differential equation
When σ = 0, dΨ = 0, which implies the conservation of Ψ with no noise, as expected. Obviously, noise introduces dissipation, in the sense that Ψ is no longer conserved. When Eq. 15 is integrated, we obtain
Clearly, when σ is very small, Eq. 16 implies Ψ changes very slowly. We accordingly rescale time as t = τ/σ2, under which dW↦σ−1dW, so that Eq. 16 becomes
For sufficiently small σ, a swimmer will complete many oscillations around a periodic orbit at a fixed Ψ (see Figure 2) before its Ψ value will have drifted appreciably. Hence, its motion may be decomposed into fast motion around the deterministic periodic orbits at fixed Ψ, and slow motion transverse to the periodic orbits, caused by noise. In the σ → 0 limit, we can thus approximate the terms on the right-hand side of Eq. 17 by averaging them over one period of the current orbit at fixed Ψ. This averaging principle is derived rigorously for general two-dimensional conservative dynamical systems perturbed by white noise in Ref. [25] and references therein. We describe the technique here and apply it to the swimmer in the Kolmogorov flow.
The averaged equations are as follows. The assumption that Ψ changes slowly implies that the integrand of the first term of Eq. 17 can be approximated by its average value over one period of the orbit with fixed Ψ, which we denote by f(Ψ). The function f is given by
where T(Ψ) is the period of this orbit, and
Meanwhile, a standard result from stochastic processes, proved in Supplementary Appendix A2, is that the stochastic integral that is the second term of Eq. 17 is given by
Eq. 20 essentially states that the effect of the prefactor in front of the noise increment dW is to rescale the time elapsed along the realization of the Wiener process W(τ) by the integrated variance of the noise increment, (∂Ψ/∂θ)2dτ. The integrand on the right-hand side of Eq. 20 can also be approximated using the averaging principle. We define the averaged diffusivity as
Thus, the righthand side of Eq. 20 can be approximated by
Substituting Eqs 19, 22 into Eq. 17, we obtain the approximate drift-diffusion process after averaging,
By time-averaging, we have reduced the two-dimensional drift-diffusion process (10) to the one-dimensional process (23), which describes how a swimmer diffuses across the deterministic orbits.
Eq. 23 is valid in each of the topologically distinct regions of phase space, i.e., the islands and the jets. However, the drift-diffusion processes in each of these regions must be stitched together with proper boundary conditions in order to describe the averaged dynamics on the full phase space [25]. After averaging, the (y, θ) phase space can be reduced to a graph, depicted in Figure 8. Each edge represents one of the topologically distinct regions of phase space, i.e., an island or a jet. Each point on the interior of an edge represents a distinct periodic orbit with a particular value of Ψ. Thus, Ψ is a coordinate along each of the edges. We denote the absolute value of Ψ on the separatrix as Ψs = |Ψ(0, 0)| and the maximum value of Ψ, occurring at the center fixed point (y, θ) = (0, π), as Ψc = Ψ(0, π). By symmetry, the separatrices occur at Ψ = ±Ψs and the centers occur at Ψ = ±Ψc. Thus, the nodes of the graph in Figure 8 at Ψ = ±Ψc represent the centers, and the nodes at Ψ = ±Ψs represent the saddles and separatrices that are on the island-jet boundaries. The graph in Figure 8 has the same structure regardless of whether the jets exhibit cross-stream drift or rotational drift.
FIGURE 8. Reduction of the swimmer phase space to a graph. Each edge of the graph corresponds to the labeled region of phase space, and the nodes correspond to fixed points and separatrices at the boundaries of each region. Ψ is a coordinate along the edges of the graph.
To each edge, we associate a time-dependent probability density pi (Ψ, τ), with i ∈ {1, 2, 3, 4}. The densities pi evolve according to the Fokker-Planck equations associated with Eq. 23,
where Ji is the probability current density
In Eqs 24, 25, fi and Di are the averaged drifts and diffusions [Eqs 18, 21 respectively] evaluated in the regions of phase space corresponding to the edges of the graph in Figure 8. At the nodes of the graph where the islands and jets meet, the local conservation of probability implies that the total probability current density entering a node must equal the total probability current density leaving the node, similar to Kirchoff’s first law for circuits. To obtain the total probability current density entering (leaving) a node, one sums over the Ji for edges i connected to that node such that the node is approached in the direction of increasing (decreasing) Ψ. The boundary conditions are thus
Equations 24, 26 thus constitute a drift-diffusion process on the graph illustrated in Figure 8, which captures the swimmer dynamics in the σ → 0 limit.
In order to investigate the steady-state behavior of the swimmer in the weak-noise limit, we seek the steady-state solution of Eqs. 24–26. Because the jets are identical to each other by symmetry, f2 = f3 and D2 = D3 in Eq. 24, and therefore p2 = p3 must be satisfied in the steady-state. This allows us to merge the distinct pi into a single steady-state density, p0(Ψ), defined over the entire range −Ψc ≤ Ψ ≤ Ψc, which satisfies
The density p0 satisfies the steady-state Fokker-Planck equation
Eq. 28 can be integrated once, which after rearrangement gives
where C1 is an integration constant and (⋅)′ ≡d(⋅)/dΨ. Eq. 29 is linear, and its solution is
where C0 is another integration constant. Due to the symmetry of the problem, the conditions
must be satisfied. We substitute Eq. 30 into Eq. 31a, which gives
where we have used Eqs 31b, 31c in evaluating the first term of Eq. 30 on the right-hand side of Eq. 31a. Eq. 32 implies C1 = 0. Hence, the steady-state Ψ distribution has the form
where we have absorbed another constant into C0, which now plays the role of normalization constant. Eq. 33 is the invariant density of the reduced drift-diffusion model (23). Physically, it is the steady-state probability density of finding a swimmer on a deterministic trajectory with a particular constant of motion Ψ in the σ → 0 limit.
We now use p0 to reconstruct the phase-space distribution P0 (y, θ) that is obtained in the σ → 0 limit. The idea is to make a change of coordinates (y, θ)↦(Ψ, s), where Ψ = Ψ (y, θ) is defined by Eq. 6 and s = s (y, θ) is the elapsed Euclidean arclength along a trajectory at fixed Ψ, defined by
where q = (y, θ) and
Here, P∗(Ψ, s) is the steady-state distribution in (Ψ, s) coordinates. Under the assumptions of the averaging principle, the fast motion along the s coordinate is decoupled from the slow motion along the Ψ coordinate, so the joint density P∗(Ψ, s) must be the product of the density of the Ψ coordinate, p0(Ψ), and the invariant density along an orbit of fixed Ψ. The latter is inversely proportional to the phase-space speed
The orbit period T is included in Eq. 36 to normalize the invariant density of the s coordinate.
The calculation of the determinant in Eq. 35 requires the partial derivatives of Ψ and s with respect to (y, θ). The partial derivatives of Ψ are straightforward to obtain from Eq. 6. The partial derivatives of s follow from Eq. 34,
where ∇ ≡ (∂/∂y, ∂/∂θ). From Eq. 37, it follows that
Thus, a straightforward calculation leads to
Combining Eqs 35, 36, 39, we finally obtain
where
Eq. 40 shows that the nonuniformity of the steady-state probability density P (y, θ) in the limit of small noise results from two effects: rotational slowdown and noise-induced drift. Regardless of where one is in phase space, the density is modulated by the factor of (1 − α cos 2θ)−1. This factor accounts for the slowdown of rotating elongated particles, caused by the orientation-dependent angular velocity in Eq. 3b. In addition to this, the density is modulated by the function g(Ψ), which weights each point (y, θ) according to the deterministic orbit it belongs to, indexed by Ψ (y, θ). The modulation by g accounts for the slow dynamics across the deterministic periodic orbits, encapsulated by the averaged drift-diffusion model (23).
4.2 Unraveling Density Variations and Depletion
We proceed by evaluating p0(Ψ) through the numerical evaluation of Eqs 18, 21, 33 for selected parameters (v0, α) [23]1. The resulting distributions are compared against the distributions p(Ψ) obtained directly from the numerical solutions P (y, θ) to the Fokker-Planck Equation 12. We compute p(Ψ) using
for small ΔΨ. Example results are presented in Figure 9. Note that in each panel of Figure 9, p is plotted for the range −Ψc ≤ Ψ ≤ Ψc. Because Ψc depends on v0 and α, the numerical ranges of p are different in each panel. For sufficiently small σ, we see excellent agreement between p0 predicted by the averaging principle and p obtained from the exact solutions of the Fokker-Planck equation (Figures 9A,C,E,F). As the intensity of noise increases, the two Ψ distributions begin to deviate from each other, as seen when comparing Figures 9A,D. The assumption underlying the averaging principle is that the noise-induced drift across orbits is slow compared to the fast motion around the deterministic periodic orbits. We estimate the time scale of the fast motion using the period T (Ψc) of small oscillations around the centers (y, θ) = (0, π) and (π, 0), which is given by
FIGURE 9. Comparison of the steady-state Ψ distributions from the averaging principle [p0(Ψ), red curves] and those obtained from solving the Fokker-Planck Equation 12 and applying Eq. 42. The parameters are the same as in Figure 4. Panels A–D are for parameters before the bifurcation from rotational to cross-stream drift, and panels E and F are after this bifurcation. Black dotted lines indicate the Ψ values of the separatrices ±Ψs. Note that p0(Ψ) does not depend on σ, so the curves plotted in panels (A) and (D) are the same.
Eq. 43 can also be expressed in terms of the dimensional quantities as
where we have dropped the numerical constants. The scaling behavior in Eq. 43 explains why for the same value of σ = 0.15, the agreement between p0 and p is better for large v0 (Figure 9E compared to Figures 9B,D) and α < 0 (Figure 9C compared to Figure 9D).
Next, we compare the probability densities P0 predicted by the reduced model and given by Eq. 40 to the densities P obtained by solving the Fokker-Planck Equation 12. In particular, we focus on assessing the extent to which the non-uniformity of P is caused by noise-induced drift, as opposed to rotational slowdown. Owing to the form of P in the σ → 0 limit given in Eq. 40, we filter out the effect of rotational slowdown by multiplying the densities by 1 − α cos (2θ). The results are shown in Figure 10 for selected parameters. Figures 10A–C show Fokker-Planck densities P with the rotational slowdown modulation removed, and they correspond to Figures 4A,C,E, respectively. For α > 0 (Figures 10A,C), the peaks at θ = 0 and π are gone, confirming that these peaks are caused by rotational slowdown. Similarly, for α < 0 (Figure 10B), the peaks at θ = π/2 and 3π/2 are gone. Furthermore, the filtered densities P (1 − α cos 2θ) from the Fokker-Planck simulations (Figures 10A–C) agree well with the corresponding filtered densities P0 (1 − α cos 2θ) obtained using the averaging technique in the weak-noise limit (Figures 10D–F). After filtering out the rotational slowdown, all that remains in P0 is g (Ψ(y, θ)). Thus, Figures 10D–F display the orbit weight-factor g at each point of phase space.
FIGURE 10. Comparison of (A–C) the phase-space distributions with the rotational slowdown effect filtered out, P (y, θ) (1 − α cos 2θ), with (D–F) the distributions P0 (y, θ) (1 − α cos 2θ) predicted by the averaging model given in Eq. 40.
Figure 11 shows g(Ψ) for several representative combinations of v0 and α. Note that g(Ψ) appears to be a smooth function, while p0(Ψ) possesses singularities at the separatrices Ψ = ±Ψs, as seen in Figure 9. It can be shown that these singularities occur because in Eq. 33, D (±Ψs) = 0. However, from Eq. 41 we have g ∝ p0/T, and T also diverges at the separatrices. Evidently, these singularities exactly cancel each other out, making g a smooth function.
FIGURE 11. Phase space weight factors g for (A) α = − 0.7 and (B) α = 0.9. The solid curves are g(Ψ), and the dotted curves represent the separatrices at ±Ψs. Different shades represent different swimming speeds v0. (A) Black is v0 = 0.02, and grey is v0 = 0.6. (B) Black is v0 = 0.02, and grey is v0 = 0.8.
The good agreement between Figures 10A–F both validates the reduced model derived in Sec. 4.1 and shows the importance of noise-induced drift in shaping the overall phase-space densities. In particular, the shape of g(Ψ) for different values of swimmer speed v0 and shape α explains many of our observations concerning the steady-state distributions P (y, θ) in Sec. 3. For example, for α < 0, g is always peaked at the boundaries ±Ψc (Figure 11A). In phase space, these features manifest as peaks on the centers, as seen in Figures 10B,E, as well as Figures 4C,F. When v0 is sufficiently small, g also has a secondary, broad peak centered on the jet region (Ψ = 0), as evidenced by the black curve (v0 = 0.02) in Figure 11A. This leads to a significant density in the jet regions, as seen in Figure 4C. As v0 increases, the height of this secondary peak relative to the peaks at ±Ψc decreases, while the local minima on either side of the central peak creep inward. For sufficiently large v0, these local minima merge and the secondary peak disappears entirely, as evidenced by the grey curve (v0 = 0.6) in Figure 11A. Thus, as v0 grows, the orbits in the jet region become increasingly suppressed for α < 0 swimmers (Figure 4F).
We observe qualitatively different behavior for g when α > 0 (Figure 11B). When v0 is small, g is always peaked near the separatrices, as illustrated by the black curve (v0 = 0.02) in Figure 11B. Note that the peaks are slightly offset from the separatrices, shifted towards the islands. This feature explains the bright bands just inside the separatrices seen in the filtered P (y, θ) distributions in Figures 10A,D and the raw distributions in Figures 4A,B,D. As v0 increases, both the separatrices and the peaks of g move inward towards Ψ = 0. For sufficiently large v0, these peaks merge and g exhibits a single peak at Ψ = 0, as illustrated by the grey curve (v0 = 0.8) in Figure 11B. Hence, at sufficiently high v0, the orbits in the jet region are weighted the highest, as seen in Figures 10C,F.
The variation of g with the swimmer parameters also clarifies the variation of the depletion index ID as a function of v0 and α. Figure 12A shows the variation of ID predicted using P0 given by Eq. 40. The quantitative agreement with the depletion index from the Fokker-Planck model is excellent for σ = 0.05 (Figure 12B) and reasonable for σ = 0.15 (Figure 12C). For α < 0, we can now ascribe the transition from low-to high-shear depletion (positive to negative ID) to the increasing dominance of the peaks of g at Ψ = ±Ψc as v0 increases. Meanwhile, for α > 0, we ascribe the transition from high-to low-shear depletion to the peaks of g tracking the separatrices, which narrowly hug the centerlines for small v0 and gradually widen as v0 increases.
FIGURE 12. Depletion index ID as a function of (v0, α) for fixed values of σ. (A) σ → 0, using the averaged drift-diffusion model and Eq. 40. The black dotted curve is the bifurcation curve from Figure 3. (B) σ = 0.05, from the Fokker-Planck model. (C) σ = 0.15, from the Fokker-Planck model.
As for the second transition to high-shear depletion that occurs for α > 0 as v0 increases further, i.e., the blue region in the upper-right corner of Figure 12A, this is actually due to the interplay between noise-induced drift and rotational slowdown. To demonstrate this, in Figure 13 we compare ID calculated using the full P0 given by Eq. 40 (Figure 13A) to a calculation of ID where we use instead the filtered density P0 (y, θ) (1 − α cos 2θ) = g (Ψ(y, θ)), which must be normalized (Figure 13B). The two calculations agree well everywhere except for α near 1 for high v0. In that region, ID decreases as v0 increases in both cases, but with rotational slowdown included (Figure 13A), the decline of ID is enhanced to the point that it becomes negative. For α near 1 and high v0, g is peaked at Ψ = 0 (Figure 11B), so that orbits in the jet region are highly weighted (Figures 10C,F). At the same time, the rotational slowdown factor 1 − α cos 2θ is accentuated, due to a larger |α|. The combination of these two effects leads to a very high concentration of P (y, θ) near the saddles, modest ridges in the jet regions, and a steep drop-off everywhere else, as seen in Figure 4E. These elements combined then lead to ID < 0. Overall, we conclude that noise-induced drift is the dominant effect underlying the (v0, α) parameter dependence of ID, while rotational slowdown is a secondary effect that becomes most important for fast and highly elongated α > 0 swimmers.
FIGURE 13. Comparison of calculations of the depletion index ID (A) using P0 (y, θ) vs. (B) using the rotational slowdown filtered P0 (1 − α cos 2θ) = g.
5 Conclusion
To conclude, we have analyzed the distribution of microswimmers in the planar, laminar Kolmogorov flow. These distributions are highly nonuniform for elongated swimmers, and the swimmer concentration may be depleted from either low- or high-shear regions of the flow [1, 10]. To explain this effect, we derived a reduced model using an averaging technique [25], which captures the slow motion of swimmers transverse to the deterministic orbits that is induced by weak diffusion. The model leads to an invariant density on the space of deterministic swimmer orbits, parametrized by the constant of motion Ψ of the deterministic dynamics. The steady-state phase-space densities predicted by the reduced model are in good agreement with those obtained from the original model, showing that the main cause of depletion is the noise-induced drift of swimmers in phase space.
The averaging technique we employed here can be applied to other problems involving noisy self-propelled particles in fluid flows. For systems with an effectively two-dimensional phase space, the technique is applicable so long as the system possesses a conserved quantity and the noise is diffusive. For example, it may also be applied to gyrotactic swimmers in a planar laminar Kolmogorov flow [14] and swimmers with both translational and rotational diffusion [19, 22]. For higher dimensional systems possessing one or multiple conserved quantities, such as swimmers in three-dimensional channel flows [13], a generalization of the two-dimensional theory exists for extracting the slow dynamics induced by weak noise [25]. While we focused here on a steady flow, the averaging technique could potentially be generalized to an unsteady flow with simple time dependence. Even though we do not expect a conserved quantity in that case, it is possible that there is a fast time scale of the motion that can be averaged out in a similar fashion. It is also an open question whether an averaging approach would be applicable to models where the rotational noise of the swimmer contains a tumbling component in addition to (or instead of) rotational diffusion [22, 26, 27].
The simplicity of swimmer motion in the Kolmogorov flow makes it a good testbed for understanding other aspects of swimmer-flow interaction. For example, it could be used to better understand the effects of chemotaxis on density variations of swimming microorganisms in fluid flows [1, 26, 27]. Our model can be modified by making the noise intensity σ vary with y and θ, in order to mimic the position- and orientation-dependent tumbling rate of chemotactic swimming bacteria in response to a chemical gradient in the flow. This would modify the averaged drift and diffusion functions f and D, defined in Eqs 18, 21, and it would potentially break the symmetry between the two distinct jet regions. However, the rest of the procedure carried out in Sec. 4 would still be applicable.
Finally, it may be possible to exploit the parameter dependence of the density of variations of swimmers in the Kolmogorov flow to sort particles with differing properties.
Data Availability Statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Author Contributions
SB conceived and designed the research. SB, KF, and NB performed the analytical and numerical calculations. All authors contributed to the interpretation of the results and data. SB wrote the first draft of the manuscript. All authors contributed to the revision of the manuscript.
Funding
This study was supported by the National Science Foundation under grants CMMI-1825 379 and DMR-1806 355.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2021.816663/full#supplementary-material
Footnotes
1Github Our code used to solve the Fokker Planck equation and compute the quantities associated with the averaged drift-diffusion model (2021) Available at: https://github.com/saberman52/KolmogorovFlowSwimmerDynamics. (Acccessed: December 3, 2021).
References
1. Rusconi R, Guasto JS, Stocker R Bacterial Transport Suppressed by Fluid Shear. Nat Phys (2014) 10:212–7. doi:10.1038/nphys2883
2. Ebbens SJ, Howse JR In Pursuit of Propulsion at the Nanoscale. Soft Matter (2010) 6:726–38. doi:10.1039/b918598d
3. Bechinger C, Di Leonardo R, Löwen H, Reichhardt C, Volpe G, Volpe G Active Particles in Complex and Crowded Environments. Rev Mod Phys (2016) 88:045006. doi:10.1103/RevModPhys.88.045006
4. Sanchez T, Chen DTN, Decamp SJ, Heymann M, Dogic Z Spontaneous Motion in Hierarchically Assembled Active Matter. Nature (2012) 491:431–4. doi:10.1038/nature11591
5. Torney C, Neufeld Z Transport and Aggregation of Self-Propelled Particles in Fluid Flows. Phys Rev Lett (2007) 99:078101. doi:10.1103/PhysRevLett.99.078101
6. Khurana N, Blawzdziewicz J, Ouellette NT Reduced Transport of Swimming Particles in Chaotic Flow Due to Hydrodynamic Trapping. Phys Rev Lett (2011) 106:198104. doi:10.1103/PhysRevLett.106.198104
7. Berman SA, Mitchell KA Trapping of Swimmers in a Vortex Lattice. Chaos (2020) 30:063121. doi:10.1063/5.0005542
8. Ariel G, Schiff J Conservative, Dissipative and Super-diffusive Behavior of a Particle Propelled in a Regular Flow. Physica D: Nonlinear Phenomena (2020) 411:132584. doi:10.1016/j.physd.2020.132584
9. Berman SA, Buggeln J, Brantley DA, Mitchell KA, Solomon TH Transport Barriers to Self-Propelled Particles in Fluid Flows. Phys Rev Fluids (2021) 6:L012501. doi:10.1103/PhysRevFluids.6.L012501
10. Barry MT, Rusconi R, Guasto JS, Stocker R Shear-induced Orientational Dynamics and Spatial Heterogeneity in Suspensions of Motile Phytoplankton. J R Soc Interf (2015) 12:20150791. doi:10.1098/rsif.2015.0791
11. Vennamneni L, Nambiar S, Subramanian G Shear-induced Migration of Microswimmers in Pressure-Driven Channel Flow. J Fluid Mech (2020) 890:A15. doi:10.1017/jfm.2020.118
12. Zöttl A, Stark H Nonlinear Dynamics of a Microswimmer in Poiseuille Flow. Phys Rev Lett (2012) 108:218104. doi:10.1103/PhysRevLett.108.218104
13. Zöttl A, Stark H Periodic and Quasiperiodic Motion of an Elongated Microswimmer in Poiseuille Flow. Eur Phys J E (2013) 36:4. doi:10.1140/epje/i2013-13004-5
14. Santamaria F, De Lillo F, Cencini M, Boffetta G Gyrotactic Trapping in Laminar and Turbulent Kolmogorov Flow. Phys Fluids (2014) 26:111901. doi:10.1063/1.4900956
15. Arguedas-Leiva J-A, Wilczek M Microswimmers in an Axisymmetric Vortex Flow. New J Phys (2020) 22:053051. doi:10.1088/1367-2630/ab776f
16. Chen H, Thiffeault J-L Shape Matters: A Brownian Microswimmer in a Channel. J Fluid Mech (2021) 916:A15. doi:10.1017/jfm.2021.144
17. Jeffery GB. The Motion of Ellipsoidal Particles Immersed in a Viscous Fluid. Proc R Soc Lond A (1922) 102:161–79. doi:10.1098/rspa.1922.0078
18. Strogatz SH. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Boca Raton: CRC Press (2018).
19. Thiffeault JL, Guo J. Shake Your Hips: An Anisotropic Active Brownian Particle with a Fluctuating Propulsion Force. (2021) arXiv:2102.11758 [cond–mat.soft].
20. Hyon Y, Marcos, , Powers TR, Stocker R, Fu HC The Wiggling Trajectories of Bacteria. J Fluid Mech (2012) 705:58–76. doi:10.1017/jfm.2012.217
21. Solon AP, Cates ME, Tailleur J Active Brownian Particles and Run-And-Tumble Particles: A Comparative Study. Eur Phys J Spec Top (2015) 224:1231–62. doi:10.1140/epjst/e2015-02457-0
22. Berman SA, Mitchell KA Swimmer Dynamics in Externally-Driven Fluid Flows: The Role of Noise (2021) arXiv:2108.10488 [physics.flu-dyn].
23.Github Our code used to solve the Fokker Planck equation and compute the quantities associated with the averaged drift-diffusion model (2021) Available at: https://github.com/saberman52/KolmogorovFlowSwimmerDynamics. (Acccessed: December 3, 2021).
25. Freidlin MI, Wentzell AD. Random Perturbations of Dynamical Systems. New York, NY: Springer (2012). doi:10.1007/978-3-642-25847-3
26. Locsei JT, Pedley TJ Run and Tumble Chemotaxis in a Shear Flow: The Effect of Temporal Comparisons, Persistence, Rotational Diffusion, and Cell Shape. Bull Math Biol (2009) 71:1089–116. doi:10.1007/s11538-009-9395-9
Keywords: self-propelled particles, Fokker—Planck equation, averaging principle, conservative dynamical systems, Brownian motion, microswimmer suspensions, microswimmers dynamics
Citation: Berman SA, Ferguson KS, Bizzak N, Solomon TH and Mitchell KA (2022) Noise-Induced Aggregation of Swimmers in the Kolmogorov Flow. Front. Phys. 9:816663. doi: 10.3389/fphy.2021.816663
Received: 17 November 2021; Accepted: 16 December 2021;
Published: 01 February 2022.
Edited by:
Giancarlo Ruocco, Italian Institute of Technology (IIT), ItalyReviewed by:
Massimo Cencini, Istituto dei Sistemi Complessi, CNR, ItalyCharles Reichhardt, Los Alamos National Laboratory (DOE), United States
Carlo Massimo Casciola, Sapienza University of Rome, Italy
Copyright © 2022 Berman, Ferguson, Bizzak, Solomon and Mitchell. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Simon A. Berman, c2ltb24uYWIuYmVybWFuQGdtYWlsLmNvbQ==; Kevin A. Mitchell, a21pdGNoZWxsQHVjbWVyY2VkLmVkdQ==
†Present address: Kyle S. Ferguson, Medical Physics Graduate Program, Duke University, Durham, NC, United States