Confinement-induced accumulation and de-mixing of microscopic active-passive mixtures

Understanding the out-of-equilibrium properties of noisy microscale systems and the extent to which they can be modulated externally, is a crucial scientific and technological challenge. It holds the promise to unlock disruptive new technologies ranging from targeted delivery of chemicals within the body to directed assembly of new materials. Here we focus on how active matter can be harnessed to transport passive microscopic systems in a statistically predictable way. Using a minimal active-passive system of weakly Brownian particles and swimming microalgae, we show that spatial confinement leads to a complex non-monotonic steady-state distribution of colloids, with a pronounced peak at the boundary. The particles’ emergent active dynamics is well captured by a space-dependent Poisson process resulting from the space-dependent motion of the algae. Based on our findings, we then realise experimentally the de-mixing of the active-passive suspension, opening the way for manipulating colloidal objects via controlled activity fields.

Evolution has enabled living systems to achieve an exquisite control of matter at the microscopic level. From the precise positioning of chromosomes along the mitotic spindle 1 to the many types of embryonic gastrulation 2 cells harness their internal and external motility to reach a predictable order despite the stochasticity intrinsic to the microscopic realm 3 . Understanding how order emerges in these active systems is a fundamental scientific challenge with the potential to bring disruptive technologies for the macroscopic control of microscopic structures. Here we address this problem within a minimal active-passive experimental model system.
From a physical perspective, living systems fall under the general category of active matter 4 , characterised by emergent phenomena including flocking 5-7 , active turbulence [8][9][10][11] and motility-induced phase separation 12,13 . As our understanding of single-species active systems progresses, attention has started to veer towards more complex cases where components with different levels of activity interact. Indeed this is often the case for biological active matter. Intracellular activity, for example, can be used for spatial organisation of passive intracellular organelles [14][15][16] ; and clustering induced by motility differentials helps bacterial swarms expand despite antibiotic exposure 17 . In order to study the emergent properties of these complex systems, an important and phenomenologically rich minimal model is one that mixes active and inert agents. Active impurities have been used to alter the dynamics of grain boundaries in colloidal crystals 18,19 and favour the formation of metastable clusters in semi-dilute suspensions 20 . Sufficiently large concentrations of both active and passive species often reveal a rich spectrum of phases depending on the interactions between constituents [21][22][23][24][25][26][27][28][29] .
Active baths, where individual passive inclusions are dispersed within an active suspension, are particularly appealing. Their conceptual simplicity makes them a natural starting point to develop a statistical theory of active transport [30][31][32][33] with potential applications to micro-cargo delivery, micro-actuation [34][35][36] and nutrient transport 37,38 . In general, passive particles in homogeneous and isotropic active baths display enhanced diffusion due to a continuous energy transfer from the active component via direct collisions and hydrodynamic interactions [39][40][41][42][43][44][45][46][47] . Designing larger objects with asymmetric shapes can then turn active diffusion into noisy active translation or rotation 30,48,49 .
However the strategies for the predictable patterning and transport of passive cargo are still very limited. This is in contrast with the exquisite level of external control possible on the motility of biological and synthetic active particles themselves [50][51][52][53][54][55] , such as photokinetic bacteria whose swim speed depends on the local light intensity; property that can be harnessed to pattern surfaces 53,54 .
Here we show that a steady gradient of activity within activepassive suspensions can be harnessed to control the fate of generic passive particles. Quasy-2D microfluidic channels are filled with a binary suspension of polystyrene colloids and unicellular biflagellate microalgae Chlamydomonas reinhardtii (see Supplementary Movie 1). Confinement induces a spatially inhomogeneous and anisotropic distribution of microswimmers as a result of wall scattering [56][57][58][59] , which translates into a space-dependent active noise for the colloids. We show that this produces complex non-monotonic colloidal distributions with accumulation at the boundaries and then develop an effective microscopic jump-diffusion model for the colloidal dynamics. The latter shows excellent analytical and numerical agreement with the experimental results. Finally we demonstrate how confinement with aptly designed microfluidic chips can fuel the de-mixing of activepassive suspensions, opening the way for manipulating passive colloidal objects via controlled activity fields.

