Lock-key microfluidics: simulating nematic colloid advection along wavy-walled channels

Abstract

Liquid crystalline media mediate interactions between suspended particles and confining geometries, which not only has potential to guide patterning and bottom-up colloidal assembly, but can also control colloidal migration in microfluidic devices. However, simulating such dynamics is challenging because nemato-elasticity, diffusivity and hydrodynamic interactions must all be accounted for within complex boundaries. We model the advection of colloids dispersed in flowing and fluctuating nematic fluids confined within 2D wavy channels. A lock-key mechanism between homeotropic colloids and troughs is found to be stronger for planar anchoring on the wavy walls compared to homeotropic anchoring on the wavy walls due to the relative location of the colloid-associated defects. Sufficiently large amplitudes result in stick–slip trajectories and even permanent locking of colloids in place. These results demonstrate that wavy walls not only have potential to direct colloids to specific docking sites but also to control site-specific resting duration and intermittent elution.

---

1 Introduction

Manipulating particle trajectories is an important feature of many microfluidic applications, including drug delivery,¹ cell sorting² and medical diagnoses.³ Accordingly, a variety of approaches have been developed to control the trajectories of microparticles.^4–7 One means of sculpting particle trajectories is by suspending them in complex carrier fluids, as in topological microfluidics,⁸ in which particles are guided by liquid crystalline carrier fluids.⁹ Liquid crystalline materials are alluring for microfluidic transport because they are highly responsive to flows,^10–12 suspended inclusions^13,14 and confining surfaces.^15–17 When colloidal particles are dispersed in liquid crystalline fluids, the anisotropic nature of liquid crystals gives rise to emergent properties^18–21 and imposed anchoring at colloidal surfaces results in topological defects in the vicinity of the colloids to ensure topological charge neutrality. Strong homeotropic or planar anchoring endows the surface with a topological charge of +1 and necessitates the existence of an accompanying −1 charge in the bulk fluid, either as two −1/2 point defects in 2D, or in 3D as defect loops (Saturn rings), −1 point defect (hyperbolic hedgehog), or surface defects (boojum defects).

Such defects mediate particle–particle^22–24 and particle–wall interactions.^25,26 Particle–wall elastic interactions allow micropatterned structures on walls to have significant, long-ranged effects on colloids.^27,28 For example, indentations in walls can govern colloid position: when the width of an indentation is comparable to colloid size there can be a strong attraction,²⁹ whereas this attraction is weak when the width is much greater or smaller than the diameter of the colloid.^27,30,31 These ideas were advanced by investigating periodic wavy boundaries, in which periodic nematic deformations near the wavy walls causes colloids with Saturn ring defects to be attracted to the troughs, and dock via a ‘lock-and-key’ mechanism, while docking location depends on the orientation of the defect for colloids with hedgehog defects.³² Variations of this setup have included imposed twist,³³ saw-tooth or crenellated substrates,³⁴ and ellipsoidal³⁵ or four-armed colloids.³⁶ In addition to setting stable docking sites, structured boundaries can control dynamic trajectories. Trajectories can arise due to diffusion or external fields,³⁷ with docking or continued motion controlled by a balance between nemato-elasticity and driving force. Periodic structures, in particular, allow precise control of colloidal transport.³⁸ In these cases, nematic colloids must respond to both elastic deformations through free energy minimization and advection through drag.

In this manuscript, we model discoidal colloids advecting in fluctuating 2D nematics flowing through microfluidic wavy channels (Fig. 1) and vary the amplitude of the undulations. We investigate this using a multi-particle collision dynamics (MPCD) algorithm to simulate thermalized nematodynamics and the inclusion of colloidal particles, and we provide a method for implementing wavy boundaries. This approach includes colloidal diffusion, nemato-elasticity and hydrodynamics within complex boundaries. Inspired by experimental studies in which docking sites are of comparable size to colloids,^32,33,35 we focus on the effect of the amplitude of the wavy boundaries on the colloidal trajectories and velocities, for different anchoring conditions, and demonstrate that these systems have the potential to trap colloids at given points or allow colloidal conveyance. For intermediate amplitudes, elution follows stick–slip dynamics, with the “sticking” duration increasing with amplitude until the colloid locks in place. The amplitude at which the transition from stick–slip to locking occurs decreases with nematic elasticity for a given pressure gradient. We consider colloids with strong homeotropic anchoring but allow the anchoring on the channel walls to vary. Our results show simple microfluidic systems can temporarily lock particles in place by combining surface structure and advective flow, allowing precise control over elution dynamics and the ability to lock colloids in place for specific times.

Fig. 1 Simulation snapshot, with planar anchoring on the wavy channel walls, and homeotropic anchoring on the mobile colloid. In this study, all colloids have homeotropic anchoring but the anchoring on the channel walls varies. The color shows local scalar order parameter field S and white dashes show the director field. Colloid radius R = 5a; wall amplitude B₀ = 2.5a; wall wavelength λ = 20a; average channel height 〈h〉_x = 20a; and pressure gradient . Snapshot from t = 20τ, early in the simulation.

2 Methods

To simulate colloids in a flowing and thermally fluctuating nematic liquid crystal confined within a complex channel, a multi-particle collision dynamics (MPCD) algorithm is used.^39–42 MPCD has been employed to simulate reaction–diffusion dynamics,^43,44 electrophoresis,^45,46 thermophoresis,⁴⁷ swimmers,^48–50 polymers,^51–54 colloidal suspensions,^55,56 binary mixtures,⁵⁷ viscoelastic fluids,⁵⁸ ferrofluids,^59,60 and dense stellar systems.^61–63 Most relevantly, MPCD has been used to simulate nematic liquid crystals.⁶⁴ In this context, MPCD has simulated nematohydrodynamics,⁶⁵ suspended colloids,^66,67 magnetic colloids,^68,69 living liquid crystals,⁷⁰ and active nematics.^71–73 For the present study, MPCD is chosen because it is ideal for moderate Péclet numbers, mobile solutes and complex boundaries.

The nematic MPCD (N-MPCD) algorithm discretizes the continuous fluid into N point particles. Each such point particle i is given a set mass m_i = m ∀ i, possesses an instantaneous velocity v̲_i(t) and unit-length orientation u̲_i(t) to model the direction of the nematogens at each instant t. The MPCD algorithm consists of two steps:^74,75 (i) streaming and (ii) collision. The streaming step (see Appendix A) updates each particle position r̲_i assuming ballistic motion over a time step δt. The collision step (see Appendix A) simulates the interactions between fluid particles via a coarse-grained collision operator. The collision operator stochastically updates the particles' velocities and orientations, while constrained to respect conservation laws. To do this, the simulation domain is split into cubic cells of size a and index c containing N_c(t) particles at any instant t. A random grid shift ensures Galilean invariance.⁷⁶ Only particles in the same cell interact and the interactions involve all particles in the same cell. The average number of particles per cell sets the fluid density ρ = m〈N_c〉/a^d, where 〈·〉 is the average over all cells and d = 2 for a two-dimensional fluid. The collision event can itself be broken into three stages (see Appendix A):

1. Stochastic momentum exchange: particle velocities are updated by a collision operator Ξ_i,c(t) for particle i in cell c at time t. An Andersen-thermostatted collision operator randomly generates velocities from a Boltzmann distribution with energy k_BT.^77,78 This operator conserves momentum and thermostats the energy.

2. External forces: a pressure gradient is applied as where the effective external acceleration g̲ is included in the collision operator.

3. Orientation exchange: in N-MPCD, each particle's orientation u̲_i is stochastically drawn from the Maier–Saupe distribution about each cell's director n̲_c as p_c(u̲_i) ∝ exp(US_c(n̲_c·u̲_i)²/k_BT), where S_c is the scalar order parameter of cell c and U is the mean field potential, which controls how strongly the orientations align. Additionally, MPCD particle orientation responds to velocity gradients through Jeffery's equation, which is parameterized by a bare tumbling parameter ξ and heuristic shear coupling coefficient χ. The χ parameter acts as a relaxation parameter, effectively allowing Jeffrey's equation to be averaged over a small number of time steps of the fluctuating hydrodynamic field. Finally, director dynamics are coupled back to the velocity field through local torques parameterized by a rotational mobility coefficient γ_R (see Appendix A for details).^64,69 A small value is chosen to keep backflow effects to a minimum. Numerical analysis has shown that N-MPCD describes a linearized nematohydrodynamic model in which viscosity and elastic effects are isotropic.⁷⁹

