Controllable particle migration in liquid crystal flows

Abstract

We observe novel positional control of a colloidal particle in microchannel flow of a nematic liquid crystal. Lattice Boltzmann simulations show multiple equilibrium particle positions, the existence and position of which are tunable using the driving pressure, in direct contrast to the classical Segré–Silberberg effect in isotropic liquids. In addition, particle migration in nematic flow occurs an order of magnitude faster. These new equilibria are determined through a balance of elastic forces, hydrodynamic lift and drag as well as order-flow interactions through the defect structure around the particle.

---

1 Introduction

Particle-laden flows are at the heart of many scientific and engineering applications, ranging from the dispersion of pollutants and rain formation in the atmosphere to pharmaceutical aerosols and industrial spray applications to geological sedimentation and combustion processes, and even planetary formation. Such flows consist of a continuous liquid host phase and dispersed particles.

At low particle concentrations these systems are generally modelled by the Navier–Stokes equations, with suitable modifications to account for the momentum exchange between the host liquid and dispersed particles.¹ In spite of their relative simplicity, particle-laden flows show non-trivial behaviour, such as the preferential migration of the particles to specific regions. Segré and Silberberg^2,3 were among the first to report this phenomenon. In capillary flow of isotropic fluids, they observed that neutrally buoyant dispersed particles migrate to a single radial position approximately 0.6 radii from the capillary axis, while being advected downstream. This effect is generally attributed to a balance between a shear-induced inertial lift force that drives the particles towards the walls and regions of higher velocity gradient, and a pressure-induced force that increases as the particles move closer to the walls and acts in the opposite direction. A qualitative explanation,^4,5 a subsequent quantitative analysis^6,7 and then an experimental extension to higher Reynolds numbers⁸ and different tube geometries⁹ have already been provided.

Much less is known about the dynamics and migration of particles in liquid crystalline host phases, whose internal structure offers additional possibilities in terms of elastic and order-driven control. Theoretical studies have considered the Stokes drag of a particle in a nematic host¹⁰ or particle aggregates,^11,12 and simulation studies have investigated flow-order tensor interactions around colloidal particles¹³ for fixed or simply advected particles. Experimental studies of nematics in microfluidic geometries have been concerned with pure phases.^14–21 However, recent experimental research has also considered particle-laden nematic flows which contain defects and/or obstacles,²² or include electro-hydrodynamic²³ or electro-osmotic flow effects.^24–26 Further studies investigated the microfluidic flow of liquid crystals in more exotic settings, for instance the sedimentation of discoidal particles in cholesteric finger textures²⁷ or modified free energy landscapes with alternating splay and bend distortions that were introduced through wavy walls.²⁸ More recently, active nematics were studied in a network of connected annular microfluidic channels.²⁹

Here, we report results of 3D simulations that demonstrate novel positional control of a colloidal particle in microchannel flow of a nematic liquid crystal, which serves as a prototype for anisotropic fluids with internal order structure. As well as multistability, we show that the equilibrium particle position in the dynamic system with respect to the channel centre or walls is tunable through the applied pressure gradient.

2 Theory

2.1 Landau–de Gennes free energy

The local order of the liquid crystal is described by a traceless and symmetric second-order tensor Q(r,t). Its largest eigenvalue q < 2/3 is referred to as the scalar order parameter and provides a measure of the liquid crystalline order at a certain position and time. The eigenvector, d, associated with q is called the director and describes the average orientation of the liquid crystal molecules at a certain position and time.

In equilibrium, the liquid crystal order is determined through minimisation of its free energy, commonly described by the Landau–de Gennes free energy functional

which includes the volume contribution f = f_b + f_g, that itself consists of a bulk contribution f_b and a gradient contribution f_g, and a surface contribution f_s. The bulk free energy density is given bywhere we use the Einstein summation convention, which implies that Greek indices that appear twice are summed over. A₀ is a constant that sets the overall energy scale and the parameter γ controls the temperature difference from the isotropic–nematic transition, and is related to a reduced temperature τ by

For γ > 3 the ordered, nematic state is the equilibrium phase, whereas for 2.7 ≤ γ ≤ 3 the nematic state is metastable. For γ < 2.7 the isotropic state is the equilibrium phase.

The gradient free energy density f_g contains the contributions of splay, bend, splay and twist deformations of the director field,

where ∂_α = ∂/∂x_α and ε_αμν is the Levi–Civita symbol in three dimensions. In principle, the elastic constants κ₀, for splay and bend deformations, and κ₁, for twist deformations, can be different. However, in our simulations we use the one elastic constant approximation κ₀ = κ₁.

The director is assumed to have a preferred normal orientation to the wall surfaces and to the surface of the colloidal particle, known as a homeotropic weak anchoring, and is described using a surface free energy term

where w is the surface anchoring strength with values w_wall and w_part at the wall and particle surfaces, respectively, and the preferred orientation Q⁰_αβ is assumed uniaxial and is given bywhere n is the surface unit normal, δ_αβ is the Kronecker delta and S₀ is the preferred surface scalar order parameter given by

2.2 Beris–Edwards model

The time evolution of Q_αβ is governed by the Beris–Edwards equation³⁰

∂_tQ_αβ + ∂_γ(u_γQ_αβ) + S_αβ(W,Q) = ΓH_αβ, (8)

