Sorting of binary active–passive mixtures in designed microchannels

Abstract

Mixtures of active and passive particles are ubiquitous at the microscale. Many essential microbial processes involve interactions with dead or immotile cells or passive crowders. When passive objects are immersed in active baths, their transport properties are enhanced and can be tuned by controlling active agents’ spatial and orientational distribution. Active–passive mixtures provide a platform to explore fundamental questions about the emergent behaviour of passive objects under simultaneous thermal and active noise and a foundation for technological applications in cargo delivery and bioremediation. In this work, we use computational simulations to study an active–passive mixture confined in microchannels designed with funnel-like obstacles that selectively allow the passage of passive particles. Active particles follow overdamped Langevin translational dynamics and run-and-tumble rotational dynamics. We find that adjusting the tumbling rate of active agents and the microchannel geometry leads to a maximum enhancement of the transport properties of the passive particles (diffusion coefficient and advective velocity) that correlates with the highest mixture sorting efficiency and the shortest response time. We demonstrate that the active drift is the cause of the observed enhanced separation of the mixture in contrast with scenarios where only thermal or active diffusion are present.

---

1 Introduction

The interaction between active and passive agents is ubiquitous and determines critical processes for microorganisms. Parasites moving in crowded environments to spread infection,¹ motile microorganisms swimming in waterbodies’ complex environments to feed and facilitate trophic fluxes,² and the swimming mechanism modulation of Euglena gracilis in crowded media³ are some examples.

From a fundamental point of view, mixtures of active and passive particles constitute a novel class of out-of-equilibrium thermodynamic systems, where the fluctuations shaking the passive species have statistical properties that can be externally tuned at the microscale. From an experimental perspective, active baths can be easily realized with suspensions of motile microorganisms like bacteria^4–7 or microalgae.^8–10 Depending on the passive colloids’ size, the self-propulsion speed of the microorganisms, and the concentration of cells in the active baths, the passive colloids’ diffusion coefficient might be orders of magnitude larger than the thermal one.^4,7,8 Recently, synthetic active–passive mixtures have been studied in the context of the effects of crowders and external fields on active colloids’ motion, showing that the complex interplay of hydrodynamics and chemical fields,^11–13 as well as collisions with passive agents,¹⁴ result in modulations of translational and rotational dynamics of active colloids.

Despite recent progress in probing the behaviour of active–passive mixtures from a theoretical,^15,16 numerical,^17–26 and experimental^4–10,27,28 perspective, understanding their phase behaviour in bulk and under confinement and the mechanism behind the enhancement of the transport properties of passive colloids, is still very much an open challenge.^29–31

The transport properties’ enhancement of passive colloids can be modified by reorganizing the active bath so that the spatial and orientational distribution of active agents are non-homogeneous. Different methods may achieve this reorganization. The first method consists in modifying the shape of passive agents. Experimental and simulation studies have shown that passive elongated agents (polymers or colloidal chains) embedded in active baths induce a non-homogeneous local packing fraction of active particles close to them, resulting, in turn, in localized distributions of the passive agents within the bath and enhanced anisotropic diffusion.^32,33 Designing the passive objects in the shape of spheres,³⁴ rods,^34,35 chevrons^34,36,37 or cogwheels^38,39 induces the accumulation of active agents in the concave regions of the objects, leading to a net translation and rotation, respectively. The net transport induced by the shape might be exploited to develop engines powered by active agents.^40,41 The second method imposes external stimuli on active agents, such as light or spatial gradients of physico-chemical properties. For example, light has been used to control the assembly of active colloids,⁴² tune the motion of active Brownian particles (ABP),⁴³ and control the flow patterns of microalgae suspensions.^44–46

Finally, the third method is confinement, a simple way to generate a space-dependent distribution of active agents, affecting fluctuations and transport of passive particles. This is the method we explore in this article. Recently, it was experimentally demonstrated that spatial confinement leads to a non-uniform steady-state distribution of microalgae that accumulates at boundaries and induces de-mixing of passive particles.⁴⁷ Studies have demonstrated that chevron- or funnel-like obstacles can effectively sort, trap, and separate active particles. For instance, the carnivorous plant Genlisea utilizes its rhizophylls’ funnel-shaped structures to rectify bacterial swimming, directing them toward its digestive vesicle.⁴⁸ At the same time, microfabricated funnel walls have been shown to concentrate swimming bacteria by guiding their motion through the funnels.⁴⁹ Trapping of self-propelled particles with alignment interactions was observed between channels composed of funnel walls,⁵⁰ as well as trapping of self-propelled rods by chevron-like obstacles.^51,52 Similarly, random porous media⁵³ and chiral flowers structures⁵⁴ have been shown effective in trapping and sorting active particles.

By using the methods described above, one can leverage the enhancement of the transport properties and the eventual induced net transport of passive particles to develop technological applications in micro cargo-delivery and bioremediation.^55–58 Most studies on active–passive mixtures consider passive colloids in bacterial baths where bacteria are smaller than passive colloids. Much less studied are passive colloids embedded in baths of eukaryotic microorganisms that can be selected to be larger than the passive colloids. Inspired by the experiments with confined mixtures of the biflagellate microalgae C. reinhardtii and colloidal beads,^8,47 we consider an active–passive mixture whose active agents are larger than the passive species. The size difference might be exploited in sorting or separation processes at the microscale powered by active matter.

It has been shown that the boundaries of either the passive species – when sufficiently large –³⁴ or of the confining space,⁵⁹ can influence the statistical properties of the active bath. Thus, it stands to reason that, in turn, they will also regulate the transport of passive agents. Here we explore the latter case by confining the active–passive mixture in a microchannel with funnel-like obstacles in its central part (Fig. 1). By systematically changing the curvature of the boundary and the diffusivity of the active species, our results reveal the existence of boundary shape parameters that are optimal for the expulsion of passive particles. This phenomenon is then rationalised in terms of the ratio between the obstacles’ size and the run distance of the active agents.