The wavy 2D channel is modeled with plane-wavy solid boundaries (see Appendix B). Each surface (both bounding walls and colloidal) is labelled by index b and represented as a surface on which . A baseline surface is defined for planes and, to account for wavy boundaries, a second term is included, which allows waves with amplitude B_b,0 and frequency B_b,1. In this study, we focus on 2D; however, this framework allows simulations of egg-carton and corrugated walls (Appendix B). Both the top and bottom wavy walls of the channel have the same amplitude B_b,0 = B₀ and wavelength λ given by for all wall boundaries, where is the system size in that direction. The boundaries defining the discoidal colloids are smooth and so B_b,0 = B_b,1 = 0 for all colloids, though wavy colloids and cylinders are possible in our framework (Appendix C). Bounce-back boundary conditions with phantom particles ensure no-slip.^80–82 Periodic boundary conditions are plane surfaces at the channel extremities.

2.1 Simulation parameters

The simulation units are particle mass m, thermal energy k_BT and cell size a. From these, the time unit is . Additional derived units include dynamic viscosity μ₀ ≡ k_BTτ/a^d,⁷⁸ stress k_BTa^−d and Frank elasticity k_BTa^2−d. In these two dimensional (d = 2) simulations, the time step is set to δt = 0.1τ and density is set to ρ = 20m/a². The mean field potential is U = 10k_BT (unless otherwise stated), placing the fluid deep in the nematic phase.⁶⁴ A mean field potential of U = 5k_BT is relatively close to the nematic–isotropic transition point, with a system-wide (averaged over all cells c) scalar parameter S ≈ 0.5, while U = 10k_BT and 20k_BT have scalar parameters much closer to S ≈ 1.⁶⁴ These large values are an idealization that ensures the system is deep in the nematic phase and that elastic constants are large enough to ensure small Ericksen numbers (Appendix D). The Frank elastic constants were previously found to obey a one constant approximation and be linear with the mean field potential as K = (113 ± 6)U for ρ = 20m/a².⁶⁴ The nematic is in the flow aligning regime with ξ = 2, shear susceptibility χ = 0.5 and rotational friction γ_R = 0.01τk_BT. The choice of small rotational friction minimizes the amount of backflow and limits any backflow-related effects in this study. The gravitational acceleration is set to g̲ = 0.001a/τ²x̂ to generate a pressure gradient . Ten repeats (n = 10) of each configuration are run for 1.5 × 10⁴τ.

The channel is modeled with plane-wavy solid boundaries along x̂. Both walls have the same amplitude B₀ and frequency B₁. Planes set 60a apart define the channel walls with average normal ±ŷ. The domain size is 60a × (20a + 2B₀) and the frequency of the wavy walls, B₁, is set to π, so that 3 full periods fit within the length of the channel, which gives a wavelength λ = 20a. Positions of x = λ/4 and x = 3λ/4 correspond to the center of the trough and crest, respectively. The amplitude of the waves B₀ is varied between 0 and 4.5a.

The colloid is simulated as a moving circular boundary with a radius R = 5a and mass ρπr². The presence of the −1/2 companion defect pair (Fig. 1) indicates that the anchoring is strong (Appendix B) and that the characteristic anchoring length is small compared to the colloid size.⁸³ Colloids are initialized with their center at r̲_c = (x, y) = (25a, 10a + B₀), which is between the wavy walls and offset to the left between a crest and a trough on the walls (Fig. 1), unless otherwise stated. Immediately after initialization, colloids are free to move in response to nematic and hydrodynamic forces. The colloid position r̲_c is a continuous variable that is not constrained to the MPCD collision lattice. These choices define a characteristic length scale L for the geometry of our system. Here, L = 10a coincides with the diameter of the colloid 2R, the average distance from the centerline to the wall 〈h〉/2, and the distance from trough to crest λ/2. This length scale results in characteristic velocities. The velocity scale of the nematoelastic response is Ṽ = K/μL^d−1 for viscosity μ and the scale due to pressure gradients is V = L²ρg_x/μ. For the parameters used here, the nematoelastic speed Ṽ ≈ 12a/τ is much larger than the advective speed V ≈ 0.2a/τ (Appendix D), which we identify to be the characteristic speed of interest in our lock-key dynamics. In this study, all colloids have homeotropic anchoring but the anchoring on the channel walls varies.

3 Results

Lock-key microfluidics result from two competing processes: nematoelastic relaxation pushes the colloids to locations that minimize the distortion free energy, while the minimum dissipation theorem leads colloids to advect with a velocity that matches the fluid flow. For this reason, we begin by simulating homeotropic colloid kinetics in the absence of externally driven flows (Section 3.1), then proceed to consider nematohydrodynamic flows in wavy channels in the absence of the colloid (Section 3.2) before putting the two aspects together to study lock-key microfluidics (Section 3.3). In all cases, each colloid has a topological charge of +1 and so necessitates the existence of companion defects in the vicinity of the colloid. In 3D, both Saturn rings and hedgehog defects are possible⁸⁴ with rings more likely in simulations.^85,86 In 2D only the quadrapolar state is expected to be stable^87,88 and the results reported here are consistent with this since they exhibit a pair of −1/2 companion defects.

3.1 Colloid trajectories with no pressure gradient

We place colloids into channels with one plane boundary wall and one wavy wall, with an amplitude of B₀ = 2.5a. The nematic anchoring on the colloid is always homeotropic but both planar and homeotropic anchoring on the channel walls is considered. No pressure gradient is applied. The colloid is initially positioned either directly above a crest or trough (Fig. 2; top row) for homeotropic anchoring on the boundary walls. The observed numerical trajectories are consistent with Luo et al.³⁷ For homeotropic anchoring on the boundary walls, the colloid is attracted to the trough, traversing away from the crest where it settles (Fig. 2a). The principle motion is in the x̂ direction with primarily stochastic motion in the ŷ direction. Although the colloid does not move substantially in the ŷ direction, careful inspection shows that it does drift upward since the repulsive forces from the different walls are not expected to balance at the same height for each x-position. When the colloid is initiated above the trough (x = λ/4), it still diffuses, but remains within the trough in the vicinity of the free energy minimum, exhibiting no deterministic forcing away from this point (Fig. 2b).

Fig. 2 Typical trajectories for a colloid with homeotropic anchoring confined between a plane wall above and wavy wall below. The start and end points are shown by magenta and cyan dots respectively. (top row) Walls have strong homeotropic anchoring. The colloid is either initiated above a crest (a) or above a trough (b). (bottom row) Walls have strong planar anchoring, while the anchoring on the colloid surface is still strong homeotropic. Colloid radius R = 5a; wall amplitude B₀ = 2.5a; wall wavelength λ = 20a; and average channel height 〈h〉_x = 20a.

Homeotropic anchoring on the colloid and planar anchoring on the walls was not investigated by Luo et al.,³⁷ but shows interesting features (Fig. 2; bottom row). In the case of planar anchoring on the walls, the colloid is stable above the crest (Fig. 2c), rather than the trough. This is consistent with experimental images of the analogous case of a planar colloid being attracted to the crests of a homeotropic wall³⁷ and with the trapping of colloidal particles at the tips of sharp protrusions whenever the colloid and surface have opposing anchoring conditions.⁸⁹ In this case, the colloid simply diffuses as though it is in a well within the free energy landscape. Likewise, a colloid initiated above the trough does not move left or right; but rather, moves directly away from the wavy wall towards the upper plane wall (Fig. 2d). While the wall above the colloid is also repulsive, the repulsion is smaller than from the wavy wall. Thus, the homeotropic colloid's equilibrium position is closer to the plane wall than the wavy wall.

Having considered the dynamics of a colloid confined by one wavy wall, we further consider how homeotropic colloids move in response to wavy boundaries on both the lower and upper wall. The colloid is initialized in the center of the channel, and at the inflection point halfway between a trough and a crest. For both homeotropic and planar anchoring on the walls, colloids are observed to move towards the troughs (Fig. 3) and generally remains centered in the trough for all amplitudes.

Fig. 3 Mean x-position of the colloid with homeotropic anchoring in relation to the trough for planar and homeotropic anchoring on wavy boundaries. The equilibrium position is above the troughs (x = λ/4) for all anchoring conditions an non-zero amplitudes B₀. Center of the trough shown as the black line and errorbars correspond to the standard deviation of the position. Colloid radius R = 5a; wall amplitude B₀ = 2.5a; wall wavelength λ = 20a; and average channel height 〈h〉_x = 20a.