where ∂_t = ∂/∂t, u is the flow velocity, H is the molecular field, Γ is a mobility parameter, and S(W,Q) denotes the response to shear and W is the velocity gradient tensor. The shear term is given bywhere and are the symmetric and antisymmetric contributions to the velocity gradient tensor, respectively, and ξ is the so-called flow-alignment parameter, a material constant representing an effective molecular aspect ratio which determines whether the liquid crystal molecules are in a flow-aligned state at the Leslie angle or tumbling state. The molecular field H is defined as a functional derivative of the free energy functional with respect to the order parameter,

The second term in eqn (10) involving the trace ensures tracelessness of the tensor order parameter as it evolves obeying eqn (8). This leads to the following molecular field:

The governing equations of hydrodynamic motion are the equation of mass conservation, also known as the continuity equation, and the Navier–Stokes equation that describes the conservation of linear momentum. In tensor notation they read

∂_tρ + ∂_α(ρu_α) = 0 (12)

and

∂_t(ρu_α) = ∂_βΠ^(LC)_αβ + ∂_βΠ^(HD)_αβ, (13)

respectively. Eqn (12) relates the local rate of change of the density ρ to the advection of mass by the fluid velocity u_α. Eqn (13) is Newton's second law of momentum change for the fluid and involves the thermotropic stress tensor Π^(LC)_αβ and the hydrodynamic stress tensor Π^(HD)_αβ. The thermotropic stress arises due to the liquid crystal and is given by

In eqn (14), σ_αβ and τ_αβ are the symmetric and antisymmetric stress contributions, respectively, defined as

where is the isotropic pressure, and

τ_αβ = Q_αμH_μβ − H_αμQ_μβ. (16)

The final term in eqn (14) is expanded as

The hydrodynamic stress tensor is defined as

Π^(HD)_αβ = −pδ_αβ − ρu_αu_β + η(∂_βu_α + ∂_αu_β) + ζ∂_μu_μδ_αβ, (18)

where η and ζ are the dynamic and bulk viscosity, respectively. The pressure p is related to the density via an ideal gas equation of state as p = c_s²ρ with c_s as lattice speed of sound as is standard in lattice Boltzmann. The last term vanishes in incompressible fluids as eqn (12) becomes ∂_αu_α = 0. No-slip and no-penetration boundary conditions are applied on the walls and particle surfaces, and the boundary conditions for Q are found from the minimisation of the free energy³¹where Q_αβ,γ = ∂Q_αβ/∂x_γ.

The total force on the particle consists of a thermotropic contribution F^(LC) and a hydrodynamic contribution F^(HD). The thermotropic contribution is the integral of the gradient of the stress tensor Π^(LC) in eqn (14) over the surface S of the particle. This can be separated into contributions from the bulk and gradient free energy:

Splitting the molecular field H in eqn (11) into terms that contain only bulk and gradient contributions,

H^(g)_αβ = κ₀∂_α∂_μQ_μβ + κ₁∂_μ(∂_μQ_αβ − ∂_αQ_μβ), (22)

we obtain together with eqn (15)–(17) for the bulk contributionand for the gradient contribution

The hydrodynamic contribution to the force on the particle is also given through the surface integral of the gradient of the hydrodynamic stress tensor Π^(HD) over the particle surface S as

See the following Section 3 on how these forces are evaluated.

3 Simulation method

3.1 Simulation setup

Fig. 1 shows a sketch of the three dimensional computational geometry, which consists of a channel of dimensions L_x × L_y × L_z = 128 × 64 × 256 lattice sites.

Fig. 1 Sketch of the computational geometry: we apply no-slip boundary conditions and homeotropic anchoring conditions at the walls in x-direction and periodic boundary conditions at the y- and z-boundaries.

We used a hybrid lattice Boltzmann scheme³² that treats the dynamics of the Q-tensor order parameter with a finite-difference scheme and applies the lattice Boltzmann method to the hydrodynamic variables. Parallel, solid walls are positioned at x = 0 and x = L_x, whereas periodic boundary conditions are applied in y- and z-direction with the z-boundaries acting as inlet and outlet of the microchannel. The pressure gradient Ψ = Δp/L_z is applied in z-direction as body force density that acts on all fluid sites besides the forces that arise from the thermotropic and hydrodynamic stresses in eqn (14) and (18), respectively. Our approach uses Guo forcing^33,34 to apply all forces on the fluid, which removes undesirable artefacts that may enter the continuity and momentum equation due to the time discretisation.

The colloidal particles are discretised as solid, mobile particles with a radius of R = 7.2 or 9.6 lattice sites. All of our results use the radius R = 9.6 apart from those presented in Fig. 10. As mentioned above, on the walls and particle surfaces no-slip and no-penetration boundary conditions are imposed through the bounce-back on links scheme.^34,35 The surface free energy in eqn (5) invokes a weak homeotropic anchoring condition with a preferred orientation of the director normal to the surfaces. The thermotropic force on the particle F^(LC) in eqn (20) is integrated using a finite-difference scheme in Q and its gradients. The hydrodynamic force F^(HD) in eqn (25) is integrated in a similar way using the hydrodynamic stress tensor Π^(HD), which is directly accessible in the lattice Boltzmann method through second order moments of the non-equilibrium distributions. The total force F = F^(LC) + F^(HD) is fed into a molecular dynamics algorithm to integrate the motion of the particles.

