Lock-Key Microfluidics: Simulating Nematic Colloid Advection along Wavy-Walled Channels

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 colloids and troughs is found to be stronger for planar anchoring compared to homeotropic anchoring 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.


I. INTRODUCTION
Manipulating particle trajectories is an important feature of many microfluidic applications, including drug delivery [1], cell sorting [2] and medical diagnoses [3].Accordingly, a variety of approaches have been developed to control the trajectories of microparticles [4][5][6][7].One means of sculpting particle trajectories is by suspending them in complex carrier fluids, as in topological microfluidics [8], in which particles are guided by liquid crystalline carrier fluids [9].Liquid crystalline materials are alluring for microfluidic transport because they are highly responsive to flows [10][11][12], suspended inclusions [13,14] and confining surfaces [15][16][17].When colloidal particles are dispersed in liquid crystalline fluids, the anisotropic nature of liquid crystals give rise to emergent properties [18][19][20][21] and imposed anchoring at colloidal surfaces results in topological defects in the vicinity of the colloids, to ensure topological charge neutrality.Strong homeotropic and 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][23][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 [29], 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 * t.shendruk@ed.ac.uk 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 [32].Variations of this setup have included imposed twist [33], saw-tooth or crenellated substrates [34], and ellipsoidal [35] or fourarmed colloids [36].In addition to setting stable docking sites, structured boundaries can control dynamic trajectories.Trajectories can arise due to diffusion or external fields [37], 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 [38].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 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 colour shows local scalar order parameter field S and white dashes show the director field n.Colloid radius R = 5a; wall amplitude B 0 = 2.5a; wall wavelength λ = 20a; average channel height ⟨h⟩ x = 20a; and pressure gradient ∇P = 0.02ma 1−d τ −2 x.Snapshot from t = 20τ , early in the simulation.over elution dynamics and the ability to lock colloids in place for specific times.
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 and unit-length orientation u i to model the direction of the nematogens.The MPCD algorithm consists of two steps [74,75]: (i) streaming and (ii) collision.The streaming step (see SI A) updates each particle position r i assuming ballistic motion over a time step δt.The collision step (see SI 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 [76].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 SI 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 B T [77,78].This operator conserves momentum and thermostats the energy.
2. External forces: A pressure gradient is applied as ∇P = ρg, 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 , 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 SI A for details) [64,69].The small value chosen keeps 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 [79].
The wavy 2D channel is modelled with plane-wavy solid boundaries (see SI B).Each surface (both bounding walls and colloidal) is labelled by index b and represented as a surface on which S b (r) = 0.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-cartoon and corrugated walls (SI B).Both the top and bottom wavy walls of the channel have the same amplitude B b,0 = B 0 and wavelength λ given by B b,1 = B 1 = 2πL/λ for all wall boundaries, where L 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 (SI C).Bounceback boundary conditions with phantom particles ensure no-slip [80][81][82].Periodic boundary conditions are plane surfaces at the channel extremities.

II.1. Simulation Parameters
The simulation units are particle mass m, thermal energy k B T and cell size a.From these, the time unit is τ = a m/k B T .Additional derived units include dynamic viscosity µ 0 ≡ k B T τ /a d [78], stress k B T a −d and Frank elasticity k B T a 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 2 .The mean field potential is U = 10k B T (unless otherwise stated), placing the fluid deep in the nematic phase [64].A mean field potential of U = 5k B T 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 B T and 20k B T have scalar parameters much closer to S ≈ 1 [64].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 (SI 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 2 [64].The nematic is in the flow aligning regime with ξ = 2, shear susceptibility χ = 0.5 and rotational friction γ R = 0.01τ k B T .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/τ 2 x to generate a pressure gradient ∇P = 0.02k B T /a 1+d x.Ten repeats (n = 10) of each configuration are run for 1.5 × 10 4 τ .
The channel is modelled with plane-wavy solid boundaries along x.Both walls have the same amplitude B 0 and frequency B 1 .Planes set 60a apart define the channel walls with average normal ±ŷ.The domain size is 60a × (20a + 2B 0 ) and the frequency of the wavy walls, B 1 , is set to π, so that 3 full periods fit within the length of the channel, which gives a wavelength λ = 20a.Posi-tions of x = λ/4 and x = 3λ/4 correspond to the center of the trough and crest, respectively.The amplitude of the waves B 0 is varied between 0 and 4.5a.
The colloid is simulated as a moving circular boundary with a radius R = 5a and mass ρπr 2 .The presence of the −1/2 companion defect pair (Fig. 1) indicates that the anchoring is strong (SI B) and that the characteristic anchoring length is small compared to the colloid size [83].Colloids are initialized with their center at r c = (x, y) = (25a, 10a + B 0 ), 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 2 ρg x /µ.For the parameters used here, the nematoelastic speed Ṽ ≈ 12a/τ is much larger than the advective speed V = 0.208(9)a/τ , 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.

III. 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 ( § III.1), then proceed to consider nematohydrodynamic flows in wavy channels in the absence of the colloid ( § III.2) before putting the two aspects together to study lock-key microfluidics ( § III.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 [84] 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.

III.1. Colloid Trajectories With No Pressure Gradient
We place colloids anchoring into channels with one plane boundary wall and one wavy wall, with an ampli-

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 0 = 2.5a; wall wavelength λ = 20a; and average channel height ⟨h⟩ x = 20a.tude of B 0 = 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.These numerical trajectories are consistent with Luo et al. [37].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 little only stochastic motion in the ŷ direction.Although the colloid does not move sub- stantially 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).
Homeotropic anchoring on the colloid and planar anchoring on the walls was not investigated by Luo et al. [37], 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 [37] and with the trapping of colloidal particles at the tips of sharp protrusions whenever the colloid and surface have opposing anchoring conditions [89].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.
While the colloid's equilibrium position is centered on the troughs in all cases, the locations of the colloidassociated defects differ (Fig. 4).When the colloid with strong homeotropic anchoring, possessing two satellite −1/2 defects, is confined between two wavy walls with strong homeotropic anchoring the 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 walls, with their planar anchoring and the colloid, with its homeotropic surface, 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).

III.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 ( § III.1).We now consider pressure-driven flow of a nematic liquid crystal through a wavy microfluidic channel.
Wavy channel walls complicate the hydrodynamics compared to plane walls (SI D), since the flow field v f has both x-and ŷ-components, which depend on both y and x 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 0 .Consistent with SI 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 de- creases 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 [90].
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 = 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 [91].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 emphasised.

III.3. Lock-key microfluidics
Having explored diffusive lock-key dynamics in wavy channels in the absence of advection ( § III.1) and flows through wavy channels without colloids ( § III.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).
We begin by considering individual trajectories of colloids.Three example trajectories are shown in Fig. 6 for small (B 0 = 1.5; left), intermediate (B 0 = 2.5; center) and large (B 0 = 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).
To better quantify the elution rate, we next consider the ensemble average x-position of the colloid over time changes with the amplitude B 0 (Fig. 7b).For B 0 = 0, the transport occurs at a constant rate of advection.As B 0 increases, the average rate decreases.By B * 0 = 2.9a, 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 0 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 = 0, the motion is completely linear and the colloids simply advect with resistance independent of x-position, as expected.For amplitudes B 0 smaller than B * 0 = 2.9a, 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 SI 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 0 = 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 B * 0 = 2.9a 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 B 0 = 3.0a-only slightly above B 0 = 2.8a, 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.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.
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 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 maximum of the flow, and is in 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 (Eq.( 3)).
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 B * 0 = 2.9a.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. 5bb).Secondly, the strength of the elastic forces increases with amplitude and so the time spent crossing the free energy barrier increases.
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 B * 0 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 B * 0 .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 am-plitude.The minimum speed occurs at different points along the channel for different anchoring conditions and for the flow.For homeotropic anchoring, the colloid has its minimum-velocity position much higher up on the side of the trough compared to both the planar anchoring and the flow.Thus, it is nematic interactions with the wavy wall that set the position of the minimum-speed.Through its no-slip boundary conditions, the colloid also impacts the average speed of the fluid.As the nematic forces push the colloid towards a free energy minimum, it 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 (SI 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 (SI 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 B * 0 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.

IV. 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 ( § III.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 ( § III.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-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 flowing 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 [94].
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 SI B to simulate three dimensional systems confined by wavy walls or study the role of colloid surface structure (SI 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.

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.
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 1 is the identity matrix, r c = r j c is the center of mass position and r ′ k = r k − r c .The spurious angular momentum can arise from the randomly generated velocities δL c (t) = m r ′ j × v j − ξ j c N c or changes in orientation δL c (t), which will be discussed below.b.2.Acceleration: Pressure gradients can be written as an effective external acceleration ∇P = ρg and the external acceleration g can be included in the collision operator [95] as an impulse v i (t + δt) = v c (t) + Ξ i,c (t) + gδt.This approach has been employed to simulate flows around obstacles [80,96], colloid/polymer transport [4,97,98], and rheotaxis of microswimmers [99].
b.3.Orientation exchange: Particle orientations are updated according to a nematic collision operator Ψ i,c (t) for particle i in cell c [64].After the collision, the updated orientation is The collision operator is a rotation operator causes the orientation of each N-MPCD particle to change as ui over the time step δt.The collision operator and resulting rerotation can be decomposed into two contributions ui = δu MS i /δt + δu J i /δt from the thermostat (δu MS i /δt) and from the flows (δu J i /δt).

A local thermostat:
The local thermostat draws random orientations from the Maier-Saupe distribution centered about n c 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 tensor 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 2E c (t T .The director dynamics are coupled back to the velocity field through the change in angular momentum δL c (t) = −γ R Nc i u i (t) × ui , 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.FIG.12: The types of wavy plane walls implemented in this study.
pendent of z and take the form The plane is normalised to standardise the value of the frequency B b,1 .This choice then fixes S b,2 (r) to be .
(B4) As 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 SI C.
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.

Boundary Rules
The b th surface boundary exists at S b (r) = 0, thus a particle is defined as being within the simulation domain if S(r i ) ≥ 0 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 ν b and perpendicular t b to the boundary (SI 3).
Keeping in mind the goal of creating periodic boundary conditions, the position r i (t) can be transformed by shifts in the normal D b,ν and/or tangential D b,t directions.After crossing the surface b, particle i's position is transformed as r i → r i + D b,ν ν b + D b,t t b .Likewise, the particle's velocity may be transformed.The velocity is altered by multiplicative parameters M b,ν and M b,t , such that upon colliding with the surface the velocity becomes 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 M b,ν = M b,t = µ b,ν = µ b,t = 1, while D b,t = 0 and D b,ν = L j , where L j is the system size in the direction j.
Solid walls are boundaries that do not apply a shift (D b,t = D b,ν = 0) but do transform the direction of the velocity.For perfect-slip boundaries, the particles are reflected and M b,ν = −1 and M b,t = 1; whereas, for no-slip conditions bounce-back rules are applied with M b,ν = M b,t = −1.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.

Strong Anchoring
For solid boundaries, the above boundary rules do not fully enforce perfect no-slip [82].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 [83].In order to generate strong anchoring conditions, all the particles in cells that intersect S b are subjected to the boundary rules described in § 2. Because the entire cell now essentially has the anchoring condition applied, this generates stronger anchoring.

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, S b,0 (r) → S b,0 (r, t).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 [4].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.

(C6)
For all surfaces, the derivative of S b,0 (r) is given by ∂S b,0 ∂r i = p b,i A b,i [A b,i (x i − q b,i )] p b,i −1 .(C7) However, the derivatives of S b,1 (r) and S b,2 (r) vary depending on the surface type.

Planes
The normals for wavy planes can be found with the derivatives of S b,1 A b,z .(C13b)

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 C ′ = CA b,z , the derivatives are ∂S b,1 ∂r i = A b,k A b,l C ′ ×      0 for r i = r j −(r l − q b,l ) for r i = r k r k − q b,k for r i = r l Flowing isotropic fluids at low Reynolds number obey plane Poiseuille flow with parabolic flow profiles [100], which is a direct consequence of the the non-slip condition at the walls of the channel.Though the coupling between flow profile, local director and anisotropic dissipation can induce complications [101], N-MPCD simulates linearized fluctuating nematodynamics with isotropic viscous dissipation [79].
The director can discontinuously transition from a homogenous alignment to deformed states [102].For strong homeotropic anchoring and sufficiently small pressure gradients, the director remains uniformly homeotropic across the channel; however, for sufficent flow the director field bows into the Bowser state [102], in which the field curves slightly in the direction of the flow (Fig. 15a).

FIG. 3 :
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 0 .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 0 = 2.5a; wall wavelength λ = 20a; and average channel height ⟨h⟩ x = 20a.
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 0 = 2.5a; wall wavelength λ = 20a; and average channel height ⟨h⟩ x = 20a.

FIG. 5 :
FIG. 5: Fluid velocity v f (r) field for a channel with wavy walls.(a) Time-averaged velocity field.Colour maps the x-component of the velocity in units a/τ , and the arrows represent the velocity vectors.Wall amplitude B 0 = 2.5a; wall wavelength λ = 20a; average channel height ⟨h⟩ x = 20a; and pressure gradient ∇P = 0.02ma 1−d τ −2 x; U = 10k B T ; 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 B T and dashed lines for isotropic phase flow.

FIG. 7 :
FIG. 7: Colloid trajectories for different wall amplitudes B 0 with homeotropic anchoring on the colloid and planar anchoring on the walls.Colours specify the amplitude B 0 .Black sinusoids on the right and left present undulations of the wavy walls.
FIG. 8: The probability density function (PDF) of finding the colloid with homeotropic anchoring on the colloid at a given position x for three different amplitudes B 0 = {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.
FIG. 9: Homeotropic colloid velocities at different positions x for two different amplitudes B 0 with planar anchoring on the wavy walls.Magenta line shows the estimated speed from the probability density.
FIG. 10: Homeotropic colloid velocity measures with varying amplitude B 0 of the wavy boundaries with planar and homeotropic anchoring.
A6) is calculated.The largest eigenvalue of Q c 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 (a) Egg carton plane wall.(b) Corrugated plane wall.
where Π b = ν b ⊗ ν b and and T b = 1 − Π b are the normal and tangent projection operators.Lastly, the orientation is also multiplicatively transformed

5 . 2 ∂S
Spheres and EllipsoidsThe surface functions for spheres and ellipsoids include arctan functions, which make the derivatives of S b,1∂S b,1 ∂z = 0 ∂S b,1 ∂x = − A b,x A b,y B (y − q y ) (C11a) ∂S b,1 ∂y = A b,x A b,y B (x − q b,x )and the derivatives of S b,(x − q b,x )(z − q b,z ) (y − q b,y)(z − q b,z ) C √ B ,where we defineB = [A b,x (x − q b,x )] 2 + [A b,y (y − q b,y )] 2(C13a) C = S b,0 + R p b,r b for r i = r k or r l .(C15)AppendixD: Flow through Plane Channels

FIG. 16 :
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.
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 0 /L = 1.The nematic case has homeotropic anchoring on the colloid and planar anchoring on the walls.