While the colloid's equilibrium position is centered on the troughs in all cases, the locations of the colloid-associated defects differ (Fig. 4). When the colloid is confined between two wavy walls with strong homeotropic anchoring, the two accompanying −1/2 defects are found immediately to the left and right of the colloid along the centerline (Fig. 4a). However, for the same colloid initialized in the trough of a wavy channel with planar anchoring, the defects take a diagonal configuration (Fig. 4b), rotating off the centerline to the points where the walls are closest to the colloid. By residing in the smallest space between the planar walls and the homeotropic-anchored colloid, the defect pair can reduce the nematoelastic free energy. Since there are four equal points for which the distance between surfaces are minimal, the orientation of the diagonal configuration is spontaneous. When the colloid is initialized off center from at the crest between planar-anchored walls (Fig. 4c), the diagonal pairing of defects is no longer the preferred configuration. Instead the system prefers to have the two defects shifted towards the closest boundary. This configuration lowers the instantaneous deformation free energy by placing the pair of defects at the two minimal separation points between surfaces but the state is not stable. In this state, the colloid moves away from the crests towards the center of the trough. In this non-diagonal configuration, the defect separation is not as large as when they are found at polar opposite locations. Once the colloid approaches the center of the troughs one of the defects is able to escape the deformation barrier cost of crossing to the opposite side. Unlike the case of one wavy and one plane wall (Fig. 2), the crest is not a stable point (Fig. 4d).

Fig. 4 Colloid with homeotropic anchoring between two wavy walls with various anchoring conditions. (a)–(c) The defect positions (green dots) for homeotropic and planar anchoring on the wavy walls. The green dots show the defect positions. (a) Homeotropic anchoring. (b) Planar anchoring with the colloid is centered within the trough. (c) Planar anchoring with the colloid offset to the right. (d) Trajectory of a colloid initiated above a crest. The magenta and cyan dots show the start and end positions respectively. Colloid radius R = 5a; wall amplitude B₀ = 2.5a; wall wavelength λ = 20a; and average channel height 〈h〉_x = 20a.

3.2 Flowing nematodynamics

The previous section investigated the diffusive dynamics of colloids through the distortion free energy landscape of the nematic in the absence of externally driven flows (Section 3.1). We now consider pressure-driven flow of a nematic liquid crystal through a wavy microfluidic channel without colloids.

Wavy channel walls complicate the hydrodynamics compared to plane walls (Appendix D), since the flow field v̲_f has both x̂- and ŷ-components due to the changing channel height h(x) (Fig. 5a). Increasing the amplitude increases the surface area of the no-slip walls, which in turn increases hydrodynamic resistance. For this reason, the elution of fluid is expected to decrease compared to plane channels for constant pressure gradient. To quantify this, we compared the average, maximum and minimum velocities along the centerline of wavy channels for different amplitudes B₀. Consistent with Appendix D, there is nearly no difference between the flow velocities of isotropic or nematic fluids through the channels (Fig. 5b). The average flow through a channel decreases monotonically with increasing amplitude for both nematic and isotropic fluids (Fig. 5b), suggesting that the hydrodynamic resistance increases proportionally in a manner that is not strongly dependent on nematicity. This is analogous to rough walls in plane channels, for which the flow through channels has been found to decrease non-linearly.⁹⁰

Fig. 5 Fluid velocity v̲_f(r̲) field for a channel with wavy walls. (a) Time-averaged velocity field. Color maps the x̂-component of the velocity in units a/τ, and the arrows represent the velocity vectors. Wall amplitude B₀ = 2.5a; wall wavelength λ = 20a; average channel height 〈h〉_x = 20a; and pressure gradient ; U = 10k_BT; strong planar anchoring. (b) The average, maximum and minimum of the velocity for flow through wavy channels of different amplitudes. Solid lines for nematic phase with U = 10k_BT and dashed lines for isotropic phase flow.

The difference between the maximum and minimum flow rates generally increases with the amplitude (Fig. 5b). Because N-MPCD is a coarse-grained algorithm for fluctuating hydrodynamics, a difference exists even for B₀ = 0. The maximum velocity occurs between the crests where the channel is narrowest and the minimum between troughs (Fig. 5a). This is consistent with previous theoretical and numerical work on flow of isotropic fluids through wavy channels.⁹¹ Since the maximum velocity occurs between the crests and the minimum between troughs (Fig. 5a), an advected colloid will speed up and slow down with the fluid at these points, in the absence of nematoelastic forces, Since the nematic forces also tend to drive the colloid towards the center of the troughs, these dynamics can be expected to be further emphasized.

3.3 Lock-key microfluidics

Having explored diffusive lock-key dynamics in wavy channels in the absence of advection (Section 3.1) and flows through wavy channels without colloids (Section 3.2), we now turn our attention to lock-key microfluidics by studying the advective dynamics of nematic colloids with homeotropic anchoring confined within wavy channels with planar anchoring. To investigate the colloidal transport rate through the channel, we track the position of a colloid through the channel (Fig. 6).

Fig. 6 Example trajectories for small (B₀ = 1.5; left), intermediate (B₀ = 2.5; center) and large amplitudes (B₀ = 4; right). Trajectories are colored by simulation time (short times in green and long in blue).

We begin by considering individual trajectories of colloids. Three example trajectories are shown in Fig. 6 for small (B₀ = 1.5; left), intermediate (B₀ = 2.5; center) and large (B₀ = 4; right) amplitudes. In the small amplitude case (Fig. 6; left), the stick slip dynamics are apparent in the x̂-component of the trajectories (Fig. 7a). Individual trajectories are seen to have periods within a trough (sticking events) punctuated by abrupt hops over the crest (slipping events). As the amplitude is increased (Fig. 6; center), the sticking duration increases substantially and slipping events become less frequent (Fig. 7a). Indeed by large amplitudes (Fig. 6; right), the colloids lock into their initial trough and are never observed to slip forward over the crest into the next trough (Fig. 7a).

Fig. 7 Colloid trajectories for different wall amplitudes B₀ with homeotropic anchoring on the colloid and planar anchoring on the walls. Colors specify the amplitude B₀. Black sinusoids on the right and left present undulations of the wavy walls. (a) Ten (n = 10) individual trajectories for each B₀ (thin lines), along with cumulative probability density function t(x) (eqn (1); Fig. 7c). (b) Ensemble averaged x-positions of colloids as a function of time. (c) Ensemble averaged cumulative probability density function of position scaled by time (eqn (1)).

To better quantify the elution rate, we next consider the ensemble average x-position of the colloid over time changes with the amplitude B₀ (Fig. 7b). For B₀ = 0, the transport occurs at a constant rate of advection. As B₀ increases, the average rate decreases. By , the colloids are no longer advected through the channel. By averaging over many simulation runs, the elution appears to be rather smooth for all B₀. However, we recognize that the average does not fully capture the dynamics seen in the individual runs, which exhibit clear a stick–slip dynamics.

To better measure the stick–slip dynamics of the colloids, we consider the probability density function (PDF) p(x) of finding a colloid at a given position x. The PDF can be integrated to give a cumulative probability, which gives the average time taken to arrive at a point

where α = δt/n is a scaling parameter converting the cumulative probability to a time for n numerical realizations and time step δt. The time t(x) better represents the stick–slip dynamics of each colloid (Fig. 7a). For an amplitude of B₀ = 0, the motion is completely linear and the colloids simply advect with resistance independent of x-position, as expected. For amplitudes B₀ smaller than the colloid moves along the channel in a stick–slip manner (Fig. 7c). Stick–slip dynamics are only observed in nematic fluids and not in flowing isotropic fluids (see Appendix E). We can see that the time the colloid spends “sticking” to the docking site of the trough center increases with amplitude, until the repulsive forces from the crests become large enough to fully stop the colloid from advancing. In the case of B₀ = 2.8a, the colloid often locks in place for the entire simulation and other times hops one or two troughs over the course of the entire simulation (Fig. 7a). The colloid permanently sticks at an amplitude of and the colloids are never observed to advance (Fig. 7c).

We have discussed colloids moving through wavy channels with planar anchoring, but colloids eluting through channels with homeotropic walls exhibit comparable stick–slip dynamics (data not shown). The homeotropic cases do not show qualitative differences from Fig. 7. Indeed, permanent locking occurs at = 3.0a—only slightly above = 2.9a, the point where it occurred for planar anchoring. Later, we will explicitly compare measures of colloid velocity for planar and homeotropic wavy walls and see that they are qualitatively similar.