Several technical limitations of our model should be noted. While the centre of mass of the particle is integrated off-grid according to Newton's equation, the particle itself is discretised using a stair-case geometry. This requires remapping of the particle onto the lattice as the particle moves. Consequently, this can entail spikes in the force at single iteration steps, although these quickly average out over a few iteration steps with no detrimental effect on the trajectories. The pressure obeys an ideal gas equation of state and is directly related to the density via p = c_s²ρ, while the effect of the constant pressure gradient is modelled through an additional body force density on the fluid. Both treatments are common in the lattice Boltzmann methodology and allow for an accurate modelling of a weakly compressible fluid, but the assumption of a constant pressure gradient represents obviously a simplification over the real situation. Thermal fluctuation have not been included as our simulations were carried out at a point in the phase diagram that is deep in the ordered state well away from the isotropic–nematic transition line where elastic forces from the anchoring of the liquid crystal dominate over thermal forces by orders of magnitude. Our approach uses the one-constant approximation to simplify elasticity where the elastic constants for splay, bend and twist deformations carry the same value. This approximation is commonly used as first approach and does not compromise our results qualitatively. However, relaxing this approximation could lead to quantitative differences, and potentially also richer phenomenology. The Beris–Edwards model uses a simplified approach to viscosities compared to the Ericksen–Leslie theory, which has six viscosity coefficients α₁,…,α₆³⁶ (only five are independent as the Parodi relation applies³⁷). The viscosities in the Beris–Edwards model are implicitly given through the isotropic dynamic shear viscosities η, the rotational diffusion constant Γ, the flow alignment parameter ξ and the scalar order parameter S₀. They can be directly related to the Ericksen–Leslie viscosities α₁,…,α₆,³² but parameterise only a subset of possible values.

All simulations were run with Ludwig, our lattice Boltzmann code for complex fluids, version 0.15.0. Typical simulations ran for 8 × 10⁵ iteration steps at various pressure gradients, each of which took approximately 16 hours to complete using a hybrid Message Passing Interface/OpenMP parallelisation with 4 MPI-tasks each running on 20 OpenMP threads. Model parameters used for the simulations is provided in Table 1. More information about the specific implementation used in this work can be also found in the Ludwig code repository³⁸ and related literature.^39,40

Table 1 Overview of simulation parameters

Bulk energy scale	A 0	0.01
Inverse temperature	γ	3.1
Elastic constants	κ 0, κ1	0.01
Wall anchoring strength	w wall	0.02
Particle anchoring strength	w part	0.01
Flow alignment parameter	ξ	0.7
Mobility parameter	Γ	0.5
Density	ρ	1.0
Dynamic viscosity	η	5/6
Bulk viscosity	ζ	5/6
Particle radius	R	7.2, 9.6

---

3.2 Parameter mapping to physical units

In order to map our simulation units to physical units we need to calibrate the units of length, time and pressure. To this end, we relate the lattice spacing Δx, the algorithmic time step Δt and the reference pressure p*, which are all unity in lattice Boltzmann units (LBU), to their values in SI units.

The calibration of the length scale is straightforward as it is simply set by considering the diameter D of the colloidal particle. Assuming the largest radius particle that we consider corresponds to a relatively small diameter of D = 0.2 μm in SI units, results in a LBU of length Δx ≙ 10⁻⁸ m = 10 nm in SI units. This length scale allows for an accurate resolution of the liquid-crystalline order structure and flow field around the particle, while keeping the necessary computational resource relatively low.

To obtain a pressure scale, we use the measurements of the Landau–de Gennes parameters⁴¹ (see Appendix D therein), which suggest

for some liquid crystals. Using A₀ = 0.01 and γ = 3.1 in our simulations leads to a reference pressure of p* = 1 LBU ≙ 10⁸ Pa in SI units.

For the timescale calibration we use the following formula, which relates the rotational viscosity γ₁ of the director, to the equilibrium scalar order parameter q and the order parameter mobility Γ:

We use Γ = 0.5 and bulk energy density parameters that give q ≈ 0.5 since it is assumed that the system is well within the nematic phase. Therefore, the rotational viscosity γ₁ = 1 LBU. Typical values for liquid crystals in SI units are γ₁ = 0.1 Pa s.³⁶ Together with 1 Pa equating to a pressure of 10⁻⁸ in LBU, we obtain for the algorithmic time step Δt ≙ 10⁻⁹ s = 1 ns.

The Reynolds number gives the ratio of inertial to viscous forces, and in our approach this can be estimated on the basis of typical velocities of colloidal particles as they are advected with the flow. For instance, in a later simulation (Fig. 6) with Ψ = 9.6 × 10⁻⁶ and R = 9.6, we observe that a colloidal particle that migrates preferentially to the attractor region at x ≃ 47 travels around 1.6 × 10⁴Δx in the z-direction during the 8 × 10⁵Δt (all simulations in Fig. 3, 6 and 10 were run for 8 × 10⁵ algorithmic steps). This leads to a velocity v ≃ 0.02 in LBU. With the values for density ρ and dynamic viscosity η as specified in Table 1 and a diameter of a colloidal particle D = 20Δx as typical length scale Λ, this results in a Reynolds number

for this pressure difference, and indicates a flow regime where viscous forces are larger than, but comparable to inertial forces.

The Ericksen number gives a ratio of viscous to elastic forces and, with the above values of the flow velocity v, dynamic viscosity η, characteristic length scale Λ and either elastic constant κ₀ or κ₁ (see Table 1), we obtain

which defines a flow regime where the director field is strongly influenced by the flow.

4 Results and discussion

Before addressing preferential particle migration in nematic host phases, we begin with presenting some general results for contextualisation and later reference.

