Emergent ordering of microswimmers in smectic liquid crystals

Self-propelled agents can interact in many different ways, including by perturbing a shared physical environment, e.g, hydrodynamic interactions. We study motile particles that are embedded into a smectic liquid crystal, locally distorting the smectic layer spacing. This results in interactions mediated by a smectic liquid crystal distortion field and corresponds to a form of “active smectic” liquid crystal. We identify several dynamical phases that emerge in different regimes of the smectic stiffness and particle reorientation time. We characterise these as (i) ballistic motion, (ii) clustering and (iii) collective motion where orientational order emerges even thought the system lacks explicit co-alignment. We further identify an order-to-disorder transition on the addition of angular noise.


Introduction
Collective motion has been observed in many natural systems across all biological scales in a wide range of environments, e.g. from the macroscopic: shoals of fish [1], flocks of starlings [2], and herds of sheep [3], to the microscopic: swarming krill [4], bioconvecting phytoplankton [5], and swimming bacteria [6]. There is contemporary interest in identifying novel in vitro systems that mimic this behaviour. With that in mind we study the behaviour of motile "microswimmer" particles, e.g. flagellated bacteria, embedded in a large layer-spacing smectic liquid crystal.
Smectic liquid crystals are comprised of molecules that are arranged into layers (lamellae) but are free to move within the plane of the layer with each layer behaving like a 2D liquid. The discrete layers of molecules are organised into lamellar stacks, with nearly uniform spacing between neighbouring layers [7]. Examples of such objects include the lamellar phase of surfactants [8] and di-or triblock copolymers [9,10], as well as many others in chemistry. Such phases can presumably be built from microbiologically compatible components although little research has yet explored this. We hope the present work may motivate interest in such phases as they may represent one of the most suitable experimental systems for studies of microswimmers in smectics. Here, the main design consideration is that the smectic layer spacing should be larger than or comparable to the swimmer size. This ensures that the swimmer does not destroy the smectic ordering locally. It seems likely that lyotropic copolymer lamellar phases can be synthesised with layer spacings in the micron range, only an order of magnitude larger than the lamellar thickness. This would suggest that such phases could then be compatible, e.g. with small flagellated or ciliated bacteria.
The presence of inclusions between the smectic layers, such as colloidal particles, can distort the ordering of the smectic layers [11]. The consequence of this is a distortion energy that depends on the relative positions of any particles, as well as microscopic parameters controlling the stiffness of the smectic. This results in an interaction potential between inclusions that generates smectic-mediated forces that can be either repulsive or attractive depending on the spatial configuration of particles. This motivates us to study the embedding of microswimmers in smectic phases that then 1 3 provide non-trivial smectic-mediated interactions between the particles. The goal of the current work is to explore the behaviour of such systems and any collective motion that results when there are many particles present. Whilst there has been recent interest in active nematics [12], active smectics remain relatively unexplored [13].
Our strategy is to first build a minimal physics model and then run dynamical simulations to explore its behaviour.

The smectic-mediated force
Colloidal particles, either active (motile) or passive, that are embedded between layers of a smectic liquid crystal may distort the local smectic layer spacing. This distortion can become large, e.g. when the particle size becomes comparable with the mean smectic layer spacing. This distortion has an associated energy cost that we now calculate using the Landau-de Gennes Hamiltonian [7,14,15]. This is a perturbative approach involving a small displacement, u( ) of the layers in the z-direction, normal to the plane representing their ground-state position in the absence of inclusions. On symmetry grounds two terms appear in the energy density of the smectic: A term B( z u) 2 accounts for local layer compression ( z u < 0 ) or expansion ( z u > 0 ), with B a compression modulus, and a term K(∇ 2 || u) 2 accounts for the bending of layers, with K a bending stiffness and ∇ 2 || the Laplacian parallel to the layers, in the x-y plane spanned by || = (x, y) with r �� = √ x 2 + y 2 . An additional energy density must be included to account for how the inclusions distort the smectic layer spacing, written z u . Here is the local density of colloids so that a single colloid corresponds to a region over which ∫ dV = 1 , and is a coupling constant controlling how much the layer spacing is expanded locally by the presence of a single particle. In most systems this likely depends on the size of the particle compared to the smectic layer spacing, with increasing with the particle size. The Hamiltonian then reads with the integral over the system volume (here considered infinite), dV a volume element and the prefactor 1/2 an unimportant convention. Minimising the Hamiltonian for two particles with a separation yields the interaction potential G(r) between them given by 1 K∕B a characteristic length comparable to the smectic layer spacing. A sketch showing two interacting particles is shown in Fig. 1 while the two-particle interaction potential shown in Fig. 2.
The smectic-mediated force is given by the spatial derivative of G( ) . We assume that particles (swimmers) are unable to move between layers, so we set the z-component of this force to zero on the assumption that any such force is balanced by a restoring force provided by the impenetrable layers. We further neglect interactions driven by deformations induced by the z-component(s) of the smectic forces as higher order in our perturbation analysis ( ∼ 4 ) and thus small. The smectic mediated-force exerted on one swimmer by another at a separation can then be computed in Cartesian coordinates:

Fig. 1
Microswimmers or other inclusions locally distort the smectic layers. This distortion has an associated potential energy, the gradient of which which results in a force, Sm , that acts to move the inclusions towards their energetic minima at in-plane separation r �� = √ 8 �z� ( r || = 4 for particles separated by two layers, as in this sketch), indicated by the dashed green radius (color figure online) 1 3 In a system of N swimmers, the force on a particular swimmer at is the sum of the forces exerted on it by all the other swimmers in the system, located at i : where ∇ acts on the variable.
In what follows we neglect hydrodynamic interactions between the swimmers, noting that these may be strongly screened by the lamellar layers.

Swimmer reorientation
An embedded micro-swimmer, in the absence of any smectic-mediated force, is assumed to travel with self-propulsion velocity = v sp and speed v sp , with = (cos , sin , 0) the swimming direction, parallel to the plane of the smectic layers, and an angle measured relative to the x-axis. Rotations about the p-axis likely occur but are neglected as the particle is assumed to be axisymmetric. We neglect rotations about an axis that is perpendicular to p and parallel to the smectic layers on the grounds that they are suppressed by steric interactions with the neighbouring layers, although these could later be incorporated if there is evidence that they occur. The only rotation we consider here is the yaw about the axis perpendicular to the smectic layers.
The swimmers will have some centre of frictional or hydrodynamic stress, located at position H [17]. This is where the hydrodynamic force arising through their motion through the fluid may be considered to be acting and depends on the particle shape. Any force acting on the swimmer through a point other than H will generate a torque and hence a rotation. The smectic-mediated interaction is, in general, such a force. The torque that it generates Sm = ∧ Sm where is the (signed) offset between the centres of smectic and hydrodynamic forces, see Fig. 3. This torque acts to increase the angle between and Sm if is negative, and to decrease it if is positive. This is to say, that it acts to align the swimmer with the direction of the smectic force if the swimmer is a pusher, and to anti-align it if a puller, see Fig. 3.
The viscous torque due to rotation that balances this smectic torque also depends on the shape of the swimmer. For a spherical swimmer it is given by Faxen's law = 8 a 3 (̇ ∧ ) . However, for a swimmer with more general axisymmetric shape, it takes the form V = c 1 (̇ ∧ ) , where is the dynamic viscosity, and c 1 is a parameter (with dimensions of volume) that depends on the shape. Following Newton we balance the torques V + Sm = 0 to find − sin cos we can then write

Centre of mass motion
The overdamped equation of motion for the particle position is assumed to be  Fig. 2 The two-body interaction potential G( ) between particles separated by z∕ ∈ [1,5] layers share the common features of a repulsive maximum near r || = | || | = 0 for particles that lie above one-another and attractive wells for particles that are offset. The lateral interactions are therefore repulsive at short range and attractive at long range, decaying exponentially in r || much faster than an algebraic decay in z. The potential landscape can be much more complex when many particles are present, not just two Fig. 3 The offset between the smectic force centre Sm and the hydrodynamic centre H , can result in a change in the agent's swimming direction, . a If the distance = + > 0 , the resultant torque will reduce the angle between and Sm , aligning the swimming direction with the smectic force (red). This is probably the natural sign for most pushers. We adopt this convention in the remainder of this paper. b If = − < 0 , the natural sign for most pullers, the swimmer will instead anti-align (blue) (color figure online) where c 2 is a parameter (with dimensions of length) that depends on the shape of the swimmer, see Fig. 4. Equations (4), (6) and (7), define our dynamical model.

Simulation details
Our simulations (in Python) are performed in a periodic box comprising ten smectic layers (and hence interlayer gaps) each 50 across, where is taken to be the distance between layers ( 50 × 50 × 10 in dimensionless units). Unless noted otherwise we study N = 100 agents, 10 agents per inter-layer gap, with initial x and y positions and orientations sampled uniformly from all allowed values. We use a dimensionless timestep 0.001. Our swimmers are modelled as point particles with no direct (steric) interaction with other agents in the same layer. Swimmers in the same layer can approach arbitrarily closely. For strong smectic forces (large S) this occurs frequently and can result in some agents persistently residing in the same local potential well. Whilst such arbitrarily close inter-particle separations could be avoided by the incorporation of some short ranged, same-layer repulsion we chose not to do this because local clustering at some small but non-zero separation will anyway persist if the smectic field favours it. We initially ran simulations for S, ∈ {0.001, 0.01, 0.1, 1, 10, 100, 1000} . Position and orientation data was collected for all agents in the simulation over 10000 time-steps following a burn-in of 150000 time-steps (sufficient time for a swimmer to traverse the periodic box several times). While this corresponded to the limits of our available computational resources it nonetheless remains difficult to tell if a some sort of steady state is reached in all of these simulations, or even if one exists at all. To examine order in the system we ran simulations for subsets of the parameter values. Running the simulation for 500000 timesteps allowed us to see how the order changed over longer timescales, and informed our choice to collect time averaged data for 200000 timesteps, following a burn-in of 300000 timesteps, when exploring the impact of noise on the system.