The probability density function of finding a colloid at a given position x along the channel explains the trajectories observed for planar anchored wavy walls. The PDFs for three different amplitudes in planar anchored channels demonstrates the three different behaviours (Fig. 8). Fig. 8a shows the typical shape for the low amplitudes, where the peak is slightly offset from the center of the trough at x = λ/4. The PDF is highly biased and skewed downstream of the trough center. In this case, the colloid continuously advects through the channel. This changes as the amplitude increases. At intermediate amplitudes, the probability densities look close to the point where colloids mostly stop moving (stick) but still can hop over crests from time to time (slip) (Fig. 8b). Both the height and location of the peak changes with the amplitude. The distribution becomes bi-modal, with both PDF peaks offset from the trough's center. At the largest amplitudes (Fig. 8c), the probability density is a single peak, demonstrating that the colloid has completely “locked” into the deep free energy minimum at the trough location of x = λ/4. In this case, the probability densities are single peaks, with narrow variance but with a mean position that is slightly offset downstream of x = λ/4, due to the pressure gradient acting on the locked colloids. The biggest difference from smaller amplitude systems is that the probability density is zero everywhere away from the peak.

Fig. 8 The probability density function (PDF) of finding the homeotropic colloid at a given position x for three different amplitudes B₀ = {1.5, 2.5, 4} with planar anchoring on the wavy walls. Magenta line shows a fit of three superimposed Gaussians. The black curve illustrates undulations of the wavy walls. (a) B₀ = 1.5. (b) B₀ = 2.5. (c) B₀ = 4.

The cumulative probabilities can be used to measure the average axial velocities of colloids v_c(x) at each position x. The cumulative probability gives the time as a function of position, and so the velocity along x is found by

where t(x) is given by the cumulative distribution function (Fig. 7b). Since the cumulative probability density can be calculated from the probability densities p(x), the velocities are found to be

v_c(x) = [αp(x)]⁻¹, (3)

where α accounts for the normalization of p(x). The colloidal velocity data is stochastic, due to thermal diffusion.

The stick–slip dynamics exhibit three different scenarios (Fig. 9). Either the colloid advects smoothly with the flow (Fig. 9a) or the dynamics are stick–slip (Fig. 9b) or the colloids locks into a trough (Fig. 9c). For the continuous advection case (Fig. 9a), the maximum velocity occurs as the channel widens. This coincides with the region where the nematoelastic forces act in the same direction as the flow. The minimum occurs slightly offset from the bottom of the trough (Fig. 8). For the locked scenario (Fig. 9c), the colloid stays in the trough and the velocity is essentially zero. In the locked case, the measured probability p(x) is zero away from the trough, with the sticking time appearing to have diverged; however, longer simulations may eventually find rare crossing events. Any such rare crossing events would necessarily be accompanied by an unlikely thermal fluctuation in colloid speed to overcome the crest (eqn (3)).

Fig. 9 Homeotropic colloid velocities at different positions x for two different amplitudes B₀ with planar anchoring on the wavy walls. Magenta line shows the estimated speed from the probability density. (a) B₀ = 1.5. (b) B₀ = 2.5. (c) B₀ = 4.

The speed profiles lead to interesting average colloid dynamics, which demonstrate that homeotropic or planar wall anchoring only change quantitative details and not the qualitative behaviour. To quantify the stick–slip translation of the colloid, we track the average speed of the colloid 〈v_c(x)〉, the maximum max[v_c(x)] and the minimum min[v_c(x)] (Fig. 10). Colloids in channels with homeotropic anchoring elute more quickly than colloids in channels with planar walls (Fig. 10a). The average speed 〈v_c(x)〉 of both decreases as amplitude increases until hitting zero at . In general, the average speed of the colloid decreases with an increasing amplitude until it reaches zero (Fig. 10a). This occurs for two reasons. Firstly, as the amplitude increases the flow through the channel decreases even in the absence of the colloid (Fig. 5b). Secondly, the strength of the elastic forces increases with amplitude and so the time spent crossing the free energy barrier increases.

Fig. 10 Homeotropic colloid velocity measures with varying amplitude B₀ of the wavy boundaries with planar and homeotropic anchoring. (a) Average colloidal velocity 〈v_c(x)〉/V. (b) Maximum colloid speed max[v_c(x)]/V. (c) Minimum colloid speed min[v_c(x)]/V.

The average elution rate is different for planar and homeotropic anchoring on the walls. While the general shapes are the same, suggesting similar behaviours, the elastic forces from the wavy boundaries appear to be weaker for homeotropic anchoring than for planar. This is likely due to the different positioning of the defects. In the planar case, when the colloid is situated directly above the troughs, there is more space for the defects, compared to when it is situated above the crest where they are shifted closer to the colloid. However, for the homeotropic case, the position of the colloid has no effect on the defect positions—they always reside directly up- and downstream (Fig. 4a). The shape of the colloid matches much better with the walls when situated above the trough, meaning it is a more stable point. Thus, the mechanisms responsible for the motions of the colloid are different in the different anchoring cases.

The maximum speed of the colloids behaves qualitatively differently from the average. It suddenly decreases to max[v_c(x)] = 0 (Fig. 10b). This sudden drop occurs at the same as the average going to zero and signals the transition to a locked state that no longer moves along the channel. This once again illustrates the two distinct states of stick–slip and locked, with a discontinuous transition between them at . Prior to the sudden stop, max[v_c(x)] increases with the amplitude, which is a consequence of Bernoulli's principle. The situation is nearly identical for homeotropic and planar anchoring (Fig. 10b). Similarly, we can consider the minimum colloidal speed (Fig. 10c), which likewise shows similar behaviour for homeotropic and planar anchoring on the wall. These results indicate that homeotropic and planar anchoring result in very similar average elution rates and as well as similar maximum and minimum colloid velocities.

However, the minimum speed decreases with the amplitude. The minimum speed occurs at different points along the channel for different anchoring conditions. For homeotropic anchoring, the colloid has its minimum-velocity position much higher up on the side of the trough compared to the planar anchoring. Thus, it is nematic interactions with the wavy wall that set the position of the minimum-speed. Through no-slip boundary conditions, the colloid also impacts the average speed of the fluid. As the nemoelastic forces push the colloid towards a free energy minimum, the colloid applies a drag force to the fluid. When the colloid sticks or slows between the trough and the crest, the nematic interactions work against the pressure gradient and the fluid slows (Appendix F). On the other hand, the colloid speeds up once it crosses the barrier due to the nematic force pushing it towards the next free energy minimum. During this slip phase of motion, the colloid drags the fluid forward and increases the fluid flow (Appendix F).

In addition to comparing planar and homeotropic anchoring on the walls, we also compared different mean field potentials U for planar anchoring on the walls. The maximum colloid speed has the same general shape as we previously found (Fig. 11). As the flow through the channel is independent of the mean field potential, this is expected. However, stronger mean field potentials increase the nematoelastic interactions and so the velocities must be scaled by U, which does not collapse the curves for small amplitudes because the fluid flow is independent of U. Most clearly, the transition amplitude from stick–slip to locked moves to smaller values as the mean field potential U increases (Fig. 11). This suggests a mechanism by which fine control of the stick–slip transition to locking can be controlled.

Fig. 11 Maximum homeotropic colloid speed max[v_c(x)]/V in wavy channels with planar anchoring. Sudden transition from stick–slip dynamics to locking. Same as Fig. 10b but for various mean field potentials U. The velocity has been multiplied by the mean field potential, U, to collapse the curves.

4 Conclusion

In this work, we have considered a simple lock-key microfluidic system. We chose a colloid size commensurate with the walls' wavelength, creating periodically repeating barriers and docking sites, which exhibit interesting dynamics arising from the competition between nematoelastic free energy minimization and nematohydrodynamic advection. To tackle these complex dynamics, we first neglected flow (Section 3.1). In wavy channels, colloids with homeotropic anchoring migrate to the widest point in the channel (the troughs), regardless of whether the anchoring conditions on the wall are homeotropic or planar (Fig. 3). However, the pair of defects associated with the colloids take a different configurations for the different anchoring (Fig. 4). Our results are not comprehensive and indicate that future work is required to fully understand how wall structure and defect configuration lead to equilibrium colloid positions. We then considered flows of nematics within wavy channels in the absence of the colloids (Section 3.2). Since the transport coefficient is constant, the flow rate through the channel is observed to decrease as the the surface area of the boundary surface increased with increasing amplitude of the wavy channel (Fig. 5).