In the quiescent state the homeotropic wall anchoring induces a nematic order with a director orientation parallel to the wall normals (z-direction in Fig. 1). When the flow is switched on, flow-alignment of the director to the appropriate Leslie angle takes place. The nematic liquid crystal can adopt two possible conformations, bend or splay, as shown in Fig. 2.

Fig. 2 Bend or H-state (left) and splay or V-state (right) at the centre of the channel.

In both states the director flow-aligns to a positive Leslie angle in the lower half of the channel (where there is positive shear) and to a negative Leslie angle in the upper half of the channel (where there is negative shear). The bend state (sometimes called the H-state) and splay state (sometimes called the V-state) differ in how the director rotates between the positive and negative Leslie angles at the centre of the channel. In the bend state the director at the centre is perpendicular to the walls and in the splay state the director at the centre is parallel to the walls.

In the absence of any liquid crystalline order, i.e. either in a classical Newtonian fluid or in a liquid crystal host phase at temperatures above the isotropic–nematic transition point, preferential migration of colloidal particles occurs to an x-position between the centre and wall of the channel. This is referred to as the Segré–Silberberg effect and shown in Fig. 3 through trajectories for all starting positions x_s at all applied pressure gradients Ψ. Two aspects of the Segré–Silberberg effect are important to understand when comparing to our results further below. Firstly, except for particles that start exactly at the centre of the channel, migration is to a single equilibrium x-position over time. Secondly, the value of x does not vary significantly at low Reynolds numbers. As seen in Fig. 3, in simulations we see exactly these aspects of the Segré–Silberberg effect: all particles migrate to an equilibrium at x-position of approximately 38 regardless of their starting position or applied pressure.

Fig. 3 Particle migration in an isotropic fluid due to the Segré–Silberberg effect. Shown are particle position x, across the channel gap, versus the particle position z, along the channel, for a variety of starting positions (color coded) and pressure gradients in the range 1.25 × 10⁻⁶ ≤ Ψ ≤ 1.75 × 10⁻⁵. All quantities are given in LBU.

It should be noted that the Segré–Silberberg equilibrium position seen in Fig. 3 is not at the well-known 0.6 tube radii (or channel half-widths in the planar case) from the centre line, but at approximately 0.4 channel half-widths. This shift of the equilibrium position towards the channel centre has been observed before in simulations^42,43 and can be attributed to our particle confinement ratio of R/L_x = 19.2/128 = 0.15. In contrast, the analytical results^6,7 were obtained for zero confinement ratio, i.e. for particles that are negligibly small compared to the tube diameter or gap width. Our study uses Reynolds numbers that are around two orders of magnitude smaller than those in previous simulation studies,^42,43 which is known to result in equilibrium positions that are again closer to the channel centre.

Using the model described in Section 2 we studied the migration behaviour of a single particle at different pressure gradients Ψ (= Δp/L_z) and start positions x_s.

Fig. 4 shows particle trajectories x(z) for starting positions in the lower portion of the channel (the channel centre is at x = 64). For low pressure gradients Ψ = 7.5 × 10⁻⁶ (red lines in Fig. 4) the particle migrates either to the channel walls (for starting position x_s ≤ 40) or towards the centre of the channel (for x_s ≥ 44). Particles that start at the centre stay in that position, while those inserted further away from the centre show a tendency to overshoot and remain in an off-centre position. At slightly larger pressure gradients Ψ = 8.125 × 10⁻⁶, 8.75 × 10⁻⁶ (blue and green lines in Fig. 4, respectively), this bi-stability remains, but with the division between migration towards the wall or centre now at around x_s ≃ 38. Some trajectories exhibit the onset of a pull-back behaviour from overshot off-centre positions towards the centre, in particular for the larger pressure gradient of Ψ = 8.75 × 10⁻⁶. It should be noted that for low pressure gradients, Ψ ≤ 5 × 10⁻⁶, we also observe off-centre equilibrium positions in form of a weak attractor that moves to the channel centre with increasing pressure gradient.

Fig. 4 Particle trajectories for particle radius R = 9.6 at sub-critical pressure gradients Ψ = 7.5 × 10⁻⁶, 8.125 × 10⁻⁶, 8.75 × 10⁻⁶, showing particle x-positions across the channel gap versus the z-distance travelled in the flow direction for various initial positions x_s. The particle Ericksen numbers Er (particle Reynolds numbers Re) are from front to back Er = 25.62 (Re = 0.37), Er = 27.91 (Re = 0.40), Er = 30.20 (Re = 0.43). All quantities are given in LBU.

A further increase in the pressure gradient leads to a sudden onset of a new type of preferential migration behaviour. In Fig. 5 we see for Ψ = 9.125 × 10⁻⁶ that a pronounced trajectory kink emerges at x ≃ 50–53, the precursor of which is also visible for lower pressure gradients in Fig. 4. For the pressure gradients Ψ = 9.6 × 10⁻⁶ (blue lines in Fig. 5) this kink transitions into a third, emergent particle attractor, in addition to the channel centre and wall. At Ψ = 9.6 × 10⁻⁶ this emergent attractor is located at x ≃ 48 and trajectories with initial positions 40 < x_s < 50 migrate towards the attractor. Increasing the pressure gradient further, leads to a movement of the emergent attractor towards the wall, and for Ψ = 1.125 × 10⁻⁵ the attractor position is x ≃ 40.