Fig. 1 The system and geometry of the microchannel. Pink particles are active, and their orientation vector is represented in dark red. Blue particles are passive, and yellow particles are frozen and constituents of the channel's walls. The funnel-like obstacles in the central region of the channel have gaps that only allow the passage of passive particles. The active–passive binary mixture is confined in the y-axis by placing walls at y = l_w and y = L_y − l_w. Periodic boundary conditions are applied along the x-axis. (a) The particles are initially placed in random positions in the top compartment. (b) At the steady state, the sorting of the mixture is achieved; the concentration of passive particles in the bottom compartment is much higher than that of the top compartment. (c) The gap between obstacles is Δ = 1.25σ₂₂ and the angle between the obstacles and the x-axis is θ_obs. The length of the obstacles is . The reference for the concave and convex regions is the top compartment of the microchannel. The blue dashed-line rectangles represent the funnels composed of two barriers labeled as “left” and “right”. For the concave gaps, the clogging region for θ_obs ≤ 30° is that consisting of the intersection of the yellow semi-circles of radius, R = σ₃₃/2 + σ₂₂ + σ₁₁/2, in the funnel k, and for θ_obs > 30° the clogging area is that consisting of the area below the straight lines parallel to the barriers at a distance R from the center of the barrier and depicted in dark orange in funnel k. For the convex gaps, the clogging region is defined as the intersection of 2 circles placed at the center of the extreme wall particles with radius, R. All the graphics of trajectories of simulations presented in this article were generated in part using the visualization software Ovito.⁶⁰

The article is organized in the following way. In Section 2, the model and the analysis tools are described. Next, in Section 3, the results are presented and discussed. Finally, in Section 4, we present the conclusions and perspectives for future work.

2 Numerical details

2.1 Model

The system under study consists of a two-dimensional binary mixture of active (type 1) and passive (type 2) particles, confined in a designed microchannel composed of frozen particles (type 3) (see Fig. 1). Active particles (in pink in Fig. 1) are simulated as run-and-tumble particles, whereas passive particles (in blue in Fig. 1) as Brownian particles.

Run-and-tumble is implemented as follows. Active particles obey a Langevin dynamics,

where m_i, r_i, n_i, v_i, are the mass, position vector, unit orientation vector, and velocity vector of particle i, respectively. F_a is the magnitude of the active force that acts along the orientation of the active particles, and γ is the friction coefficient of the implicit solvent in which particles are embedded. For active particles, we set F_a = 10, while F_a = 0 corresponds to the case of passive particles (blue particles in Fig. 1). k_B is the Boltzmann constant and T is the temperature of the underlying thermal bath. We set k_BT = 1 in all studied cases. This generates a thermal force W_i = (W^x_i,W^y_i) acting on the ith particle, given by the standard Gaussian white noise with zero mean and correlation 〈W^x_i(t)W^y_i(t′)〉 = δ_ijδ_xy(t − t′). We choose particles to interact repulsively via a WCA⁶¹ potential,

Here, the indexes α and β run over the three particle types: active, passive, and wall particles. We set σ₁₁ = 2, σ₂₂ = 1, and ε₁₁ = ε₂₂ = 1. The Lorentz–Berthelot mixing rules,^62,63 σ_αβ = (σ_αα + σ_ββ)/2 and , are used when particles of different type interact. The particles constituting the microchannel's walls (type 3, in yellow in Fig. 1) do not interact with each other but interact with the active (type 1) and passive (type 2) particles via the same WCA potential and following the mixing rules with ε₃₃ = 1 and σ₃₃ = 1. Throughout the article, all quantities will be expressed in reduced units using k_BT, σ₂₂ and m₂ as energy, distance, and mass units, respectively.

Since the active and passive particles have different sizes but are embedded in the same solvent, we can write the following relation as a consequence of Stokes’ law, γ₁/γ₂ = σ₁/σ₂ = 2. Thus, we set γ₁ = 10 and γ₂ = 5. With these values of the friction coefficients, the system is close to the overdamped regime: even though inertial effects are present at short times, they are small enough not to be the most important contribution to the system's dynamics (see Fig. S1). Besides, solving the full Langevin's eqn (1) is a more general approach that applies to active–passive mixtures in a broader context where inertial effects are present,^64,65 even though the latter is not the scope of the presented work.

Regarding the rotational dynamics, active particles follow a run-and-tumble motion, a common swimming phenotype observed in bacteria such as E. coli,^66,67 B. subtilis,^68,69 and Enterobacter sp. SM3,⁷⁰ but also, to a first approximation, in eukaryotic microorganisms such as C. reinhardtii.⁷¹ The following discrete stochastic process determines the frequency of tumbling events

where θ_i is the angle formed between the orientation vector n_i of particle i and the positive x-axis. α = t_trial/τ_r, is the tumbling rate, t_trial being the time between trials, and τ_r the mean time between tumbling events or reorientation time. This model uses t_trial = dt, i.e. tumbling trials are performed at every step. ζ ∈ [0, 1] is a uniformly distributed random number and Φ ∈ [0, 2π] is a uniformly distributed random angle.

The system is composed of a total number of N = 3000 particles with N₁ = χ_activeN and N₂ = N − N₁, where χ_active is the fraction of active particles in the mixture, here fixed at 0.50. In Fig. 1, we present the system and the geometry of the obstacles considered in the article. The snapshots are portions of the system to facilitate visualization. As shown in the figure, the simulation box has a rectangular shape with height L_y = 40σ₂₂ along the y-axis and width along the x-axis chosen such as to maintain the packing fraction of the mixture, . The system is initialised with all types 1 and 2 particles randomly placed in the top compartment of the microchannel (see Fig. 1a). We impose periodic boundary conditions along the x-axis, while along the y-axis the system is confined by two parallel walls at y = l_w and y = L_y − l_w. Here l_w = 0.12 is the separation between particles forming the walls (see Fig. 1).

In the central region, the microchannel has funnel-like obstacles featuring gaps that only allow the passage of the passive particles. Under these conditions, the microchannel is elongated in the x-direction and short in the y-axis, thus favoring the particles’ encounters with the obstacles. Fig. 1c illustrates the geometry of the obstacles: their length is set to , and they form an angle θ_obs with respect to the x-axis, which will be varied in the range [0°, 75°] in steps of 15°. The width Δ of the gaps between the obstacles is kept constant in all simulations, at Δ = 1.25σ₂₂. This size is chosen to allow the flow of only one passive particle at a time. The motion of active particles induces a flow of passive particles towards the bottom compartment such that, once steady state is reached, the concentration of passive particles is much higher at the bottom than at the top of the channel (see Fig. 1b).

Simulations are performed using the open-source molecular dynamics package LAMMPS⁷² in which we have implemented the run-and-tumble motion described in eqn (3). The Langevin's thermostat and the NVE integrator are used with dt = 0.001τ. Simulations are run for 20 000τ simulation time units, being . Averages of all computed quantities are taken at steady state. Videos of the system for different values of θ_obs are presented in the SI Video.