Both these results suggest the qualitative dynamics of colloid transport slows as the amplitude to the wavy walls increases, which is indeed demonstrated in simulations (Fig. 7). We find the elution slows because of stick–slip dynamics, with colloids “sticking” to the docking sites for longer durations as a function of boundary amplitude (Fig. 7 and 8). Despite this net decrease in the transport rate, the maximum velocity, of when the colloid hops over the crests from docking site to docking site, increases with amplitude until the colloid suddenly locks into the troughs and no longer moves (Fig. 9). Modifying the material parameters of the nematic fluid can shift this lock transition point (Fig. 11), with nematoelastic forces increasing with the mean field potential to nonlinearly shift the locking transition to smaller amplitudes. The transition between locking and slipping is essentially a type of elastic ‘snap-through’ event,^92,93 which would suggest a critical slowing and may allow simple-but-accurate models of the stick–slip dynamics. Such considerations might lead to the development of quantitative methods to classify different trajectories, as well as techniques to accurately measure the fraction of time a colloid spends sticking or slipping. We hope these results will motivate analytical work on the nematic interactions between mobile colloids and complex confining walls that go beyond interactions between colloids and plane walls.⁹⁴

This study has considered the effect of wavy wall amplitude on the stick–slip dynamics of eluting colloids in 2D. In addition to amplitude there are many more parameters that are expected to impact elution rate. The colloid size, channel height and wall wavelength are all length scales whose combinations will modify elution rate. The geometry could be further complicated by allowing the phase of the top and bottom wall to differ. Additionally, the strength of the wall/colloid interactions could also be tuned by modifying the Frank elastic constant of the nematic or the anchoring strength on the walls or colloids, just as the driving force could be changed by varying the pressure gradient. State diagrams for stick–slip or locking behaviors could be constructed for each pairing of such parameters. Future studies can exploit Appendix B to simulate three dimensional systems confined by wavy walls or study the role of colloid surface structure (Appendix C).

Our results show that wavy walls not only have potential to direct colloids to specific docking sites but also to control site-specific resting duration and intermittent elution, allowing autonomic temporal control. They suggest how experimental systems can be engineered to hold particles in place for a given time by combining surface structure and advective flow. By coupling site-specific resting duration with position-dependent pressure gradients, future studies could explore feedback loops to reinforce or diminish stick-and-slip motion. More complex surface structures could be designed such that particles stick for longer or shorter periods at different docking sites. Such site-specific resting times could be used to illuminate or irradiate particles for a given time without switching the light source. By combining these lock-key microfluidics with different microreactors at each docking site, efficient series chemical reactions could be designed to occur for autonomous-but-responsive durations.

Conflicts of interest

There are no conflicts to declare.

Appendix

A Numerical methods

This study utilizes a nematic multi-particle collision dynamics approach, which is chosen for its ability to simulate moderate Péclet numbers, mobile colloids, complex boundaries and fluctuating nematohydrodynamics.

A.1 Multi-particle collision dynamics algorithm. The MPCD algorithm essentially consists of two steps:^74,75 (i) streaming and (ii) collision.
A.1.1 Streaming step. The streaming step updates each particle position assuming ballistic motion

r̲_i(t + δt) = r̲_i(t) + v̲_i(t)δt, (4)

where r̲_i is the particle position, t is the current time and δt is the time step.
A.1.2 Collision step.
A.1.2.1 Momentum exchange. The particle velocities are updated by the collision operator Ξ_i,c(t) for particle i in cell c at time t, which is used to update the velocitieswhere v̲_c = 〈v̲_i〉_c is the center of mass velocity of the cell and 〈·〉_c is the average within cell c. A modified Andersen-thermostatted collision operator^77,78 is employedwhere each is a randomly generated velocity with magnitude drawn from a Boltzmann distribution with energy k_BT. To conserve momentum and keep the average velocity of the cell constant, the average of the random velocities is subtracted. The last term ensures that angular momentum is conserved by removing any spurious angular velocity generated by the collision operation. The moment of inertia of the cell is , where is the identity matrix, r̲_c = 〈r̲_j〉_c is the center of mass position and . The spurious angular momentum can arise from the randomly generated velocities or changes in orientation which will be discussed below.
A.1.2.2 Acceleration. Pressure gradients can be written as an effective external acceleration and the external acceleration g̲ can be included in the collision operator⁹⁵ as an impulse . This approach has been employed to simulate flows around obstacles,^80,96 colloid/polymer transport,^4,97,98 and rheotaxis of microswimmers.⁹⁹
A.1.2.3 Orientation exchange. Particle orientations are updated according to a nematic collision operator Ψ_i,c(t) for particle i in cell c.⁶⁴ After the collision, the updated orientation is

u̲_i(t + δt) = Ψ_i,cu̲_i(t). (7)

The collision operator is a rotation operator causes the orientation of each N-MPCD particle to change as u̲_i over the time step δt. The collision operator and resulting re-rotation can be decomposed into two contributions from the thermostat (δu̲^MS_i/δt) and from the flows (δu̲^J_i/δt).

1. A local thermostat: the local thermostat draws random orientations from the Maier–Saupe distribution centered about n̲_c

p_c(u̲_MS) ∝ exp(US_c(n̲_c·u̲_MS)²/k_BT), (8)

where U is the mean field potential, which determines how strongly the orientations align. This operation updates all particles orientations without changing the overall cell director. To determine each cell's director n̲_c, the tensoris calculated. The largest eigenvalue of is the scalar order parameter S_c and the corresponding eigenvector is parallel to the local director n̲_c for the cell. The change in orientation is δu̲^MS_i = u̲_i(t) − u̲^MS_i.

2. Rotations due to flows: the nematodynamic response models the nematogens' response to velocity gradients through Jeffery's equation

for a bare tumbling parameter ξ and heuristic shear coupling coefficient χ in a flow with shear rate tensor and vorticity tensor . The director dynamics are coupled back to the velocity field through the change in angular momentum , where γ_R is a rotational friction coefficient.^64,67,69 By choosing a sufficiently small value of γ_R backflow effects can be kept small, even in the limit of small Ericksen numbers.

B Boundary conditions

The main text considers simulations performed in 2D. In this section, we present the boundary conditions generalized to 3D and state the simplifying condition for the 2D surfaces used in the simulation results. To create the wavy channel, particles are subjected to boundary conditions, each of which is composed of two parts: the boundary surface and boundary rules.

B.1 Boundary surfaces. Each boundary surface labelled by index b is represented as a mathematical surface where . The baseline form employed for can represent planes, ellipsoids and squircleswhere r̲ = (x, y, z) is the position on the surface, q̲_b = (q_b,x, q_b,y, q_b,z) defines the position of the bth boundary, and A̲_b = (A_b,x, A_b,y, A_b,z), p̲_b = (p_b,x, p_b,y, p_b,z, p_b,R) and R_b define the shape.

In this study, the surface has been extended to allow for wavy boundaries by adding a second term, which allows for waves along two perpendicular directions with amplitude B_b,0 and frequencies B_b,1 and B_b,2. We represent these surfaces as

where defines the basic surface given by eqn (11), while and define the orientations of the waves, and thus depend on . For 2D simulations, B_b,2 = 0.

B.1.1 Wavy planes. For planes, p̲_b = 1̲, A̲_b is the normal vector and R_b = 0. For wavy planes, the functions and can define (i) an egg carton surface (Fig. 12a) or (ii) a corrugation (Fig. 12b). To get the egg carton-like surface, must be mutually perpendicular. A decision has to be made to make independent of z and take the formThe plane is normalized to standardize the value of the frequency B_b,1. This choice then fixes to beAs previously stated, B_b,1 and B_b,2 specify the frequency of the waves. If they are equal, the pattern in Fig. 12a emerges with even waves in both directions. However, if they are not equal, oblong ‘bumps’ appear. The corrugated plane emerges if B_b,2 = 0 for waves in the other direction. In those cases, the other B_b,1 then defines the frequency of the waves. Although this study will focus on wavy planar walls, wavy colloids and cylinders are also possible, as detailed in Appendix C.

Fig. 12 The types of wavy plane walls that can be implemented. (a) Egg carton plane wall. (b) Corrugated plane wall.

While the wavy surfaces are generalized to work in 3D, they work without modification in 2D. For simplicity, it is assumed that 2D simulations occur in the xy-plane, such that A_b,z = 0 for all boundaries b. Additionally, only one cosine wave is needed and so we set B_b,2 = 0 for all cases. Thus, B_b,2 = 0 is not discussed in the main text.

B.2 Boundary rules. The bth surface boundary exists at , thus a particle is defined as being within the simulation domain if for particle position r̲_i. If a particle is found to be outside the boundary, a set of boundary rules define how the generalized coordinates (position r̲_i(t), orientation u̲_i(t) and velocity v̲_i(t)) are transformed. The boundary rules are applied to the components normal v̲_b and perpendicular t̲_b to the boundary (Appendix C.3).