Fig. 5 Particle trajectories for particle radius R = 9.6 at pressure gradients Ψ = 9.125 × 10⁻⁶, 9.6 × 10⁻⁶, 1.125 × 10⁻⁵, showing particle x-positions across the channel gap versus the z-distance travelled in the flow direction for various initial positions x_s. The emergence of the particle attractor is shown for Ψ = 9.6 × 10⁻⁶, 1.125 × 10⁻⁵. The particle Ericksen numbers Er (particle Reynolds numbers Re) are from front to back Er = 31.58 (Re = 0.45), Er = 33.33 (Re = 0.48), Er = 39.57 (Re = 0.57). All quantities are given in LBU.

Fig. 6 shows the complete set of trajectories for a particle with radius R = 9.6 and pressure gradients ranging from Ψ = 1.25 × 10⁻⁶ to 1.75 × 10⁻⁵. The first row of Fig. 6 shows results for low pressure gradients Ψ = 1.25 × 10⁻⁶, 2.5 × 10⁻⁶, 5 × 10⁻⁶, indicating typical behaviour at low pressure gradients of migration toward the wall or toward the weak attractor region that moves closer to the channel centre with increasing pressure gradient. In all cases the director remains in the bend state, as indicated by blue lines. The second row of Fig. 6 shows trajectories for intermediate pressure gradients Ψ = 7.5 × 10⁻⁶, 8.75 × 10⁻⁶, 9.125 × 10⁻⁶, for which we can observe trajectories with overshooting and pull-back behaviour as well as a transition from the bend (blue lines) to the splay (red lines) state, through an intermediate state (black lines). For the lowest of these intermediate pressure gradients, Ψ = 7.5 × 10⁻⁶, overshooting can be seen.

Fig. 6 Particle trajectories in a nematic liquid crystal host phase for particle size of R = 9.6 and applied pressure gradients ranging from Ψ = 1.25 × 10⁻⁶ to 1.75 × 10⁻⁵. Blue lines indicate that the director structure is in a bend state, whereas red lines indicate that the director has transitioned to the splay state. Black lines mark a transition state between bend and splay. The particle Ericksen numbers Er and particle Reynolds numbers Re are given in each sub-plot. All quantities are given in LBU.

Overshooting followed by pull-back to the centre is seen for Ψ = 8.75 × 10⁻⁶, but the overshoot has disappeared at the higher pressure gradient of Ψ = 9.125 × 10⁻⁶. In all cases of overshoot and pull-back there is a bend to splay transition. The precursor of the emergent attractor is clearly seen in the trajectories for the 9.125 × 10⁻⁶ case, the right most image of the second row of Fig. 6. For the highest intermediate pressure gradient, Ψ = 9.6 × 10⁻⁶, the third attractor state emerges. The last two rows of Fig. 6 cover the higher range of pressure gradients, for which the emergent attractor state is possible, albeit for sufficiently high pressure gradients. It is clear that the position of the emergent attractor state moves towards the wall with increasing pressure gradient. To further check the consistency of our findings, we also ran a number of additional simulations with R = 9.6 at a lower shear and bulk viscosity η = ζ = 1/6, which led to the same equilibrium position and confirmed these results.

The main plot in Fig. 7 summaries these results, depicting a phase diagram of the equilibrium position as a function of the starting position x_s and the applied pressure gradient Ψ. For particles initially close to the walls, migration towards the walls occurs for all pressure gradients. Similarly, particles initially close to the centre of the channel migrate to the centre of the channel. Initial particle positions that are around one or two particle radii away from the centre give rise to a much more complicated behaviour: at low pressure gradients the particle migrates to an off-centre weak attractor position; for a narrow range of pressure gradients around Ψ ≃ 9 × 10⁻⁶ centre positions are again favoured; and above a critical pressure gradient Ψ ≃ 9.5 × 10⁻⁶ the phase diagram is increasingly dominated by the emergent attractor states.

Fig. 7 Phase diagram of the preferential migration of the colloidal particle in a nematic host phase. Coloured regions show the equilibrium particle position (wall, weak attractor, emergent attractor, centre) as a function of the initial x-position of the particle and the applied pressure gradient Ψ. Equilibrium director bend states are marked with dots, whereas equilibrium splay states are shown with crosses. Subplots (a)–(e) show the scalar order parameter (green isosurface showing the low order region) and director field (short coloured lines) around the particle for typical equilibrium states: (a) bend state, centre position; (b) bend state, wall position; (c) bend state, emergent attractor position; (d) splay state, centre position; (e) splay state, weak attractor position. The centre of the channel is marked by the horizontal blue line. See also Movies 1–5 in the ESI† for the dynamic evolution to the steady states (a), (b), (c), (d) and (e), respectively.

Also shown in Fig. 7, through the subplots (a)–(e), are the director and scalar order parameter structures within the channel. Previous theoretical and experimental work has shown that the director structure within channel flow of a nematic can exhibit either predominately bend or splay deformation^17,19 for low and high flow speeds, respectively. Both states are observed in our simulations, and both exhibit flow alignment at the appropriate Leslie angle, positive (negative) angles in the lower (upper) half of the channel where the shear gradient is positive (negative). As mentioned above, the states differ in their transition between positive and negative Leslie angle at the centre of the channel, with the bend state exhibiting director alignment along x at the centre of the channel (for instance in Fig. 7(a)), and the splay state exhibiting director alignment along z at the centre (Fig. 7(d)). For the pressure gradients we consider here, the bend state would normally be maintained. However, we observe a novel particle-induced mechanism for switching from the bend to splay state (see Fig. 2) with the migration of the particle to the centre being the initiating event for a range of pressure gradients (see Fig. 7, indicated by crosses, and ESI† Movies 4 and 5).