2.2 Analysis tools

To monitor the dynamics of the sorting, we calculate the instantaneous fraction of passive particles in the bottom compartment, ,

Here, n₂ is the number of passive (type 2) particles whose position, x₂(τ) at time τ is below the obstacles.

To quantify the passive particles sorted from the mixture and thus the separation efficiency, we compute the average fraction of passive particles in the bottom compartment at steady state:

We sample n_conf = 12 000 configurations to estimate the steady state value of , which are separated by τ = 1 simulation time units (1000 simulation time steps).

One of the factors determining the transition of passive particles across a gap is the likelihood of it being closed by an active particle. We will therefore define a gap as closed whenever an active particle is in a position that, geometrically, would not allow the passage of passive particles between the top and bottom sides of the chamber. The criterion depends on the intersection of exclusion zones defined around each one of the slanted obstacles forming the gaps. An exclusion zone is simply defined as the set of points that are closer than σ₁₁/2 + σ₂₂ to the surface of the obstacle. Whenever the centre of an active particle is within the exclusion region, passive particles cannot fit within the gap between the surface of the active particle and the obstacle. As a consequence, whenever the centre of an active particle is within the area at the intersection of the two exclusion zones of adjacent obstacles, or clogging area A_c, it is impossible for a passive particle to cross the corresponding gap (see Fig. 1c). That gap is then considered as clogged. Examples of A_c's for different values of θ_obs can be found in Fig. S2. Notice that, for large angles, the clogging regions for concave gaps are large enough to host more than one active particle. Since only one active particle is required to clog a gap, we consider the occupied gaps, not the active particles inside the clogging regions. With this information, we can calculate the fraction of open concave and convex gaps, and , as well as the total fraction of open gaps, .

3 Results

Fig. 2 presents the time evolution of the fraction of passive particles in the bottom compartment of the microchannel, , calculated using eqn (4), keeping constant either the tumbling rate (panel a) or the angle between obstacles (panel b).

Fig. 2 The sorting process follows first-order kinetics. The simulation data for (see eqn (4)) are represented with circles. Solid lines represent the fits to of the first-order kinetic model (see eqn (8)). The curves are represented in a semi-log scale in the x axis and have a positive shift of 0.20 in the y-axis for clarity. (a) The fraction of passive particles in the bottom compartment, as a function of time. The active particles’ tumbling rate, α = 1 × 10⁻⁴, is kept constant, whereas the angle of the obstacles θ_obs is varied. (b) as a function of time when the angle of the obstacles, θ_obs = 45° and α is changing.

As shown by the simulation data of the time-evolution of (see circles in Fig. 2), it is safe to assume that the process follows a simple first-order kinetics resulting from the mass balance of the passive particles at the top and bottom compartments of the microchannel:

where is the fraction of passive particles at the top/bottom compartment ; k⁻ is the kinetic constant determining the flow of passive particles from bottom to top and k⁺ from top to bottom. The solution of the system of eqn (6) readswhere the characteristic constant is given by τ_d = 1/(k⁺ + k⁻).

To prove this assumption, we fit the numerical data to eqn (8) using k⁺ and k⁻ as fitting parameters. In Fig. 2, the circles correspond to data from the simulations, and the solid lines are fits to eqn (8), in excellent agreement with the simulation data. Fig. 2a shows the variation of the dynamic response upon changes in the angle of the obstacles, θ_obs, while the tumbling rate is kept constant at α = 1 × 10⁻⁴. Conversely, in Fig. 2b the angle θ_obs = 45° is kept constant while α is varied. Clearly, for both sweeps, the fits capture the system's kinetics quite well. A detailed analysis of the passive–passive case is reported in Fig. S3.

Although the simple first-order kinetic model retains no information on either the geometry of microchannels or the active bath, it is enough to describe the dynamics of the passive particles flowing through the funnels. The contribution of both activity and geometry is encapsulated within the kinetic constants k⁺ and k⁻. Note that due to the asymmetric nature of the system, k⁺ determines the flow of passive particles due to the enhancement of their transport properties when embedded in an active bath, while k⁻ determines the flow of passive particles only due to thermal motion. Still, as we will see, k⁻ depends on the behaviour of active particles, through the changes in opening and closing of the gaps. The results of the kinetic constants k⁺ and k⁻ are reported in Fig. 3a and b. We observe that k⁺ features a maximum that shifts to the left as θ_obs increases, from α ≈ 4 × 10⁻⁴ for θ_obs = 0°, to α ≈ 2 × 10⁻⁴ for θ_obs = 75° (see Fig. 3a). The latter is also the global maximum of k⁺, meaning that the flow of passive particles to the bottom compartment is maximized for those values of α and θ_obs. Conversely, as seen in Fig. 3b, k⁻ increases monotonically with both α and θ_obs.

Fig. 3 The dependence of the kinetic parameters on the activity and geometry. (a) The positive kinetic constant k⁺, (b), the negative kinetic constant, k⁻, (c), the average fraction of passive particles in the bottom compartment at steady state, and (d) the time constant of the sorting process as functions of the tumbling rate of the active particles, α, and the angle of the obstacles θ_obs. In (c) solid lines are predictions of the model in eqn (6) given by and in (a), (b) and (d) are guides for the eye. The inset in (c) is as a function of the persistence length, . The inset in (d) is f(θ_obs,α) as presented in eqn (10). Figure (a)–(c) are presented in a semi-log scale in the x-axis and (d) in a log–log scale.

After briefly discussing the behaviour of k⁺ and k⁻, we now present the results of the steady state passive particles distribution within the two compartments of the microchannel and the time constant of the first-order kinetics. The results in Fig. 3 show, , calculated from simulations according to eqn (5) (panel c) and τ_d (panel d) as a function of α and θ_obs. The rate constants can be used to calculate the fraction of passive particles at steady state, as well as τ_d. Note that the solid lines presented in Fig. 3c correspond to the prediction from the first-order kinetic model in eqn (6) at steady state, . In Fig. 3d, solid lines are just guide for the eye. The fraction displays a maximum for all θ_obs, in the range 1 × 10⁻⁴ ≤ α ≤ 2 × 10⁻⁴, with only a minor dependence on θ_obs. The global maximum, obtained for θ_obs = 75°, is approximately , meaning that the active–passive mixture initially placed at the top compartment is almost completely separated. This represents an improvement of around 30% compared to the passive–passive case, where only steric interactions play a role (; see Fig. S3b). When plotting as a function of the persistence length, , as presented in the inset of Fig. 3c, we observe that the maximum takes place for the range 5σ₂₂ ≤ l_p ≤ 10σ₂₂, which is compatible with the length of the obstacles . This suggests that the maximum sorting is achieved when active particles move a distance of order before tumbling. As mentioned above, the position of the maximum slightly depends on θ_obs. However, the curves overlap for all values of the angle. Thus, we may conclude that the steady state partitioning of passive particles between the microchannel's compartments is largely independent of the shape of the intervening partition.