Periodic boundary conditions can be achieved by shifting the position r̲_i(t) in the normal direction by and/or in the tangential direction by . After crossing the surface b, particle i's position is transformed as . Likewise, the particle's velocity may be transformed. The velocity is altered by multiplicative parameters and , such that upon colliding with the surface the velocity becomes where and are the normal and tangent projection operators. Lastly, the orientation is also multiplicatively transformed which is then rescaled back to a unit vector.

The simplest boundaries employed in this study are periodic boundary conditions. Particles that cross a periodic boundary are simply shifted by the length of the simulation domain in the direction normal to the boundary, without change to their velocity or orientation. So while and where L_j is the system size in the direction j.

Solid walls are boundaries that do not apply a shift but do transform the direction of the velocity. For perfect-slip boundaries, the particles are reflected and and ; whereas, for no-slip conditions bounce-back rules are applied with . In order for the particle to not violate the boundary, the particle is rewound to the location where it crossed the boundary, and from there it streams for the rest of the time step using its new velocity.

Controlling the orientation when a particle crosses a boundary allows for different nematic anchoring conditions for solid boundaries. If both μ_b,ν and μ_b,t are set to 1 the particle's orientation does not change and no anchoring conditions are enforced. However, by setting the tangential component to μ_b,t = 0 homeotropic anchoring is achieved. Similarly, if μ_b,ν = 0 planar anchoring is enforced on the boundary.

B.3 Strong anchoring. For solid boundaries, the above boundary rules do not fully enforce perfect no-slip.⁸² This is because when the boundaries intersect cells, they have a lower than average density since part of the cell is excluded, which lowers the viscosity of the cell. To combat this, ‘phantom’ particles are added to make up the correct density. These particles are only used when generating the new velocities during the collision stage.

In a similar fashion, the boundary rules for anchoring generate weak anchoring because only particles that have crossed the boundary have the anchoring condition applied, and are then mixed with particles that have not.⁸³ In order to generate strong anchoring conditions, all the particles in cells that intersect are subjected to the boundary rules described in Section 6.2. Because the entire cell now essentially has the anchoring condition applied, this generates stronger anchoring.

B.4 Colloids as moving boundaries. Hard colloids can be simulated as spherical (or discoidal in 2D) boundaries. However, the boundaries must be able to move as a response to the surrounding fluid, . Firstly, when a particle collides with the boundary, the impulse of the collision must be applied to the moving boundary to conserve both linear and angular momentum.⁴ Additionally, when anchoring is used, the change in particle orientation applies a torque on the colloid. In this study, the anchoring condition on the colloids are always strong homeotropic anchoring. Moving boundaries also needs to be able to collide with other boundaries, both moving and fixed. For spherical or discoidal colloids of radius R_b, the algorithm checks if it has collided with another solid boundary and applies the same rules as when colliding with an MPCD particle.

C Wavy spheres, ellipsoids and cylinders

In addition to corrugated and egg-shell structured planar walls, wavy spheres, ellipsoids and cylinders are possible viaeqn (11)–(14).

C.1 Spheres and ellipsoids. For spheres, p̲_b = 2̲, A̲_b = 1̲ and R_b is the radius of the bth sphere. For ellipsoids, A̲_b is the inverse of the bth ellipsoid's focii. There are two different directions in which a sphere can be wavy: longitudinal and latitudinal (Fig. 13a and b). These can also be combined to create a sphere with egg carton waves on it (Fig. 13c). The longitudinal and the latitudinal waves correspond to and respectively, which take the formsTo create longitudinal waves on a sphere B_b,2 = 0, and conversely B_b,1 = 0 for latitudinal waves. Neither is 0 for the egg carton shapes on the sphere. The B_b,i correspond to the number of waves around the sphere. Thus, due to the closed surface nature of spheres, all B_b,i must be even numbers in order to have complete waves. The surface will still be closed for non-integer and odd B_b,i but they will result in sharp points. These equations work for ellipsoids as well, however the waves are not perfectly even.

Fig. 13 Three types of wavy spheres. (a) Longitudinal. (b) Latitudinal. (c) Egg Carton.

C.2 Cylinders. Cylinders are a straightforward extension to spheres, where the surface is independent of one of the directions. Thus, either A_b,x, A_b,y or A_b,x is 0 but otherwise follows the same form as spheres and ellipsoids. There are three types of wavy cylinders: waves around the cylinder, waves along the cylinder, and an egg carton pattern on the cylinder which combines both (Fig. 14). The two functions and for cylinders have a cyclic property so if A_b,x = 0, then (i, j, k) = (x, y, z), if A_b,y = 0: (i, j, k) = (y, z, x); and if A_b,z = 0 then (i, j, k) = (z, x, y). So then for any choice A_b,i = 0 the wavy surfaces are given byHere, defines the waves around the cylinder (Fig. 14b) and defines waves along the direction of the cylinder (Fig. 14c). By setting the opposite B_b,i to 0, these can be achieved. If both B_b,i ≠ 0, the egg carton pattern emerges on the cylinder (Fig. 14a).

Fig. 14 The three types of wavy cylinders. (a) Egg carton. (b) Polar undulations. (c) Axial undulations.

C.3 Calculating the normal of the boundary. For some of the boundary rules, the normal and tangential directions of boundary b are needed. As the component of a vector tangential to the surface is just the vector minus the normal component, only the normal direction v̲_b(r̲) needs to be calculated asTaking the derivative of eqn (12) givesFor all surfaces, the derivative of is given byHowever, the derivatives of and vary depending on the surface type.

C.4 Planes. The normals for wavy planes can be found with the derivatives of and the derivatives of where we conveniently define

C.5 Spheres and ellipsoids. The surface functions for spheres and ellipsoids include arctan functions, which make the derivatives of and the derivatives of where we define

C.6 Cylinders. Cylinders are the most complicated boundary considered here due to the possibility of several orientations. The derivatives have a cyclic property so if A_b,x = 0, then (j, k, l) = (x, y, z), if A_b,y = 0: (j, k, l) = (y, z, x); and if A_b,z = 0 then (j, k, l) = (z, x, y). Letting the derivatives are

D Flow through plane channels

Flowing isotropic fluids at low Reynolds number obey plane Poiseuille flow with parabolic flow profiles,¹⁰⁰ which is a direct consequence of the the no-slip condition at the walls of the channel. Though the coupling between flow profile, local director and anisotropic dissipation can induce complications,¹⁰¹ N-MPCD simulates linearized fluctuating nematodynamics with isotropic viscous dissipation.⁷⁹

The director can discontinuously transition from a homogenous alignment to deformed states.¹⁰² For strong homeotropic anchoring and sufficiently small pressure gradients, the director remains uniformly homeotropic across the channel; however, for sufficient flow the director field bows into the Bowser state,¹⁰² in which the field curves slightly in the direction of the flow (Fig. 15a). The Bowser state transitions (Fig. 15b) to the the Dowser state (Fig. 15c) at sufficiently strong flows, for which the director fully aligns with the flow at the center (Fig. 15d). The Bowser and Dowser states arise from the competition between nematoelastic free energy and flow strength.¹⁰³ For the pressure gradient and channel height used in this study, the director is found to be in the bowed Bowser state. The planar anchoring case is not as interesting since the boundary conditions and strain rate coalign, though at sufficiently stronger flow the near-wall director may tilt inward as expected.¹⁰⁴

Fig. 15 Nematic director profiles for different flows in a plane channel of height h = 30a (characteristic length L = h/2). (a) Bowser state for weak flow (g_x = 0.0002a/τ²). (b) Transition state (g_x = 0.002a/τ²). (c) Dowser state for strong flows (g_x = 0.02a/τ²). (d) Director profile n_x = x̂·n̲ across the channel for the Bowser, transition and Dowser states.

Since the flow profiles remain parabolic (Fig. 16a), the viscosity can be calculated directly from the average velocity of the nematic fluid 〈v_f〉 and the plane Hagen–Poiseuille equation to be

To measure the viscosity, the pressure gradient is kept sufficiently low to be in the Bowser state and viscosity is found to be independent of the values of mean field potential (for U = {5, 10, 20}k_BT) and anchoring conditions (Fig. 16b), as expected since the nematogens are modeled as point-particles with a direction, rather than actual rod-like particles. An isotropic viscosity of μ = [9.6 ± 0.4]μ₀ is found. This measurement gives the characteristic velocities due to nematicity and pressure to be Ṽ = K/μL = [11.8 ± 0.8]a/τ and V = L²ρg_x/μ = [0.208 ± 0.009]a/τ, respectively. This establishes the dimensionless Ericksen number for this study to bewhich indicates elastic forces dominate viscous forces. Likewise, the dimensionless Péclet number is approximately Pe = 〈v_f〉L/D ∼ 10², which reflects that advection is more significant than diffusion in plane channels. The small Er/large Pe regime is where we expect topological microfluidics to be most interesting. In these simulations, elastic forces dominate over both velocity and diffusion, Er/Pe ∼ 3 × 10⁻⁵. This can be seen explicitly by comparing trajectories in a nematic solvent, which are near the middle line, to trajectories in an isotropic solvent, which can diffuse away from the center (see Appendix E).