For wall equilibrium positions a bend state occurs with a tilted, but otherwise regular, Saturn ring defect around the particle (see Fig. 7(b)). For centre equilibrium positions a bend state occurs for Ψ < 1.2 × 10⁻⁵ and a splay state for Ψ > 1.2 × 10⁻⁵. For an initial particle position away from the centre, the bend state occurs for the weak attractor, centre and emergent attractor equilibrium positions, for both low and high pressure gradients. However, at intermediate pressure gradients 6 × 10⁻⁶ ≤ Ψ ≤ 8.5 × 10⁻⁶ we observe a transition to the splay state (Fig. 7(d) and (e)). We note that the transition from the bend to the splay state does not have a determining effect on whether the particles migrate to centre or off-centre positions as we see the same behaviour for initial positions close by and/or at lower pressure gradients. However, it affects the nature of the Saturn ring defect around the particle: for bend equilibrium states the Saturn ring defect is approximately horizontal (Fig. 7(a) and (c)) and for splay equilibrium states the Saturn ring defect is approximately vertical (Fig. 7(d) and (e)). For bend equilibrium states, the defect structure also develops a pronounced lip- or cap-shaped region of low liquid-crystalline order at the bow of the particle (Fig. 7(a) and (c)).

In comparison to the Segré–Silberberg effect in isotropic fluids (Fig. 3) we notice several fundamental differences. Firstly, in isotropic fluids, particles migrate to a single equilibrium position, whereas in nematic host phases migration to one of multiple equilibrium positions can occur, depending on particle starting positions and the applied pressure gradient. Secondly, while for the Segré–Silberberg effect the location of the attractor state between the wall and centre depends only weakly on the flow velocity^3,7 (in fact our simulations with an isotropic host phase, see Fig. 3, show no appreciable change in equilibrium position over the entire range of applied pressure gradients), the position of the emergent attractor state in the nematic system depends much more sensitively on the imposed pressure gradient, and therefore on the Reynolds number. For instance, the Reynolds numbers for which there are emergent attractor states in Fig. 5 are Re ≃ 0.48 and Re ≃ 0.57 and, even with this small increase in Re, the attractor position moves by almost one particle radius. Finally, the preferential migration in a nematic host phase happens more than an order of magnitude faster than in isotropic fluids. For instance, for a pressure gradient Ψ = 1.125 × 10⁻⁵ (Fig. 5, green lines), where the position of the emergent attractor in the nematic host phase and the Segré–Silberberg equilibrium position in the isotropic host phase almost coincide, particles in a nematic host reach their equilibrium positions by the time they have travelled around z ≃ 5 × 10³ along the channel, and in an isotropic host phase (see Fig. 3) take, depending on the start positions x_s, at least an order of magnitude longer (in time or distance along the channel).

With regard to Ericksen numbers Er, our results indicate that the observed pattern of preferrential migration occurs in a regime where viscous forces are larger than elastic forces, and where consequently the director field is strongly affected by the flow field. At the lower end of this regime, for 3.43 ≤ Er ≤ 16.55, the weak attractor state exists, as shown in the top row of the Fig. 6. The emergent attractor state occurs for Er ≃ 30 and above. With increasing Er, viscous forces begin to dominate over elastic forces. This means the flow behaviour at higher Ericksen numbers is more akin to that of an isotropic fluid, which shows the classical Segré–Silberberg effect. However, the fact that a strongly flow-aligned liquid crystal forms the host phase leads always to certain qualitative differences. For instance, the trajectories shown in the two bottom rows of Fig. 6 feature the emergent attractor state at 33.33 ≤ Er ≤ 62.12, do show some similarities to the classic Segré–Silberberg effect. However, the movement of the emergent attractor region towards the walls with increasing Ericksen number, the attraction to the walls, or the existence of stable trajectories at the channel centre are all features that arise due to the anisotropic nature of the host phase and the interaction of flow-aligned director field with the defect structure around the particle, as is visible in the snapshots shown in Fig. 7(a)–(e).

In order to investigate the presence of this emergent attractor, we consider the total force and the contributions to the total force on the particle in the steady state for particles that have equilibrated (see Fig. 8). An analysis of the individual force contributions shows that, for a centre equilibrium position as in Fig. 7(d) all three force components vanish, as expected from symmetry, and for a wall equilibrium position as in Fig. 7(b) all three force components are negative, forcing the particle to remain at the wall. For particle migration to an emergent attractor state as in Fig. 7(c), our simulation show that the particle feels a force towards the channel centre from the gradient, i.e. elastic, terms and forces towards the wall from both the bulk term and the hydrodynamic force component. This is also the case for particle migration to the centre with an asymmetric defect as in Fig. 7(a) and for particle migration to the weak attractor as in Fig. 7(e), although the individual force contributions are significantly smaller. Therefore, a delicate force balance exists between these relatively large contributions leading to a zero total force at equilibrium.

Fig. 8 Total force and the contributions to the total force on the particle in the steady state for particles that have equilibrated to lie at (a) the centre with asymmetric defect, (b) at the wall, (c) at the emergent attractor, (d) at the centre with symmetric defect, and (e) at the weak off-centre attractor, for particle radius R = 9.6. The error bars indicate one standard deviation of the force data during a run over 10⁴ algorithmic steps. All quantities are given in LBU.