This shape, however, affects the speed at which the steady state is reached. Fig. 3d shows that the characteristic time τ_d decreases monotonically with increasing θ_obs for all α, and displays a flat minimum in the range 1 × 10⁻⁴ ≤ α ≤ 4 × 10⁻⁴ for all θ_obs. This corresponds to the range of α's in which maximum sorting is achieved. Thus, selecting a tumbling rate in the range 1 × 10⁻⁴ ≤ α ≤ 4 × 10⁻⁴ and θ_obs = 75° provides both maximum sorting and the fastest dynamics. The angle θ_obs influences the timescale τ_d both directly and indirectly. The direct effect has a geometric origin: the fraction of boundary that the passive particles can use to transit between top and bottom halves, depends on θ_obs as . This ratio simply corresponds to the availability of the interphase's boundary for the passive particles to move through it. In other words, is the ratio between the length of the gaps (available to the passive particles to transit between compartments) and the total length of the interphase's boundary. If the timescale τ_d(θ_obs,α) would be only affected by geometrical contributions, the expression, , should exactly hold. Since we observe deviations from the purely geometric case, we propose a factor that accounts for such deviations. The timescale can therefore be written as

where the function f(θ_obs,α) includes all the non-geometrical contributions to the difference between τ_d(θ_obs,α) and τ_d(0,α). The inset in Fig. 3d shows the values of f(θ_obs,α) from simulations. We see that, for the set of parameters explored here, f(θ_obs,α) ≃ 0.0542 ± 0.1157, a small mean value with a large error associated to positive and negative deviations. This is a measure of the indirect influence of θ_obs and α on τ_d beyond the case θ_obs = 0, through their effect on the dynamics of active particles.

One of the main ways by which active particles influence transport of passive ones is by blocking the gaps, hindering the flow of passive particles. To quantify this effect, we calculate the fraction of open gaps as a function of the tumbling rate of active particles. The results are presented in Fig. 4. Note that the steady-state distribution of the active particles within the top compartment of the microchannel is reached in a much shorter time scale than that of the system, which is determined by the steady distribution of passive particles within the microchannel (see Fig. S4), and thus the steady-state fraction of open gaps we are presenting also holds for the transient of the system when a downward net flow of passive particles takes place. For θ_obs = 0, there is no distinction between concave and convex gaps, and thus the curves corresponding to this angle in Fig. 4a–c are the same. These curves set respectively the upper/lower limits of the fractions (Fig. 4(a)) and (Fig. 4(b)) for the other values of θ_obs. As we see in Fig. 4, both and are monotonically increasing functions of the tumbling rate α, showing that active particles with shorter persistence length always obstruct the gaps less than long-persistence ones, regardless of the value of θ_obs. Interestingly, all curves show inflection points for 1 × 10⁻⁴ ≤ α ≤ 2 × 10⁻⁴, marking a transition from a regime in which the gaps are blocked most of the times, to another where the gaps are mostly free. This range corresponds with maximum sorting (see Fig. 3c), highlighting the connection between gaps’ dynamics and passive particle segregation. From this behavior, we may already infer that the most effective sorting relies more on the way active particles enhance the transport of passive particles than on the availability of gaps for the passive particles to flow towards the bottom compartment.

Fig. 4 Dependence of the available cavities for the passive particles to flow on the activity and geometry. (a) The fraction of open concave cavities, . The inset contains the clogging area, A_c, as a function of θ_obs. Details on how to estimate the clogging area can be found in Fig. S2. (b) The fraction of open convex cavities , and (c) the total fraction of open cavities, as functions of the tumbling rate, α and the angle of the obstacles, θ_obs.

Fig. 4a shows that, in contrast to the dependence on tumbling rate, is a monotonically decreasing function of θ_obs, a behaviour correlated with the corresponding increase in the clogging area A_c (Fig. 4a, inset). This is to be expected, as increasing θ_obs – and therefore A_c – leads to a larger set of configurations for single active particles to be able to clog concave pores. This is compounded by the tendency of active particles to accumulate at the concave regions, forming a multilayered structure that becomes more pronounced as θ_obs increases and α decreases (see Fig. S5). The situation is opposite for , which is a monotonically increasing function of both α and θ_obs (Fig. 4b), despite the fact that the clogging area of convex gaps is independent of θ_obs. The cause is the increase in the fraction of active particles within the concave regions as θ_obs increases. On the one hand, this leads to fewer active particles available to close the convex gaps. On the other, the higher concentration of active particles in concave regions leads to more frequent interactions capable of removing active particles from within the clogging areas of convex regions. Altogether this leads to a relative increase in with θ_obs. It is interesting to notice that, for α ≳ 3 × 10⁻³, ceases to depend on either α or θ_obs. This is due to the rapid change in orientation and the fact that, contrary to the concave case, active particles do not form layered structures around convex gaps. Finally, Fig. 4c shows that, contrary to and , the total fraction of open gaps, is largely independent of θ_obs. Deviations become more prominent for θ_obs = 75°, mostly due to the increased tendency of active particles to form multilayered structures at concave gaps. The small deviations we do observe can be either positive or negative depending on α, with a crossover around α ≃ 6 × 10⁻⁵, corresponding to a persistence length, l_p ≈ 17σ₂₂, which correlates with the size of the top compartment of the microchannel, L_y/2 = 20σ₂₂. An alternative way of understanding and is by considering the average time convex and concave gaps remain clogged (t^c_concave and t^c_convex) and unclogged (t^u_concave and t^u_convex). The expressions are and , as shown in Fig. S6.