Fig. 16 Flow profiles in planar channels. (a) The flow profiles for different pressure gradients ρg_x, where g_x is the component of the external acceleration in the x̂ direction. The flow obeys parabolic Poiseuille flow for all pressure gradients considered here in a plane channel of height h = 30a (with L = h/2). (b) The viscosity of the fluid for different mean field potentials and wall anchoring in a plane channel of height h = 20a.

E Absence of stick–slip in isotropic fluid flows

Simulations of a colloid in an isotropic fluid advecting through a wavy channel do not exhibit stick–slip dynamics (Fig. 17). While the nematic case exhibits clear periods of movement (slip) and cessation (stick), the isotropic case exhibits a steady elution rate (Fig. 17a). Furthermore, nematic repulsion between the colloids and the walls keeps the colloid close to the channel centerline (Fig. 17b). On the other hand, the colloid in an isotropic fluid is free to diffuse away from the centerline, sampling different local velocities and exhibiting greater Taylor dispersion. In this way, the lock-key dynamics investigated here differ from the Brownian dynamics of a driven diffusive colloid moving through a periodic landscape.

Fig. 17 Example trajectories for a colloid suspended in a nematic compared to in an isotropic fluid. The wall amplitude in both cases is B₀/L = 1. The nematic case has homeotropic anchoring on the colloid and planar anchoring on the walls. (a) x̂-component of the colloid trajectories as a function of time. (b) Colloid trajectories.

F Impact of colloid dynamics on fluid flow

The colloid elutes due to the drag force that the fluid exerts on it, but it also moves in response to the nematoelastic forces. Therefore, at any instance the colloid may resist or assist the fluid flow. The spatially averaged flow of the fluid is not constant with time during the stick–slip motion. When the colloid sticks, the combination of colloid drag and nematoelastic forces slow the fluid. When the colloid slips, the nematoelastic forces drive fluid flow. We measure the total average fluid velocity 〈v_f〉 as a function of the colloid position (Fig. 18). A minimum average fluid velocity occurs just past the trough, where defects oppose elution. Here, nematic interactions work against the pressure gradient and the fluid slows. During the slip phase, the colloid speeds up due to nematic forcing and the fluid flow reaches a maximum as the colloid is forced over the crest by elastic interactions with the walls. As the wall amplitude is increased, crossing events become rarer and faster, causing the statistics to worsen for x/L > 1.

Fig. 18 Average fluid velocity 〈v_f〉 as a function of instantaneous colloid position x for wall amplitudes B₀ = {1.5, 2.5}a.

Acknowledgements

This work was supported by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement no. 851196). For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising from this submission.

Notes and references