Fig. 9 shows the time evolution of the different contributions to the force on the particle during its approach to an emergent attractor state, for pressure gradient Ψ = 1.125 × 10⁻⁵ and particle radius R = 9.6. For these values of Ψ and R the attractor state is located, approximately, at the Segré–Silberberg equilibrium x-position (namely that seen in Fig. 3). Trajectories for three starting positions are shown, for x = 34.5 (light grey), x = 39.5 (medium grey) and x = 59.5 (dark grey), which cover a range of migration patterns to the attractor from below and above. Also shown are the contributions to the total force on the particle which arise from gradient (red), bulk (purple) and hydrodynamic (blue) terms with the hue (light, medium, dark) corresponding to the equivalent starting position (light grey, medium grey, dark grey).

Fig. 9 Time evolution of force contributions in the emergent attractor state with pressure gradient Ψ = 1.125 × 10⁻⁵. Trajectories for three starting positions are shown, for x = 34.5 (light grey), x = 39.5 (medium grey) and x = 59.5 (dark grey). Also shown are the contributions to the total force on the particle which arise from gradient (red), bulk (purple) and hydrodynamic (blue) terms with the hue (light, medium, dark) corresponding to the equivalent starting position (light grey, medium grey, dark grey). All quantities are given in LBU.

From Fig. 9 we observe initial transient behaviour before the particle reaches the same equilibrium position for all three starting positions. For all starting positions and for all time, the gradient contribution is positive, so that the elastic forces always act to move the particle towards the channel centre. For both bulk and hydrodynamic force contributions are negative for all starting positions and (almost all) time, so that they act to move the particle towards the wall. The exception is that for starting position x = 34.5 the bulk force contribution is positive for a very short initial period, i.e. for very small z-displacements of the particle from the initial position. Interestingly, all three force contributions are almost balanced, with a total force of zero, for all time. Variations in the total force are indiscernible on the same force scale used in Fig. 9.

This situation is in contrast to the Segré–Silberberg effect, in which the inertial component of the hydrodynamic force acts to move the particle across the shear gradient towards the wall, while the increased pressure caused by the particle moving towards the wall leads to a force acting to move the particle towards the centre of the channel. At the Segré–Silberberg equilibrium position, the total hydrodynamic force vanishes. In a nematic host material, the gradient terms act with the inertial hydrodynamic forces to move the particle towards the wall, allowing much faster migration of the particle and the appearance of an attractor state at a much smaller Reynolds number.

The presence of a particle also leads to a distorted director structure, rather than a uniform director at the Leslie angle and so increased pressure gradients (equivalent to an increased Reynolds number) can align the director around the particle more closely to the Leslie angle, thus adapting the elastic force on the particle, and therefore control of the attractor equilibrium x-position through changes in pressure gradient is possible.

Additional results for smaller particles with R = 7.2 show a similar pattern of preferential migration, and are depicted Fig. 10 as trajectories. There are quantitative differences, for instance overshooting off-centre trajectories or trajectories with overshooting and pull-back to the channel centre are much less pronounced compared to larger particles with R = 9.6. This can be attributed to smaller inertia and reduced anchoring forces, which scale down with the surface area of the particle. For smaller particles the attractor state emerges at slightly larger pressure gradient around Ψ = 1.125 × 10⁻⁵, but shows the same characteristics, in particular the fast migration of the particle to the equilibrium regions, as well as the movement of the emergent attractor position towards the walls and its prominence in the phase diagram with increasing pressure gradient. These results suggest that this new effect depends to a certain extent also on particle size, and therefore on inertia. For higher inertia we expect the emergent attractor region in the phase diagram Fig. 7 and states as in Fig. 7(c) to extend towards lower pressure gradients. The middle regions with states as in Fig. 7(d) and (e) are also likely to appear at lower pressure gradients and to grow in size. This would occur at the cost of the region with wall attraction and states as in Fig. 7(b), which is likely to have a smaller extent.

Fig. 10 Particle trajectories in a nematic liquid crystal host phase for particle size R = 7.2 and various applied pressure gradients from. Blue lines indicate that the director structure is in a bend state, red lines indicate the splay state and black sections denote the transition. The particle Ericksen numbers Er and particle Reynolds numbers Re are given in each sub-plot. All quantities are given in LBU.

5 Conclusions

In summary, we observe multiple equilibrium particle positions and a new pressure-controllable particle attractor state for a colloidal particle with nematic liquid crystalline host phase. At low pressure gradients particles migrate either to the channel centre or the walls but at higher pressure gradients a third attractor state emerges spontaneously, whose position in the channel depends sensitively on the pressure gradient. These results are in striking contrast to the classical Segré–Silberberg effect in isotropic fluids, where the equilibrium position is also reached more slowly. The discovery of these new and controllable attractor positions opens up interesting routes for tailored particle separation. While our results were obtained in pressure-driven flow, we expect them to hold as well in flux-driven flow as long as there is no significant drag between particle and fluid. However, it is likely that, as well as particle size, the confinement ratio and anchoring type and strength offer additional mechanisms to control the particle migration. These points will be addressed in future studies.

Data availability