The enhancement of the transport properties of passive particles within the active bath is a direct consequence of the encounters between active and passive particles. It can be characterised by looking at the time evolution of the probability distribution function (PDF) of passive particles displacements. We hereby will focus on the transport along the y-direction, which is the one relevant for particle sorting. More information on the dynamics along the x-direction – the channel axis – can be found in the SI. We calculate the PDF of displacements Δy within a time interval of Δt, in the top and the bottom sections, P_top/bottom(Δy,Δt) (Fig. S7) and extract both the first and second moments of these distribution as a function of time and position. The passive particles’ effective diffusivity in the top/bottom halves, D^top/bottom_y, can be calculated from the variance of P_top/bottom(Δy,Δt), which is the (de-drifted) mean-squared displacement of the particles (Fig. S7). This should be considered at times intermediate between the typical interval between interactions with active particles and the time required to diffuse a distance of order L_y.

Fig. 5 shows the ratio D^top_y/D^bottom_y, where the bottom dynamics is due to thermal diffusion of confined passive particles. We observe that the ratio of diffusion coefficients in the y-direction reaches a maximum in the range 1 × 10⁻⁴ ≤ α ≤ 2 × 10⁻⁴, which is compatible with those reported for maximum sorting (Fig. 3c), fastest time response (Fig. 3d) and inflection point of (Fig. 4c). This confirms that the best separation performance of the system correlates with the maximization of the activity-induced diffusivity along the y-direction within the top compartment of the microchannel.

Fig. 5 Enhancement of the diffusion coefficient of passive particles in the active bath. The ratio between the diffusion coefficients in the top and bottom compartments along the y direction, D^top_y/D^bottom_y, as a function of the tumbling rate, α, and the angle of the obstacles, θ_obs.

Besides the variations in diffusivity, the interactions with active particles in the top section generate a surprising net drift of passive particles, v^top_y, directed downward. This drift can be derived from the first moment of P_top(Δy,Δt) and it is a consequence of the asymmetry between the two boundaries confining the active particles, as confirmed by simulations (Fig. S7–S9).

Fig. 6 shows the dependence of |v^top_y| on both α and θ_obs. |v^top_y| reaches a maximum at θ_obs = 75° and α ∈ [1 × 10⁻⁴, 2 × 10⁻⁴], corresponding to the maximum sorting and shortest characteristic time τ_d.

Fig. 6 Drift velocity in the y-direction induced by the concomitant effect of activity and confinement. The drift velocity in the top compartment of the microchannel as a function of the tumbling rate, α, and the angle of the obstacles, θ_obs.

Therefore, active particles not only enhance the diffusion coefficient of passive ones but also induce an active drift that – as we will demonstrate – plays a crucial role in enhancing the flow of passive particles through the gaps.

3.1 A minimal one-dimensional advection–diffusion model

So far, we have focused our attention on studying the effect of activity and geometry on the active and passive particles sorting. To microscopically understand the transport mechanisms of the passive particles in the confined system (due to the presence of the active ones), we present a simple advection–diffusion (AD) model similar to the one proposed in ref. 47. We consider a one-dimensional system in the domain y ∈ [0,ρL_y], with the factor ρ fixed a posteriori as indicated below. Within this domain, the concentration of passive particles, c(y,t), evolves in time according to

Here D(y) is the position-dependent diffusion coefficient and v(y) is the effective drift velocity of the passive particles. As mentioned earlier, their values in the top/bottom sections, D_top/bottom and v_top/bottom, are calculated from the first two moments of the displacements PDFs for the passive particles. In the following, we will always consider v_bottom = 0. The two values will then be heuristically interpolated using the sigmoids.

with k = 100, which provides a steep enough sigmoid function without compromising the convergence of the numerical method used to solve the eqn (11).

Particle concentration is ensured through no-flux boundary conditions at both ends of the domain:

The initial concentration field c(y,0) is localised in the top part of the container and it is normalised to 1.

Fig. 7(a)–(c) reports the time evolution of the concentration profile calculated with a conservative explicit finite difference method.^73,74 When both type 1 and type 2 particles are passive (Fig. 7a), corresponding to D_top/D_bottom = 0.75, simulations show that (see Fig. S3b and g). Thus, a smaller diffusion coefficient at the top in the passive–passive case is related to the space occupied by the big passive particles reducing the available space for the small passive particles to diffuse. To account for the effect of the reduced space for the small particles to diffuse in the top compartment due to the presence of the bigger particles in a simple way, thus reproducing the value of from simulations with the AD model, we set ρ = 0.73. Note that for ρ = 1, the solution of the AD model would predict an even distribution of small particles of 50% in each compartment. In this case, although D_top ≠ D_bottom, we assume that ∂D/∂y = 0 since there is no activity in the top compartment (see SI, eqn (S1)–(S6), and Fig. S10).

Fig. 7 Advection–diffusion model. In all the results ρ = 0.73 and θ_obs = 75°. The time evolution of the concentration profile: (a) pure thermal diffusion case (D_top/D_bottom = 0.75, v_top = 0, ∂D/∂y = 0); (b) thermal diffusion + active diffusion (D_top/D_bottom = 1.30, v₀ = 0, ∂D/∂y ≠ 0); (c) thermal diffusion + active diffusion + active advection for v_top = −0.056, ∂D/∂y ≠ 0. (d) The best fitted values of the top drift velocity, v^opt_top, as a function of the tumbling rate, α. (e) Comparison between the steady state concentration of passive particles in the bottom compartment calculated from simulations and the A–D model. (f) Comparison of the Péclet number in the top compartment calculated in the simulations and the A–D model according to . (g) The steady state concentration profiles obtained in the A–D model with the best fitted drifts for all the values of α considered. (h) The time response of the concentration of passive particles in the bottom compartment, calculated from the A–D model, , following first order dynamics as in the simulations.

In the active–passive case (Fig. 7b and c), for each value of the tumbling rate α we start by fixing D_bottom = 1 while D_top is set to the value obtained from simulations (see Fig. 5). In this case, ∂D/∂y is given by eqn (S3) of the SI. Note that we also rescale L_y = 1.0, which jointly with the D_bottom = 1.0 define the time scale τ_AD = L_y²/D^bottom_y. We then determine the advection velocity that minimizes the discrepancy between the steady state concentration in the bottom side of the compartment evaluated with the model and the reference value obtained from simulations. The results of the optimal drift velocity, |v^opt_top|, as a function of α are presented in Fig. 7d. Fig. 7e shows that, for an appropriate choice of v_top, the model in eqn (11) reproduces well the steady state separation observed in simulations. The key nondimensional parameter here is the ratio Pe_top = |L_yv_top/2D_top|, which is akin to a Péclet number for the active displacements of the type 2 particles in the top container. Comparing the values of Pe_top from simulations and for the advection–diffusion model (Fig. 7f) we see a discrepancy by a factor ≲4.5, which is reasonable for a minimal model that distills the key dynamics of a complex process in a simple and intuitive way. The discrepancy might be related to some contributions to the system's dynamics that we are not taking into account in the AD model such as the clogging of the obstacles gaps, the actual shape of the interphase given by that of the obstacles, the diffusion in a crowded media experienced by the passive particles in the top compartment due to the presence of the active particles, and the physical extension of the particles.