1. Y. Liu, G. Yang, Y. Hui, S. Ranaweera and C.-X. Zhao, Small, 2022, 18, 2106580.
2. M. Sivaramakrishnan, R. Kothandan, D. K. Govindarajan, Y. Meganathan and K. Kandaswamy, Curr. Opin. Biomed. Eng., 2020, 13, 60–68.
3. A. Burklund, A. Tadimety, Y. Nie, N. Hao and J. X. Zhang, Advances in Clinical Chemistry, Elsevier, 2020, vol. 95, pp. 1–72.
4. T. N. Shendruk, R. Tahvildari, N. M. Catafard, L. Andrzejewski, C. Gigault, A. Todd, L. Gagne-Dumais, G. W. Slater and M. Godin, Anal. Chem., 2013, 85, 5981–5988.
5. J. Zhang, S. Yan, D. Yuan, G. Alici, N.-T. Nguyen, M. Ebrahimi Warkiani and W. Li, Lab Chip, 2016, 16, 10–34.
6. H. Zhang, D. Freitas, H. S. Kim, K. Fabijanic, Z. Li, H. Chen, M. T. Mark, H. Molina, A. B. Martin and L. Bojmar, et al. , Nat. Cell Biol., 2018, 20, 332–343.
7. P. Zhang, H. Bachman, A. Ozcelik and T. J. Huang, Annu. Rev. Anal. Chem., 2020, 13, 17–43.
8. A. Sengupta, C. Bahr and S. Herminghaus, Soft Matter, 2013, 9, 7251–7260.
9. O. D. Lavrentovich, Liquid Crystals Rev., 2020, 8, 59–129.
10. K. Fedorowicz, R. Prosser and A. Sengupta, Soft Matter, 2023, 19(37), 7084–7092.
11. J. Dalby, Y. Han, A. Majumdar and L. Mrad, SIAM J. Appl. Math., 2023, 83(6), 2284–2309.
12. A. Sengupta and M. G. Mazza, Liquid Crystals at Interfaces and Under Flow: Recent Advances and Trends, 2020, ch. 6, pp. 183–226.
13. I. I. Smalyukh, Ann. Rev. Condens. Matter Phys., 2018, 9, 207–226.
14. H. Stark, Phys. Rep., 2001, 351, 387–474.
15. Y. Han and A. Majumdar, J. Nonlinear Sci., 2023, 33, 24.
16. Y. Han, B. Shi, L. Zhang and A. Majumdar, IMA J. Appl. Math., 2023, 88(5), 645–676.
17. H.-Q. Chen, X.-Y. Wang, H. K. Bisoyi, L.-J. Chen and Q. Li, Langmuir, 2021, 37, 3789–3807.
18. Z. Dogic, P. Sharma and M. J. Zakhary, Ann. Rev. Condens. Matter Phys., 2014, 5, 137–157.
19. I. I. Smalyukh, Rep. Prog. Phys., 2020, 83, 106601.
20. H. Mundoor, S. Park, B. Senyuk, H. H. Wensink and I. I. Smalyukh, Science, 2018, 360, 768–771.
21. J. Paget, M. G. Mazza, A. J. Archer and T. N. Shendruk, Nat. Commun., 2023, 14, 1048.
22. Y. Yuan, Q. Liu, B. Senyuk and I. I. Smalyukh, Nature, 2019, 570, 214–218.
23. B. Senyuk, Q. Liu, P. D. Nystrom and I. I. Smalyukh, Soft Matter, 2017, 13, 7398–7405.
24. I. Muševič, Eur. Phys. J. Special Topics, 2019, 227, 2455–2485.
25. H. Yoshida, K. Asakura, J.-I. Fukuda and M. Ozaki, Nat. Commun., 2015, 6, 7180.
26. M. Nikkhou, M. Škarabot, S. Čopar, M. Ravnik, S. Žumer and I. Muševič, Nat. Phys., 2015, 11, 183–187.
27. F. R. Hung, B. T. Gettelfinger, G. M. Koenig, N. L. Abbott and J. J. de Pablo, J. Chem. Phys., 2007, 127(12), 124702.
28. C. P. Lapointe, D. H. Reich and R. L. Leheny, Langmuir, 2008, 24, 11175–11181.
29. N. M. Silvestre, P. Patrício and M. M. Telo da Gama, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 6.
30. D. L. Cheung and M. P. Allen, J. Chem. Phys., 2008, 129(11), 114706.
31. Z. Eskandari, N. M. Silvestre, M. M. Telo da Gama and M. R. Ejtehadi, Soft Matter, 2014, 10, 9681–9687.
32. Y. Luo, F. Serra and K. J. Stebe, Soft Matter, 2016, 12, 6027–6032.
33. G. Boniello, Y. Luo, D. A. Beller, F. Serra and K. J. Stebe, Soft Matter, 2019, 15, 5220–5226.
34. Óscar A. Rojas-Gómez, J. M. Romero-Enrique, N. M. Silvestre and M. M. T. da Gama, J. Phys.: Condens. Matter, 2016, 29, 064002.
35. Y. Luo, T. Yao, D. Beller, F. Serra and K. Stebe, Langmuir, 2019, 35, 9274–9285.
36. T. Yao, v Kos, Q. X. Zhang, Y. Luo, F. Serra, E. B. Steager, M. Ravnik and K. J. Stebe, Adv. Funct. Mater., 2022, 32, 2205546.
37. Y. Luo, D. A. Beller, G. Boniello, F. Serra and K. J. Stebe, Nat. Commun., 2018, 9, 1–11.
38. K. Chen, O. J. Gebhardt, R. Devendra, G. Drazer, R. D. Kamien, D. H. Reich and R. L. Leheny, Soft Matter, 2018, 14, 83–91.
39. R. Kapral, Adv. Chem. Phys., 2008, 140, 89.
40. G. Gompper, T. Ihle, D. M. Kroll and R. G. Winkler, in Multi-Particle Collision Dynamics: A Particle-Based Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids, ed. C. Holm and K. Kremer, Springer Berlin Heidelberg, Berlin, Heidelberg, 2009, pp. 1–87.
41. G. W. Slater, C. Holm, M. V. Chubynsky, H. W. de Haan, A. Dubé, K. Grass, O. A. Hickey, C. Kingsburry, D. Sean, T. N. Shendruk and L. Zhan, Electrophoresis, 2009, 30, 792–818.
42. M. P. Howard, A. Nikoubashman and J. C. Palmer, Curr. Opin. Chem. Eng., 2019, 23, 34–43.
43. A. Sayyidmousavi and K. Rohlf, Phys. Biol., 2018, 15, 046007.
44. S. Y. Reigh, M.-J. Huang, H. Löwen, E. Lauga and R. Kapral, Soft Matter, 2020, 16, 1236–1245.
45. O. A. Hickey, T. N. Shendruk, J. L. Harden and G. W. Slater, Phys. Rev. Lett., 2012, 109, 098302.
46. T. N. Shendruk, M. Bertrand and G. W. Slater, ACS Macro Lett., 2015, 4, 472–476.
47. J. Burelbach, D. B. Brückner, D. Frenkel and E. Eiser, Soft Matter, 2018, 14, 7446–7454.
48. A. Zöttl and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2018, 41, 1–11.
49. J.-T. Kuhr, F. Rühle and H. Stark, Soft Matter, 2019, 15, 5685–5694.
50. J. Clopés, G. Gompper and R. G. Winkler, Soft Matter, 2020, 16, 10676–10687.
51. A. Lamura, R. G. Winkler and G. Gompper, J. Chem. Phys., 2021, 154(22), 224901.
52. J. Clopés Llahí, A. Martín-Gómez, G. Gompper and R. G. Winkler, Phys. Rev. E, 2022, 105, 015310.
53. K. Choi, J. Lee, W. J. Choi and J. S. Myung, Polym. Eng. Sci., 2023, 63, 1032–1040.
54. Z. Wang, Z.-G. Wang, A.-C. Shi, Y. Lu and L. An, Macromolecules, 2023, 56, 2447–2453.
55. Y. M. Wani, P. G. Kovakas, A. Nikoubashman and M. P. Howard, J. Chem. Phys., 2022, 156(2), 024901.
56. R. Chen, R. Poling-Skutvik, M. P. Howard, A. Nikoubashman, S. A. Egorov, J. C. Conrad and J. C. Palmer, Soft Matter, 2019, 15, 1260–1268.
57. T. Eisenstecken, R. Hornung, R. G. Winkler and G. Gompper, Europhys. Lett., 2018, 121, 24003.
58. D. Toneian, G. Kahl, G. Gompper and R. G. Winkler, J. Chem. Phys., 2019, 151(19), 194110.
59. P. Ilg, J. Chem. Phys., 2022, 156(14), 144905.
60. P. Ilg, Phys. Rev. E, 2022, 106, 064605.
61. P. Di Cintio, M. Pasquato, L. Barbieri, H. Bufferand, L. Casetti, G. Ciraolo, U. N. di Carlo, P. Ghendrih, J. P. Gunn and S. Gupta, et al. , Proc. Int. Astron. Union, 2020, 16, 134–140.
62. P. Di Cintio, M. Pasquato, H. Kim and S.-J. Yoon, Astron. Astrophys., 2021, 649, A24.
63. P. Di Cintio, M. Pasquato, A. Simon-Petit and S.-J. Yoon, Astron. Astrophys., 2022, 659, A19.
64. T. N. Shendruk and J. M. Yeomans, Soft Matter, 2015, 11(25), 5101–5110.
65. S. Mandal and M. G. Mazza, Phys. Rev. E, 2019, 99, 063319.
66. D. Reyes-Arango, J. Quintana-H., J. C. Armas-Pérez and H. Híjar, Phys. A, 2020, 547, 123862.
67. H. Híjar, Phys. Rev. E, 2020, 102, 062705.
68. J. Armendáriz and H. Híjar, Int. J. Modern Phys. B, 2021, 35, 2150124.
69. J. Armendáriz and H. Híjar, Materials, 2021, 14(11), 2886.
70. S. Mandal and M. G. Mazza, Eur. Phys. J. E: Soft Matter Biol. Phys., 2021, 44, 64.
71. T. Kozhukhov and T. N. Shendruk, Sci. Adv., 2022, 8, eabo5788.
72. T. Kozhukhov, B. Loewe and T. N. Shendruk, arXiv, 2024, preprint, arXiv:2401.17777 DOI:10.48550/arXiv.2401.17777.
73. R. R. Keogh, T. Kozhukhov, K. Thijssen and T. N. Shendruk, Phys. Rev. Lett., 2024 , https://journals.aps.org/prl/accepted/35078Y2eR2c1468de9065e368120c71a35d9fdb0b, accepted.
74. A. Malevanets and R. Kapral, J. Chem. Phys., 1999, 110, 8605–8613.
75. A. Malevanets and R. Kapral, J. Chem. Phys., 2000, 112, 7260–7269.
76. T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 020201.
77. H. Noguchi, N. Kikuchi and G. Gompper, Europhys. Lett., 2007, 78, 10005.
78. H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 016706.
79. H. Híjar, Condens. Matter Phys., 2019, 22, 13601.
80. A. Lamura, G. Gompper, T. Ihle and D. M. Kroll, Europhys. Lett., 2001, 56, 319.
81. J. K. Whitmer and E. Luijten, J. Phys.: Condens. Matter, 2010, 22, 104106.
82. D. S. Bolintineanu, J. B. Lechman, S. J. Plimpton and G. S. Grest, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 066703.
83. L. C. Head, Y. A. G. Fosado, D. Marenduzzo and T. N. Shendruk, arXiv, 2024, preprint, arXiv:2404.09368 DOI:10.48550/arXiv.2404.09368.
84. S. Grollau, N. L. Abbott and J. J. de Pablo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 011702.
85. D. Andrienko, G. Germano and M. P. Allen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 041701.
86. D. M. Sussman and D. A. Beller, Front. Phys., 2019, 7 DOI:10.3389/fphy.2019.00204.
87. J. Fukuda and H. Yokoyama, Eur. Phys. J. E: Soft Matter Biol. Phys., 2001, 4, 389–396.
88. M. Tasinkevych, N. Silvestre, P. Patricio and M. Telo da Gama, Eur. Phys. J. E: Soft Matter Biol. Phys., 2002, 9, 341–347.
89. N. M. Silvestre, Q. Liu, B. Senyuk, I. I. Smalyukh and M. Tasinkevych, Phys. Rev. Lett., 2014, 112, 225501.
90. Y. Liu, J. Li and A. J. Smits, J. Fluid Mech., 2019, 876, 1129–1145.
91. N. F. Okechi and S. Asghar, Chin. J. Phys., 2021, 74, 144–156.
92. M. Gomez, D. E. Moulton and D. Vella, Nat. Phys., 2017, 13, 142–145.
93. B. Radisson and E. Kanso, Phys. Rev. Lett., 2023, 130, 236102.
94. S. B. Chernyshuk and B. I. Lev, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 011707.
95. E. Allahyarov and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 036702.
96. A. Lamura and G. Gompper, Eur. Phys. J. E: Soft Matter Biol. Phys., 2002, 9, 477–485.
97. C. Prohm, M. Gierlak and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 1–10.
98. A. Liu, Z. Yang, L. Liu, J. Chen and L. An, Macromolecules, 2020, 53, 9993–10004.
99. K. Qi, H. Annepu, G. Gompper and R. G. Winkler, Phys. Rev. Res., 2020, 2, 033275.
100. H. Bruus, Theoretical microfluidics, Oxford University Press, 2007, vol. 18.
101. F. Brochard, Mol. Cryst. Liq. Cryst., 1973, 23, 51–58.
102. S. Čopar, V. Kos, T. Emeršič and U. Tkalec, Nat. Commun., 2020, 11, 59.
103. T. G. Anderson, E. Mema, L. Kondic and L. J. Cummings, Int. J. Non Linear Mech., 2015, 75, 15–21.
104. V. M. O. Batista, M. L. Blow and M. M. Telo Da Gama, Soft Matter, 2015, 11, 4674–4685.

---