Data underpinning this publication are openly available from the University of Strathclyde KnowledgeBase research information portal.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We acknowledge support from EPSRC under Grant No. EP/N019180/2 (O. H.) and EP/R513349/1 (M. L.). M. L. acknowledges funding from the Mac Robertson Postgraduate Travel Scholarship. We would like to thank Uroš Tkalec and Timm Krüger for very helpful discussions. For the purpose of complying with UKRI's open access policy, the author has applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising from this submission. This work used the ARCHIE-WeSt High Performance Computer (https://www.archie-west.ac.uk) based at the University of Strathclyde and the Cirrus UK National Tier-2 HPC Service (https://www.cirrus.ac.uk).

References

1. G. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, 2000.
2. G. Segre and A. Silberberg, Nature, 1961, 189, 209–210.
3. G. Segre and A. Silberberg, J. Fluid Mech., 1962, 14, 136–157.
4. P. G. Saffman, J. Fluid Mech., 1965, 22, 385–400.
5. K. Gotoh, Nature, 1970, 225, 848–850.
6. B. P. Ho and L. G. Leal, J. Fluid Mech., 1974, 65, 365–400.
7. J. A. Schonberg and E. J. Hinch, J. Fluid Mech., 1989, 203, 517–524.
8. J. P. Matas, J. F. Morris and É. Guazzelli, J. Fluid Mech., 2004, 515, 171–195.
9. D. Di Carlo, D. Irimia, R. G. Tompkins and M. Toner, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 18892–18897.
10. H. Stark and D. Ventzki, Europhys. Lett., 2002, 57, 60–66.
11. S. Mondal, I. M. Griffiths, F. Charlet and A. Majumdar, Fluids, 2018, 3, 39.
12. S. Mondal, A. Majumdar and I. M. Griffiths, J. Colloid Interface Sci., 2018, 528, 431–442.
13. T. Stieger, M. Schoen and M. G. Mazza, J. Chem. Phys., 2014, 140, 054905.
14. A. Sengupta, U. Tkalec, M. Ravnik, J. M. Yeomans, C. Bahr and S. Herminghaus, Phys. Rev. Lett., 2013, 110, 048303.
15. A. Sengupta, C. Bahr and S. Herminghaus, Soft Matter, 2013, 9, 7251–7260.
16. A. Sengupta, S. Herminghaus and C. Bahr, Liq. Cryst. Rev., 2014, 2, 73–110.
17. T. Anderson, E. Mema, L. Kondic and L. Cummings, Int. J. Non-Linear Mech., 2015, 75, 15–21.
18. L. Giomi, Ž. Kos, M. Ravnik and A. Sengupta, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E5771–E5777.
19. S. Čopar, v Kos, T. Emeršič and U. Tkalec, Nat. Commun., 2020, 11, 59.
20. P. Steffen, E. Stellamanns and A. Sengupta, Phys. Fluids, 2021, 33, 072005.
21. S. Čopar, M. Ravnik and S. Žumer, Crystals, 2021, 11, 956.
22. K. Chen, O. J. Gebhardt, R. Devendra, G. Drazer, R. D. Kamien, D. H. Reich and R. L. Leheny, Soft Matter, 2017, 14, 83–91.
23. Y. Sasaki, H. Hoshikawa, T. Seto, F. Kobayashi, V. S. R. Jampani, S. Herminghaus, C. Bahr and H. Orihara, Langmuir, 2015, 31, 3815–3819.
24. O. D. Lavrentovich, I. Lazo and O. P. Pishnyak, Nature, 2010, 467, 947–950.
25. O. D. Lavrentovich, Soft Matter, 2014, 10, 1264–1283.
26. C. Peng, T. Turiv, Y. Guo, Q.-H. Wei and O. D. Lavrentovich, Liq. Cryst., 2018, 45, 1936–1943.
27. K. Chen, L. P. Metcalf, D. P. Rivas, D. H. Reich and R. L. Leheny, Soft Matter, 2015, 11, 4189–4196.
28. Y. Luo, D. A. Beller, G. Boniello, F. Serra and K. J. Stebe, Nat. Commun., 2018, 9, 3841.
29. J. Hardouin, J. Laurent, T. Lopez-Leon, J. Ignes-Mullol and F. Sagues, Soft Matter, 2020, 16, 9230–9241.
30. A. N. Beris and B. J. Edwards, Thermodynamics of flowing systems, Oxford University Press, New York, 1994.
31. M. Škarabot, M. Ravnik, S. Žumer, U. Tkalec, I. Poberaj, D. Babič, N. Osterman and I. Musevič, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 1–8.
32. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 031921.
33. Z. Guo, C. Zheng and B. Shi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 046308.
34. T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. M. Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer International Publishing, Switzerland, 2017.
35. A. Ladd, J. Fluid Mech., 1994, 271, 285.
36. P.-G. de Gennes and J. Prost, The Physics of Liquid Crystals, Oxford University Press, New York, 1993.
37. O. Parodi, J. Phys., 1970, 31, 581–584.
38. Ludwig GitHub repository, https://github.com/ludwig-cf/ludwig.
39. J. C. Desplat, I. Pagonabarraga and P. Bladon, Comput. Phys. Commun., 2001, 134, 273–290.
40. R. Adhikari, K. Stratford, M. E. Cates and A. J. Wagner, Europhys. Lett., 2005, 71, 473–479.
41. D. C. Wright and N. D. Mermin, Rev. Mod. Phys., 1989, 61, 385–432.
42. A. Kilimnik, W. Mao and A. Alexeev, Phys. Fluids, 2011, 23, 123302.
43. I. Lashgari, M. Ardekani, I. Banerjee, A. Russom and L. Brandt, J. Fluid Mech., 2017, 819, 540–561.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2sm00707j

---