Having determined the model parameters that best represent the simulation steady state, we go back to the full AD model and use it to explore which features of the active displacements of passive particles contribute most to particle segregation. We consider the case corresponding to maximum segregation (α = 1.25 × 10⁻⁴, θ_obs = 75°), leading to 〈F_bottom〉 = 0.92 at steady state. Fig. 7a–c show the time evolution of c(y,t) from eqn (11) for the cases of thermal diffusion (v(y) = 0, D_top/D_bottom = 0.75, ∂D/∂y = 0, in blue in Fig. 7h), active diffusion (v(y) = 0, D_top/D_bottom = 1.30, in red in Fig. 7h), active diffusion and active drift (v_top = −0.056, D_top/D_bottom = 1.30, in orange in Fig. 7h). Although in all three cases the time evolution of follows a first-order kinetics (see Fig. 7h), the levels of passive particles segregation achieved at steady state are substantially different. In the first case (thermal diffusion, blue curve), at steady state we obtain , which agrees with the simulation results of the passive–passive case as a consequence of our choice of ρ = 0.73. Upon considering the space-dependent diffusion profile (active diffusion, red curve), is increased to 0.732 due to the flux induced by the diffusivity gradient. Finally, considering the full model with v_top = −0.056 (active diffusion and active drift, orange curve), leads to . This last value is in good agreement with the sorting measured in simulations (see Fig. 3).

These results suggest that the main physical origin of the improved sorting compared with the passive–passive case is the advection induced by the active bath rather than the difference in diffusivity between top and bottom containers. This is confirmed also by ad hoc simulations in which the obstacles are replaced by a simple flat reflective wall affecting only active particles (Fig. S8). Although the system generates a non-uniform diffusivity profile, the symmetry between walls affecting active particles leads to the absence of a net advection (v(y) = 0, see Fig. S9). In the presence of two flat walls, the active particles local density profile is symmetric within the top compartment, with two peaks of the same height close to the walls. Thus, the currents of active particles moving upwards and downwards balance out and, in turn, there is no net drift in the passive particles motion. The results show a trapping efficiency that varies only within ∼5% of the passive–passive case, in line with the above findings for the active diffusion (see Fig. 7b and h).

The role of advection can be further rationalised by taking a closer look at the steady state solution c(y) of the AD model.

where A is a normalisation constant. The fraction will therefore be equal to

Here the rightmost expression has been calculated by approximating the profiles in eqn (12) and (13) as piecewise constant. Using the simulation results in Fig. 5 and 6 it is easy to show that, in the present system, the curve is almost entirely dictated by v(α) (Fig. S11, SI). The peak in separation performance occurs therefore as a consequence of a peak in the advection velocity. In turn, this peak corresponds to the range of tumbling rates α for which the persistence length of the active particles is of the order of the length of the obstacles, , as well as to the inflection point in the fraction of open gaps at steady state (Fig. 4c), which separates the regime where the gaps are mostly closed (lower α) from the one where they are mostly open (larger α).

Altogether, the peak in advection speed can be related to a balance between two competing effects. On the one hand, increasing the tumbling rate increases the rate at which active and passive particles interact, due to a lower accumulation of active particles at the boundaries. On the other hand, decreasing the tumbling rate increases the expected displacement of passive particles due to their interaction with active ones. The peak happens at an optimal value that balances frequency of interactions and active transport per interaction event. The advection speed in the top compartment induced by the concomitant effect of confinement and activity constitutes the principle of an active-pumping mechanism that facilitates the separation of the passive species and might be exploited in technological applications such as cargo delivery and micro-separators.

4 Conclusions

We have studied by means of computer simulations the sorting dynamics of a two-dimensional active–passive mixture within designed microchannels comprising confinement in the y direction and funnel-like obstacles in their central region. The funnels feature gaps which only allow the passage of passive particles, so they behave like a semi-permeable membrane to that effect. Simulations were carried out in the overdamped regime to neglect inertial effects but solving the complete Langevin equations for the translational dynamics, for the sake of generality. Active particles followed a run-and-tumble motion for the rotational dynamics. We have observed that the dynamics of the sorting process follows a first-order kinetic model, allowing us to study the dependence of the kinetic constants with the activity and the geometry of the microchannels.

We have observed that there is a small range of tumbling rates that guarantees the highest fraction of passive particles separated from the mixture and the fastest dynamics. The existence of the most favorable range can be explained in terms of an interplay of (1) maintaining a reasonable fraction of open gaps for passive particles to flow (see Fig. 4c), maximizing the active kicks experienced by passive particles, which in turn enhance passive particles’ diffusion coefficient (see Fig. 5) and induce a negative advective drift (see Fig. 6 and Fig. S6) and the (2) geometry of the funnels which allows the maximum relative space for the flow of passive particles (see Fig. 3b).

Solving a minimal advection–diffusion model, we determined that the flow of passive particles towards the bottom compartment and their subsequent sorting requires a negative active drift velocity in the top compartment. Such sorting is neither observed only by virtue of the difference of the diffusivity between compartments nor by the active diffusion associated to the drift term as the gradient of the diffusivity (see Fig. 7).

Finally, we rationalise the position and the existence of the peak in the advective speed, and in turn the system's peak performance, by the balancing of two competing effects determined by the tumbling rate: the increase of the frequency of encounters between active and passive particles (number of active-kicks) as α increases and the active transport per interaction event (the length of active displacements) as α decreases. The interplay of these effects is the basis of an active-pumping mechanism that might be useful in the development of cargo delivery systems and separators of micro-contaminants.

Future work will consider including hydrodynamic interactions to understand its role on the sorting dynamics and an experimental realization of the system using mixtures of C. reinhardtii and passive beads confined in microfluidics devices with funnel-like geometries. The results presented here might help developing strategies based on microorganisms’ taxis to induce a synchronized motion that might lead to a further increase of the advective velocity (synchronous active-pumping) improving even further the sorting of passive particles in the designed microchannels. Additionally, it might be interesting to study the effects of χ_active on the system's dynamics to discern the role of the multi-layered structure formed on the gaps in their clogging.