Order parameter
To characterise the amount of order in the system, we monitor two order parameters in our simulations: (i) the average smectic deformation energy per particle in units of S, written H, and (ii) the co-alignment of swimmers characterised by the order parameter where i is the swimming direction of swimmer i. A disordered system, with agents swimming in random directions, will have ≈ 0 , with = 1 when all agents are perfectly aligned. Note that under this (conventional) definition particles that move persistently in random directions (the "ballistic state" discussed below) and particles that rapidly

Results and discussion
The phase diagram shown in Fig. 5 summarises the different behaviours observed in our simulations while Fig. 6 shows snapshots of the underlying dynamics. The emergence of collective motion in these simulations is very noteworthy. Our model lacks any explicit co-alignment of agents, like that incorporated in the Vicsek model [18]. Rather, the collective motion seen here arises due to forces associated with the smectic-distorting presence of swimmers in other layers. These generate torques that give rise to alignment. It is this effect that generates collective motion. This alignment process requires the smectic to be sufficiently stiff to generate adequate forces (S large enough) and also that the particles re-align fast enough ( not too large) and can be understood as follows: A group of agents that each reside at or near a local minimum in the smectic distortion potential can remain near these minima if they are perfectly aligned, as the swarm undergoes only a collective translation. If any agent becomes mis-aligned it will drift out of its local minimum and the smectic will act to restore it, both directly, through the action of an advective Sm but also by realignment, due to this advection. Figure 7 shows how the orientational order (t) and distortion energy per particle H(t) change with time.
We see from these results that the choice of parameters S, has an impact on the time it takes the system to relax and/or the state that it ultimately targets. Rapid reorientation times lead to more stable (negative) values of total energy H. These systems remain frustrated as a result of the strong short-range smectic interaction making it very difficult for deterministically driven swimmers to escape potential wells and find new minima; once they find a relaxed state, they rarely leave. This motivates us to carry out a preliminary study of the role of added noise, omitted elsewhere in this study, see Fig. 8. As in classical models, with explicit coalignment, we observe an order-to-disorder transition on increasing noise [18,19].

Summary
We have demonstrated that collective motion can emerge in systems of swimmers embedded in a smectic liquid crystal. This arises as a consequence of indirect, smectic-mediated interactions. Orientationaly ordered states can emerge and these undergo an order-to-disorder transition on the addition of orientational noise. We hope that this study will motivate experimental interest in corresponding systems of active smectic liquid crystals. Further studies might focus on introducing hydrodynamic or direct particle-particle interactions. "Ballistic motion" refers to particle dynamics substantially similar to non-interacting particles. This occurs when the smectic field is sufficiently weak S ≲ 1 ; no order develops in the system. "Collective motion" refers to coaligned groups with quasi-regular spacing; high global order eventually develops. "Clustering" refers to closely grouped particles that remain non-aligned, due to slow reorientation. These groups then drift slowly, if at all. "B/C" refers to situations where the agents are on the transition from Ballistic to Collective motion. We speculate that even longer simulations may result in the line separating Ballistic and Collective motion moving to the left (smaller S), and the line between Clustering and Collective motion moving up (larger ). We speculate that collective motion might ultimately be expected to have greater stability in the long-term in the absence of added noise  The orientational order (t) (left panels) and distortion energy per particle H(t) (right panels, inverted scale) change during 5 × 10 5 microscopic time-steps (colours show four different realisations from different random initial conditions). Each row corresponds to a different combination of the smectic stiffness S and reorientation time .
Note that for < 1 , a high degree of orientational order and energetic order are found, whereas this occurs only very slowly for = 10 , and no change in order is seen at all for = 100 over the duration of the simulations. Increases in the energetic order are accompanied by temporary drops in orientational order (color figure online) Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. The same microswimmer-smectic system undergoes an orderto-disorder transition on increasing the strength of rotational noise . The order depends on the noise strength. In these simulations a noise term was added to the right hand side of Eq. 11. This corresponds to the addition of a stochastic angle to (t) every timestep that is itself the product of a random scalar chosen uniformly on [− , ] and the timestep magnitude t = 10 −3 . Shown are simulations with ( , S) = (0.1, 1) orange, (1, 10) red, (1, 1) blue (1, 10) green; an order-disorder transition is clear in all but the last of these. Each simulation was run for 3 × 10 5 time-steps as a burn-in period, with the order parameters then time-averaged over a subsequent 2 × 10 5 time-steps. Error bars show the standard error in the mean (color figure online)