Colloid concentration is non-homogeneous under confinement
A full description of the experimental procedures (culturing, video microscopy and data analysis) can be found in the Methods section and Supplementary Materials (data and codes also available on Zenodo 60 ). Figure 1a shows a detail of the main section of the first type of microfluidic setup used. These devices, between 14 μm and 20 μm thick, are composed of two reservoir chambers connected by long and straight channels of constant width 2W ranging from 50 to 200 μm. These are filled with a dilute mixture of Chlamydomonas reinhardtii (CR; strain CC-125, radius R~4-5 μm) and weakly Brownian colloids (radius a = 5 ± 0.5 μm) at surface fractions ϕ CR ≈ ϕ col ≈ 2-3% (bulk concentrations~2−3 × 10 7 particles/mL). The design ensures a steady concentration along the connecting channels for both species (see Supplementary Movie 1). Figure 1a shows the coordinate system employed, which has been symmetrized with respect to the channel axis.
As typical of self-propelled particles, the algae tend to accumulate at the boundaries due to their interactions with the side walls [56][57][58] . This appears as a significant peak in their time-averaged concentration profile at a position y CR ≈ 15 μm ( Supplementary Fig. 1), roughly equivalent to the sum of the cell radius and the flagellar length. Unexpectedly, we find that also the colloids explore the available space in an inhomogeneous way. Their steady-state distribution (Fig. 1b) shows a clear peak at about one particle radius (y col = 5.9 ± 1.5 μm), followed by a depleted region between 10 μm and 20 μm from the wall, which roughly corresponds to the peak in algal density, before plateauing to a uniform concentration further inside the channels. This effect is also observed within circular chambers ( Supplementary Fig. 2) suggesting it is a robust feature of the system. These distributions are in stark contrast with the equilibrium case (i.e., without microswimmers) for which the colloids are expected to be uniformly distributed despite spatial variations in colloidal diffusivity due to hydrodynamics 61 (Supplementary Fig. 3).
As a first step to gain insight into the colloids' experimental distributions, we characterise their dynamics within our active bath. Figure 1c shows a typical colloidal trajectory (561 s) colour-coded for the average frame-to-frame speed v. The dynamics can be understood as a combination of periods of slow diffusive-like displacements (v~5 μm/s;) and fast longer jumps (v~30 μm/s; see Supplementary Section 1). The latter is reminiscent of hydrodynamic entrainment events reported for micron-sized colloids 47,62 , although for these larger particles the prominent role appears to be played by direct flagellar interactions (see Supplementary Movie 2). Irrespective of their specific origin, jump events dominate the active transport of colloids within the system. Following 47,63 , we use time-correlation of displacements along individual trajectories and an estimate for the expected magnitude of frame-to-frame diffusive displacements to identify active colloidal jumps and extract their statistical properties as a function of starting position y (we assume translational invariance along the channels' axis). Figure 2a shows the space-dependent distributions of waiting times between successive jumps, with colours representing the distance from the boundary. Regardless of the position, all the distributions are exponential above~3 s but deviate from it at shorter times. This deviation from a simple Poisson process is due to the large colloidal size influencing the motility of nearby algae as reported also in the case of bacteria 64 . Here, this can lead to a rapid succession of interactions with the same cell (see Supplementary Movie 3). Nevertheless, for our purposes, the waiting dynamics can still be approximated with a single effective rate λ(y) (  Figure 2b shows that these rates increase monotonically with increasing distance from the wall. Notice that this curve does not mirror the algal accumulation at the boundary ( Supplementary Fig. 1), a reflection of the fact that proximity to the wall curtails the range of possible swimming directions that algae can have when interacting with the colloid.
Next, we look at the modulation in jump orientation and magnitude. Figure 2c shows the distribution of jump directions P(θ, y), with motion towards or away from the boundary corresponding to negative and positive values of θ respectively. Approaching the wall, the distribution develops a marked peak at θ = 0 indicating a strongly anisotropic active motion preferentially parallel to the boundary. This feature reflects the anisotropy in algal dynamics that results from the interaction with the wall 56 and that can be measured up tõ 100 μm from the boundary because of the persistence in cells' trajectories 59 . The distribution functions of jump magnitudes are similar to those already reported for micron-sized particles 47,62 , with an exponential decay above~4 μm ( Supplementary Fig. 4a). These can be used to calculate the average jump length 〈l(y)〉 (Supplementary Fig. 4b) which decreases by~30% from the bulk level within The two species are highlighted with, respectively, green ellipses and red circles. The schematic illustrates the coordinate system used throughout. b Steady-state colloidal probability distribution functions across the channel (edge to midpoint) for different W values. Boundary accumulation is followed by a depleted adjacent region which gives way to a uniform distribution within the channel. c Example of a typical colloidal trajectory inside a 100 μm-width channel. Colours represent the mean frame-to-frame speed. The dynamics is composed of slow diffusive-like motion (blue sections) interspersed with fast straighter jumps (red sections). the first 10−15 μm from the boundary (first 5−10 μm accessible to the beads). A decrease is indeed expected due to the obvious limitations to the active movement of the colloids imposed by a nearby boundary.
As we are interested in the passive particles' dynamics across the channels, it is useful to reduce the full two-dimensional jump distribution functions to their projection q y (y J ) along the y axis. Here y J is the y-coordinate of the active displacement and the subscript indicates the distance from the nearest wall (see Fig. 1a for the frame of reference). Figure 2d shows that these distributions are generally well approximated by the combination of two exponentials, one each for positive and negative directions (respectively away from and towards the boundary). The exception is for the positive tails of the distributions closest to the wall, which decay slower than expected from the exponential fit. They will not be considered in the following. The distributions q y (y J ) can then be approximated as where L ± (y) are the characteristic lengths of the exponential fits to positive and negative jumps respectively (see Methods for details on how the characteristic lengths were extracted). As shown in Fig. 3-inset (red and blue symbols and solid lines) these are identical in the core of the channel (L + = L − ≃ 4 μm) and decrease to the same small value close to the wall (~0.5 μm) as a consequence of the increasing polarisation of the active displacements along the boundary (Fig. 2c). However, the length scale for the transition between boundary and core values is shorter for L + (y) than for L − (y) (~12 μm vs.~25 μm, see Fig. 3-inset). This behaviour, common to all values of W explored ( Supplementary Fig. 5), gives rise to a net drift L + (y)−L − (y) > 0 towards the bulk of the system (green symbols and solid line, Fig. 3-inset) which, as discussed below, is responsible for the depleted region observed in the experimental colloidal distributions. The difference between the L's is likely due to a combination of local anisotropy in the distribution of swimming directions of the algae, and the stopping of colloids' active displacements by the wall. The focus on the active jumps displayed by the passive particles can only be justified if this part of the dynamics is indeed sufficient to capture the experimental steady-state distributions. This was first tested within a 1D numerical simulation of a weakly Brownian particle (diffusivity D 0 ), moving in y ∈ [0, 2W] and subject to a spacedependent Poisson noise of rate λ(y) and value drawn from q y (y J ) (see Methods Sec. 4.5 for details on the integration scheme). Figure 3 shows that the resulting steady-state distribution (black solid curve) is in excellent agreement with the experimental one (black circles) (2W = 100 μm). Numerical simulations also provide a convenient way to explore which elements of the effective colloidal dynamics play the most prominent role. For example, fixing λ(y) to the bulk value everywhere leads to a very minimal change in the spatial distribution (Fig. 3, teal dashed line), showing that the spatial dependence of the jump frequency close to the wall is not a major factor in the present case. Similarly, removing completely the background diffusion by setting D 0 = 0 leaves the distribution unchanged (Fig. 3, purple dashed line) consistent with the fact that, in our experimental system, colloidal transport is dominated by the Poissonian jumps. On the other hand, simulations that remove the spatial dependence of the jump distributions q y (y J ) and use their isotropic bulk value everywhere, show a boundary accumulation of colloids that is narrower and higher than the experimental curve and has no intermediate depletion (Fig. 3, orange solid line). This confirms that space-dependent jump anisotropy plays a key role in determining the experimental colloidal distribution.

Minimal continuum model for the dynamics of the colloids
The numerical validation of the jump-diffusion dynamics motivates a simple analytical model for the evolution of P t (y), the probability  density of finding a colloid at position y at time t: λðy À y J ÞP t ðy À y J Þq yÀy J ðy J Þdy J : In this master equation, local changes in P t (y) are due either to diffusion (first term in the r.h.s) or to the balance between active jumps from the current position or towards it from elsewhere (second and third terms in the r.h.s., respectively). This continuum equation can also be derived more formally from a stochastic description of the dynamics of single colloids, following Denisov and Bystrik 65 (see Supplementary Section 4). It is worth noting also that Eq. (2) assumes infinitely fast jumps ("teleportation"). While this might seem like a drastic approximation, it is not expected to impact the resulting spatial distributions as long as the typical jump duration is sufficiently shorter than the average waiting time between jumps. This is indeed the case in the current system (~1.2 s vs. ≳3.5 s respectively; see Fig. 2b and Supplementary Fig. 6).
The non-local term in Eq.
(2) means that the model is nontractable for generic kernels λ(y)q y (y−y J ). In order to make progress, we, therefore, perform a Kramers-Moyal expansion 66 and truncate the series at second-order, providing an approximation of the system at the drift-diffusion level. This type of approximation has already been used successfully to capture the essential features of run-and-tumble bacteria 67 , active Brownian particles 68 and active Ornstein-Uhlenbeck particles 69 . This approximation is expected to hold if the characteristic jump size is small compared to the length scale of the heterogeneities in the dynamics. In our case we have 〈l bulk 〉 ≈ 5 μm (Supplementary Figure 4b) while a length scale for the heterogeneity can be estimated for example from the L ± curves (≳ 12 μm).
Notice that D eff is similar to the asymptotic diffusivity obtained in 47 for homogeneous and isotropic entrainment-dominated transport, with a contribution from the jump events given by the product of jump frequency and variance of jump size. Enforcing no-flux at the boundaries one obtains a closed form for the steady-state solution: where B is the normalisation constant. Without a net asymmetry in the dynamics (i.e., when m 1 (y) ≡ 0), the solution is that of a statedependent diffusive process with Itô convention for multiplicative noise integration. In our case, however, the first moment of q y takes significant values in the first~30 μm from the boundaries (green solid line Fig. 3-Inset) and should not be neglected. Comparison between Eq. (4) and the experiments is facilitated by having analytical expressions for λ, m 1 and m 2 . Figure 2b shows that λ can be described heuristically by a simple exponential relaxation from the wall to the bulk values (black solid line). As for m 1 and m 2 , the description of q y (y J ) given by Eq. (1) implies that m 1 ðyÞ = L + ðyÞ À L À ðyÞ Given that L ± (y) themselves are well approximated by exponentials, this affords analytical approximations also for m 1,2 (y) (Fig. 2c-Inset). Armed with this description we can now compare the experimental curves with the prediction from the approximate jump-diffusion model. The olive-coloured curve in Fig. 3 shows that the theoretical distributions recapitulate extremely well the experimental ones. On one hand, this provides a justification a posteriori for the modelling approach taken above. On the other it reveals that, at a continuum level, the dynamics of colloids in an active bath can be reduced to two quantities: the first and the second moments of the active displacements. This provides an intuitive way to conceptualise the complex dynamics of colloidal particles within an active suspension.

Induced de-mixing of the passive particles
What we learned so far on the colloids' emergent active dynamics can be harnessed to induce the de-mixing of our active-passive suspensions, as a first step towards more complex control of microscale cargo transport. The idea is to include a confining boundary which acts as a kind of selective membrane, letting the beads cross but not the algae. This is achieved by considering circular chambers decorated with 190 μm-long side channels whose width d (Fig. 4a) is smaller than the swimmer diameter but larger than the bead diameter (here taken as 2a = 6 μm < d = 7−8 μm < 2R = 8−10 μm). Due to the active dynamics of the colloids, we expect them to be pushed inside the side channels by the algae, therefore depleting the particles in the main chamber and spatially separating the active and passive species. As shown Fig. 4b (see also Supplementary Movie 4) this is exactly what we obtain, with the number of particles in the circular chamber N c (red solid line) quickly decreasing and the number in the side channels N s (cyan solid line) increasing. Notice that the total number of particles N t = N c + N s in the system (purple solid line) remains constant over the duration of the experiment. The process follows a first order kinetic with rates k in and k out for transitions towards or out-of the side channels respectively (see Supplementary Section 6). The occupancy numbers follow an exponential relaxation with time-scale 1/(k in + k out ) = 77 ± 4 min and stationary values N c =N t = k out =ðk in + k out Þ = 0:64 ± 0:02 and N s =N t = k in =ðk in + k out Þ = 0:36 ± 0:02, leading to estimates for the two rates of k in = (7.8 ± 0.6) × 10 −5 s −1 and k out = (13.9 ± 0.8) × 10 −5 s −1 . The latter results from the passive dynamics of the colloids within the side channels, which depends on their thermal diffusivity, potential electrostatic interactions with the boundaries, short-range hydrodynamic interactions between the beads 70 , channel design, etc 71 . The former, instead, is the consequence of active colloidal displacements. Assuming no colloid/wall attraction within the side channels (as confirmed by our control experiment, see below), it encapsulates the intrinsic outof-equilibrium nature of the system, a property which is immediately clear from the significant difference in steady-state densities between the side channels (ρ s = 0.0295 ± 0.007 beads/μm 2 ) and in the main chamber (ρ c = 0.0019 ± 0.0001 beads/μm 2 ) (Fig. 4c). The value of k in can be compared with an estimate derived from the effective colloidal dynamics from Eq. (3) (see Supplementary Section 7). The mean first passage time of a colloid to the entrance of any of the side channels returns a rate of k V ,D in = ð7:408 ± 0:1Þ × 10 À5 s −1 which agrees very well with the experimental value. This shows that the simplified advectiondiffusion model for the passive particles is a good description not just for time-independent properties like the steady-state distribution of Fig. 3, but also for dynamic ones like the rate of active cargo transport to the side channels. It is instructive also to use the model to explore the contributions of the different parts of the dynamics. Hence we see that setting V eff ≡ 0 and D eff ≡ 3.55 μm 2 s −1 , the effective diffusivity far from the boundaries, one obtains a rate k D in = ð16:2 ± 4:1Þ × 10 À5 s −1 . The stark difference between k V ,D in and k D in suggests that, in our case, it is not sufficient to know the behaviour of the active-passive suspension in the bulk of the chamber. The behaviour close to the boundary, where the inhomogeneities are found, is essential to grasp the dynamics of the de-mixing process. Repeating the same experiment only with colloids shows indeed no appreciable filling of the side channels within the duration of the experiment (~10 h; see Supplementary Movie 5) and the particle density in the main chamber remaining constant (Fig. 4b  red dashed line).
Although the microfluidic chip in Fig. 4 leads to a much higher concentration of colloids in the side channels than in the main chamber, this still involves only~36% of the colloids. A different design, however, could maximise the fraction of de-mixed colloids instead of their concentration. For example, modifying the current design to have two identical chambers connected by a narrow gap of size d = 7−8 μm should produce a fraction of de-mixed colloids given by The design could then be tailored to different needs.

Discussion
In this study we used an active-passive system as a minimal model to investigate the active transport of passive microscopic objects. Leveraging the well-known boundary accumulation of microswimmers, we saw that the emerging interactions between active and passive components are also heterogeneous and anisotropic. In turn, they lead to a complex spatial distribution of colloids and to the possibility of demixing the suspension. These properties are well described by a simple advection-diffusion model for the passive species, based only on the first and second moments of their ensuing active displacements. This type of approximation, then, provides a degree of universality in the description of the dynamics of active systems including purely active and active-passive species. From the targeted delivery of chemicals, biomarkers or contrast agents at specific locations in the body to the remediation of water and soils, artificial and biological microscopic selfpropelled particles hold a great potential to address critical issues in areas ranging from personalised medical care to environmental sustainability 72,73 . Development in these areas will depend on our ability to use and control micro/nano-swimmers to perform essential basic functions, such as sensing, collecting and delivering passive cargoes in an autonomous, targeted and selective way 74 . Within this perspective, we have shown here that heterogeneous activity fields can be employed to control the fate of colloidal objects in a manner similar to biological systems, where the spatial modulation of the activity can control, for example, the positioning of organelles within cells [14][15][16] . Although this was achieved here through simple confinement, many microswimmers possess -or can be designed to have-the ability to respond to a range of external physico-chemical cues or stimuli [51][52][53][54] . We expect that similar phenomenologies should also emerge whenever the active particles in a mixed system autonomously generate a space-dependent concentration profile. Our work paves the way to using these to dynamically manipulate the fate of colloidal cargoes by externally altering microswimmers' dynamics. Advance in this area will require to understand how to predict the active displacements of passive cargoes directly from the dynamics of microswimmers 75 , a theoretical development that we leave for a future study.

Microfluidics and microscopy
The mixed colloid-swimmer solution was injected into polydimethylsiloxane-based microfluidic channels, previously passivated with a 0.5% w/w Pluronic-127 solution to prevent sticking of the particles to the surfaces of the chips. The microfluidic devices were then sealed with Vaseline to prevent evaporation. The channels of varying widths were observed under bright-field illumination on a Zeiss AxioVert 200M inverted microscope. A long-pass filter (cutoff wavelength 765 nm) was added to the optical path to prevent phototactic responses of the cells. For each channel width stacks of 80,000 images at ×10 magnification were acquired at 10 fps with an IDS UI-3370CP-NIR-GL camera.

Data analysis: trajectories
Particle Trajectories were digitised using a standard Matlab particle tracking algorithm (The code can be downloaded at http://people. umass.edu/kilfoil/downloads.html). After the particle trajectories were reconstructed, the large jumps were isolated using a combination of move size and directional correlation along the trajectory (see detailed protocol in 47,63 ).
The experimental curves L + (y) and L − (y) were obtained from the experimental jump size distributions q y (y J ) by fitting the logarithm of each side (i.e., positive and negative y J ) with affine functions. Because these distributions deviate from simple exponentials as we get closer to the boundary (Fig. 2d), the fits were restricted to |y J |< limðyÞ in order to extract the characteristic length scales of the initial decay of the distributions. The resulting fits are shown in Fig. 2d (solid lines). The limits for the fits of the experimental distributions were fixed as in Table 2.

Numerical simulations
One hundred trajectories composed of 10,000 time-steps (with δt = 0.1 s, equal to the acquisition period of the experiments) were simulated for each condition. At each time step an acceptance-rejection scheme is performed to decide whether the particle undergoes a jump: (i) a (uniform) random number is first drawn from the interval [0,1], (ii) if this number sits within the interval [(1−λ(y)δt)/2; (1 + λ(y)δt)/2], the particle undergoes a jump whose size is taken from the experimental distribution q y (y J ) (using inverse transform sampling), (iii) otherwise the particle simply undergoes Brownian motion with diffusivity D 0 = 0.0439 μm 2 s −1 (the theoretical bulk diffusivity for a 5 μm-radius particle at room temperature). Upon reaching a boundary jump moves undergo a stopping boundary condition, stopping at ±5 μm (the particle's radius) from the boundary (this is to be as faithful with our experimental observations as possible). For simplicity, diffusive moves that cross the boundary undergo a reflective boundary condition. More accurate boundary implementations for the diffusive motion do not lead to different results as the impact of thermal diffusivity in the channel experiments is negligible.

Data availability
The main data used in this study are available in the Zenodo database with https://doi.org/10.5281/zenodo.6866462.

Code availability
The codes used in this study are available in the Zenodo database with https://doi.org/10.5281/zenodo.6866462.