Another possible interesting study would be addressing the effects of varying the size ratio between particle types on the sorting dynamics. A previous numerical work in bulk²¹ reports sub-diffusive behavior of the passive particles at large size ratios and, in general, smaller effective diffusion coefficients as the size ratio increases, probably due to caging of the small passive particles in between the big active ones. We expect that under confinement, the caging effect on the passive particles would be amplified, so hindering their transport properties and possibly negatively affecting the sorting.

Finally, we could also address the effects of inertia in the sorting kinetics to extrapolate and generalize our results. Our preliminary numerical results suggest that inertia improves sorting in steady state and decreases the time scale by approximately one order of magnitude (see Fig. S12, SI).

Author contributions

Horacio Serna: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, validation, visualization, writing – original draft, writing – review & editing. C. Miguel Barriuso G.: conceptualization, data curation, formal analysis, investigation, methodology, software, validation, visualization, writing – review & editing. Ignacio Pagonabarraga: conceptualization, formal analysis, investigation, funding acquisition, writing – review & editing. Marco Polin: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, supervision, validation, writing – review & editing. Chantal Valeriani: conceptualization, formal analysis, funding acquisition, investigation, methodology, project administration, resources, supervision, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm00806a.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgements

M. P. acknowledges funding from MINECO through projects CNS2022-135975 and PID2023-146578NB-I00, and acknowledges that the fact that IMEDEA is an accredited ‘María de Maeztu Excellence Unit’ (grant CEX2021-001198, funded by MCIN/AEI/10.13039/501100011033). C. V. acknowledges fundings IHRC22/00002 and PID2022-140407NB-C21 funded by MCIN/AEI/10.13039/501100011033 and by FEDER, UE. H. S. and C. V. acknowledge the received funding from the European Union's Horizon research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 101108868 (BIOMICAR). H. S. acknowledges enriching discussions with José Martín-Roca, Juan Pablo Miranda López and Sujeet Kumar Choudhary. I. P. acknowledges support from Ministerio de Ciencia, Innovación y Universidades MCIU/AEI/FEDER for financial support under grant agreement PID2021-126570NB-100 AEI/FEDER-EU, and Generalitat de Catalunya for financial support under Program Icrea Acadèmia and project 2021SGR-673.

Notes and references

1. N. Heddergott, T. Krüger, S. B. Babu, A. Wei, E. Stellamanns, S. Uppaluri, T. Pfohl, H. Stark and M. Engstler, PLoS Pathog., 2012, 8, e1003023.
2. D. J. Montagnes, A. B. Barbosa, J. Boenigk, K. Davidson, K. Jürgens, M. Macek, J. D. Parry, E. C. Roberts and K. Simek, Aquat. Microb. Ecol., 2008, 53, 83–98.
3. G. Noselli, A. Beran, M. Arroyo and A. DeSimone, Nat. Phys., 2019, 15, 496–502.
4. X.-L. Wu and A. Libchaber, Phys. Rev. Lett., 2000, 84, 3017.
5. C. Valeriani, M. Li, J. Novosel, J. Arlt and D. Marenduzzo, Soft Matter, 2011, 7, 5228–5238.
6. A. Jepson, V. A. Martinez, J. Schwarz-Linek, A. Morozov and W. C. Poon, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 041002.
7. A. E. Patteson, A. Gopinath, P. K. Purohit and P. E. Arratia, Soft Matter, 2016, 12, 2365–2372.
8. R. Jeanneret, D. O. Pushkin, V. Kantsler and M. Polin, Nat. Commun., 2016, 7, 12518.
9. K. C. Leptos, J. S. Guasto, J. P. Gollub, A. I. Pesci and R. E. Goldstein, Phys. Rev. Lett., 2009, 103, 198103.
10. H. Kurtuldu, J. S. Guasto, K. A. Johnson and J. P. Gollub, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 10391–10395.
11. F. Hauke, H. Löwen and B. Liebchen, J. Chem. Phys., 2020, 152, 014903.
12. J. Stürmer, M. Seyrich and H. Stark, J. Chem. Phys., 2019, 150, 214901.
13. I. P. Madden, L. Wang, J. Simmchen and E. Luijten, Small, 2022, 18, 2107023.
14. K. Singh, A. Yadav, P. Dwivedi and R. Mangal, Langmuir, 2022, 38, 2686–2698.
15. R. Wittkowski, J. Stenhammar and M. E. Cates, New J. Phys., 2017, 19, 105003.
16. S. C. Takatori and J. F. Brady, Soft Matter, 2015, 11, 7920–7931.
17. S. Ye, P. Liu, F. Ye, K. Chen and M. Yang, Soft Matter, 2020, 16, 4655–4660.
18. T. Dhar and D. Saintillan, Sci. Rep., 2024, 14, 11844.
19. J. Stenhammar, R. Wittkowski, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2015, 114, 018301.
20. J. Zhang, T. Huang, G. Xu and Y. Chen, Commun. Theor. Phys., 2022, 74, 075601.
21. V. Semwal, A. Kumar, J. P. Singh and S. Mishra, Eur. Phys. J.: Spec. Top., 2024, 233, 3185–3192.
22. C. Maggi, M. Paoluzzi, L. Angelani and R. Di Leonardo, Sci. Rep., 2017, 7, 17588.
23. J. Shea, G. Jung and F. Schmid, Soft Matter, 2022, 18, 6965–6973.
24. D. R. Rodriguez, F. Alarcon, R. Martinez, J. Ramrez and C. Valeriani, Soft Matter, 2020, 16, 1162–1169.
25. P. Dolai, A. Simha and S. Mishra, Soft Matter, 2018, 14, 6137–6145.
26. Y. Wang, Z. Shen, Y. Xia, G. Feng and W. Tian, Chin. Phys. B, 2020, 29, 053103.
27. L. Angelani, C. Maggi, M. Bernardini, A. Rizzo and R. Di Leonardo, Phys. Rev. Lett., 2011, 107, 138302.
28. H. Seyforth, M. Gomez, W. B. Rogers, J. L. Ross and W. W. Ahmed, Phys. Rev. Res., 2022, 4, 023043.
29. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
30. G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp and L. Alvarez, et al., J. Phys.: Condens. Matter, 2020, 32, 193001.
31. G. Gompper, H. A. Stone, C. Kurzthaler, D. Saintillan, F. Peruani, D. A. Fedosov, T. Auth, C. Cottin-Bizonne, C. Ybert, E. Clément, T. Darnige, A. Lindner, R. E. Goldstein, B. Liebchen, J. Binysh, A. Souslov, L. Isa, R. di Leonardo, G. Frangipane, H. Gu, B. J. Nelson, F. Brauns, M. C. Marchetti, F. Cichos, V.-L. Heuthe, C. Bechinger, A. Korman, O. Feinerman, A. Cavagna, I. Giardina, H. Jeckel and K. Drescher, J. Phys.: Condens. Matter, 2025, 37, 143501.
32. M. S. Aporvari, M. Utkur, E. U. Saritas, G. Volpe and J. Stenhammar, Soft Matter, 2020, 16, 5609–5614.
33. S. M. Mousavi, G. Gompper and R. G. Winkler, J. Chem. Phys., 2021, 155, 044902.
34. S. Mallory, C. Valeriani and A. Cacciuto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 032309.
35. R. Ni, M. A. Cohen Stuart and P. G. Bolhuis, Phys. Rev. Lett., 2015, 114, 018302.
36. A. Kaiser, A. Peshkov, A. Sokolov, B. Ten Hagen, H. Löwen and I. S. Aranson, Phys. Rev. Lett., 2014, 112, 158101.
37. C. Wang and H. Jiang, Soft Matter, 2020, 16, 4422–4430.
38. R. Di Leonardo, L. Angelani, D. Dell’Arciprete, G. Ruocco, V. Iebba, S. Schippa, M. P. Conte, F. Mecarini, F. De Angelis and E. Di Fabrizio, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 9541–9545.
39. A. Sokolov, M. M. Apodaca, B. A. Grzybowski and I. S. Aranson, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 969–974.
40. P. Pietzonka, É. Fodor, C. Lohrmann, M. E. Cates and U. Seifert, Phys. Rev. X, 2019, 9, 041032.
41. É. Fodor and M. E. Cates, Europhys. Lett., 2021, 134, 10003.
42. F. Schmidt, B. Liebchen, H. Löwen and G. Volpe, J. Chem. Phys., 2019, 150, 094905.
43. I. Buttinoni, G. Volpe, F. Kümmel, G. Volpe and C. Bechinger, J. Phys.: Condens. Matter, 2012, 24, 284129.
44. J. Arrieta, M. Polin, R. Saleta-Piersanti and I. Tuval, Phys. Rev. Lett., 2019, 123, 158101.
45. X. Garcia, S. Rafa and P. Peyla, Phys. Rev. Lett., 2013, 110, 138106.
46. S. Yang, M. Huang, Y. Zhao and H. Zhang, Phys. Rev. Lett., 2021, 126, 058001.
47. S. Williams, R. Jeanneret, I. Tuval and M. Polin, Nat. Commun., 2022, 13, 4776.
48. J. Martín-Roca, C. M. Barriuso G, R. Martínez Fernández, C. Betterelli Giuliano, R. Zhang, C. Valeriani and L. G. Wilson, Proc. Natl. Acad. Sci. U. S. A., 2025, 122, e2409510121.
49. P. Galajda, J. Keymer, P. Chaikin and R. Austin, J. Bacteriol., 2007, 189, 8704–8707.
50. R. Martinez, F. Alarcon, J. L. Aragones and C. Valeriani, Soft Matter, 2020, 16, 4739–4745.
51. A. Kaiser, H. Wensink and H. Löwen, Phys. Rev. Lett., 2012, 108, 268307.
52. N. Kumar, R. K. Gupta, H. Soni, S. Ramaswamy and A. Sood, Phys. Rev. E, 2019, 99, 032605.
53. O. Chepizhko and F. Peruani, Phys. Rev. Lett., 2013, 111, 160604.
54. M. Mijalkov and G. Volpe, Soft Matter, 2013, 9, 6376–6381.
55. O. Yasa, P. Erkoc, Y. Alapan and M. Sitti, Adv. Mater., 2018, 30, 1804130.
56. F. Martinez-Pedrero and P. Tierno, Phys. Rev. Appl., 2015, 3, 051003.
57. Y. Yang and M. A. Bevan, Sci. Adv., 2020, 6, eaay7679.
58. Y. Fu, H. Yu, X. Zhang, P. Malgaretti, V. Kishore and W. Wang, Micromachines, 2022, 13, 295.
59. T. Ostapenko, F. J. Schwarzendahl, T. J. Böddeker, C. T. Kreis, J. Cammann, M. G. Mazza and O. Bäumchen, Phys. Rev. Lett., 2018, 120, 068002.
60. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2009, 18, 015012.
61. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
62. H. A. Lorentz, Ann. Phys., 1881, 248, 127–136.
63. D. Berthelot, Compt. Rendus, 1898, 126, 15.
64. A. R. Sprenger, L. Caprini, H. Löwen and R. Wittmann, J. Phys.: Condens. Matter, 2023, 35, 305101.
65. A. Deblais, T. Barois, T. Guerin, P.-H. Delville, R. Vaudaine, J. S. Lintuvuori, J.-F. Boudet, J.-C. Baret and H. Kellay, Phys. Rev. Lett., 2018, 120, 188002.
66. H. C. Berg, Sensory perception and transduction in aneural organisms, Springer, 1985, pp. 19–30.
67. C. Kurzthaler, Y. Zhao, N. Zhou, J. Schwarz-Linek, C. Devailly, J. Arlt, J.-D. Huang, W. C. Poon, T. Franosch and J. Tailleur, et al., Phys. Rev. Lett., 2024, 132, 038302.
68. G. W. Ordal, H. M. Parker and J. R. Kirby, J. Bacteriol., 1985, 164, 802–810.
69. L. Turner, L. Ping, M. Neubauer and H. C. Berg, Biophys. J., 2016, 111, 630–639.
70. S. Johnson, B. Freedman and J. X. Tang, Phys. Rev. E, 2024, 109, 064402.
71. M. Polin, I. Tuval, K. Drescher, J. P. Gollub and R. E. Goldstein, Science, 2009, 325, 487–490.
72. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
73. R. J. LeVeque, Finite volume methods for hyperbolic problems, Cambridge University Press, 2002, vol. 31.
74. K. W. Morton and D. F. Mayers, Numerical solution of partial differential equations: an introduction, Cambridge University Press, 2005.

---