Deformable active nematic particles and emerging edge currents in circular confinements

Abstract We consider a microscopic field theoretical approach for interacting active nematic particles. With only steric interactions the self-propulsion strength in such systems can lead to different collective behaviour, e.g. synchronized self-spinning and collective translation. The different behaviour results from the delicate interplay between internal nematic structure, particle shape deformation and particle–particle interaction. For intermediate active strength an asymmetric particle shape emerges and leads to chirality and self-spinning crystals. For larger active strength the shape is symmetric and translational collective motion emerges. Within circular confinements, depending on the packing fraction, the self-spinning regime either stabilizes positional and orientational order or can lead to edge currents and global rotation which destroys the synchronized self-spinning crystalline structure. Graphical abstract Supplementary Information The online version contains supplementary material available at 10.1140/epje/s10189-022-00162-w.


Introduction
Active matter systems take energy from their environment and drive themselves out of equilibrium. This can lead to novel collective phenomena and provides hope to uncover the physics of living systems and to find new strategies for designing smart devices and materials. We refer to Ramaswamy [1], Marchetti et al. [2], Cates and Tailleur [3], Menzel [4], Bechinger et al. [5] and Gompper et al. [6] for various reviews. An important example of active matter is constituted by natural and artificial objects capable of self-propulsion. A fundamental challenge is to understand how such objects interact and lead to collective phenomena. Most of the microscopic modelling approaches in this field consider active particles which have a fixed symmetric shape, and movement is defined along a symmetry axis. This leads to motion along a straight line just perturbed by random, e.g. Brownian fluctuations. Both assumptions, on shape and symmetry, are restrictive, as shape deformations as well as deviations from symmetry destabilize any straight motion and make it chiral, which would result in circular motion. As most systems are imperfect this should be the general case. While attempts exist to generalize active particle models in this direction, see, e.g. [7][8][9][10] for imposed alignment mechanisms, Denk et al. [11] and Bär et al. [12] for anisotropic particle shapes and Ohta and Ohkuma [13] and Menzel and Ohta [14] for shape deformations, multiphase-field models, e.g. a e-mail: veit.krause@tu-dresden.de b e-mail: axel.voigt@tu-dresden.de (corresponding author) [15][16][17][18], where each object is modelled by a phase field variable, which implicitly described the shape of the object, and thus naturally allow for shape deformability and also provide the possibility to incorporate asymmetry to enforce chirality. It has already been demonstrated that collisions of deformable objects can lead to alignment [14,[19][20][21]. As a result, these multiphase-field models do not require any explicit alignment interactions. A drawback of such models is the huge computational effort for large numbers of interacting objects. We here consider an intermediate modelling approach. The approach considers a particle density and combines it with internal nematic structure. We are only interested in relatively dense systems and study the influence of activity in unconfined and confined domains.
The paper is organized as follows: In Sect. 2 we postulate a minimal model which is capable of shape deformations and broken symmetry with respect to the direction of motion. The model is termed nematic active phase field crystal model. Besides the motivation, the evolution equations are explained and the numerical approaches for solving and postprocessing are sketched. Section 3 analyses the model for a single object and identifies three different regimes: resting, circular or spinning motion and translation. Section 4 considers the emerging collective behaviour in unconfined and confined geometries, and Sect. 5 discusses these results and relates the observed phenomena to that of other theoretical and experimental investigations.

Motivation
Microscopic field theoretical approaches for active system can be considered as a compromise between the full details of multiphase-field models and active particle models. They first have been introduced in [22] for active crystals and consider a local particle density variation field ψ, a local polar particle orientation field P and a self-propulsion strength v 0 . The model combines a phase field crystal model for freezing [23,24] with a Toner-Tu model for self-propelled particles [25]. More recently this approach was also considered on surfaces [26] and has been extended by an active torque, and the interplay of self-propulsion and self-spinning of crystallites was investigated in [27]. A different path was followed in [28,29] where the underlying phase field crystal energy was modified to consider independent active particles [30][31][32][33]. The approach allows to simulate a transition from a resting particle to a moving state by increasing the self-propulsion strength. Within this transition, the particle deforms and elongates perpendicular to the direction of motion. Other phenomena, considered for more particles, are cell-cell collisions, oscillatory motion in confined geometries, collective migration and cluster formation in homogeneous systems [28], as well as the rich dynamics of heterogeneous systems of active and passive particles, ranging from highly dilute suspensions of passive particles in an active bath to interacting active particles in a dense background of passive particles [29]. A common characteristic of these models is the presence of an underlying interaction potential for the density variation field ψ but no enforcement of aligning the polar particle orientation field P. Alignment results solely from inelastic particle deformations through their interaction [19].
For a single particle the interplay between a splay (or bent) instability of the polar particle orientation field P, the particle shape and the strength of the selfpropulsion v 0 has been discussed [28] and corresponds to the same mechanism as in phase-field models for active droplets, see, e.g. [34,35]. As a result of the vertical anchoring of the polar particle orientation field P at the particle boundary, one +1 defect forms within the particle. With no further interaction, due to the isotropic properties of the particle, the resulting shape of the particle is symmetric with respect to the direction of motion, see Fig. 1 (left). To incorporate chirality thus requires an additional active force , as in [27], or a different particle orientation field. Adapting approaches of active nematic droplets [36,37], we propose a microscopic field-theoretical approach, which couples a local particle density variation field ψ, a local nematic particle Q-tensor field Q and a self-propulsion strength v 0 . Similar mechanisms, as described above, also follow for this model, but now the nematic properties lead to the presence of two +1/2 defects, which allows to break the symmetry and induces chirality, see Fig. 1 (middle). This property allows to consider only one active parameter, the self-propulsion strength v 0 , to tune the rich dynamics of the model.

Evolution equations
The proposed minimal model reads with particle density variation field ψ, nematic particle Q-tensor field Q and self-propulsion strength v 0 . The first equation considers conserved dynamics for the free energy the Swift-Hohenberg energy [23,24,38], with parameter r related to an undercooling, and an additional penalization term, with parameter H > 0. The penalization enforces the density variations to remain positive. This modifies the particle interaction and allows to phenomenologically describe independent particles [30][31][32][33]. A detailed derivation of F P F C and its relation to classical density functional theory can be found in [39][40][41]. The variational derivative reads δFvP F C δψ = (r + 1)ψ + 2∇ 2 ψ + (∇ 2 ) 2 ψ + ψ 3 + 3H(ψ|ψ| − ψ 2 ). The parameter M 0 sets a mobility and is responsible for the deformability of the density peaks. The active contribution is considered in analogy to the polar model [28], with Q∇ψ playing the role of the polar particle orientation field P. The second equation considers unconserved dynamics of a Landau-de Gennes energy in its one-constant approximation with elastic constant L and entropic parameters b = 0 (as only two-dimensional systems are considered) and a = −c. The active component is constructed to ensure the Q-tensor properties and the last term restricts, in analogy to Alaimo et al. [28], the nematic particle Qtensor field Q to be different from zero only within the particles, with β > 0. The model is considered in a nondimensional setting.

Numerical approach
The coupled equations are reformulated as a set of second-order equations and solved using an operator splitting approach for ψ and Q in a semi-implicit manner. Discretization in space is done by finite elements Fig. 1 (Right) Particle density variation field ψ for a single density peak (colour coding) with the 0.01 level set indicating the shape of the particle. With the 0.01 level set we consider 95% of the mass of the density peak. The internal polar (left) and nematic (middle) structure is visualized by the director field. The figures further show the direction and strength of motion (arrows) and the location of topological defects (black points). For the polar model (left) one +1 defect is located on the symmetry axis, and for the nematic model (middle) two +1/2 defects break the symmetry and make the motion chiral [42,43], and adaptive refinement is considered to ensure a fine discretization within the particles. The approach is implemented in AMDiS [44,45]. We consider a square domain Ω = [−6d, 6d] 2 , with periodic boundary conditions, where d = 4π √ 3 the lattice distance of the phase field crystal model. A circular confinement is enforced using an interaction potential to be added to F vP F C , which reads Bψ 2 ϕ B dr with B > 0 and ϕ B a tanh-approximation of 1 Ω\Ωc , with The model parameters are fixed as r = −0.9, M 0 = 20, L = 0.2, c = 0.1, H = 10 5 , β = 10 and B = 10 5 . The self-propulsion strength v 0 will be varied and specified below. Numerical parameters concerning grid resolution, time step and tanh-approximation are chosen to guarantee mesh-independency and stable behaviour.
As initial condition we specify with prefactor A such that ψ 0 dr = Nd 2 /|Ω| (−48 − 56r)/133 and particle initial positions r i for i = 1, . . . , N with N the number of particles. As initial Q-tensor field we consider a symmetric field with one +1 defect in the centre of each particle and vertical anchoring at the particle boundary. The initial Q-tensor field is perturbed by white noise to break the symmetry. The +1 defects are unstable and immediately split into two +1/2 defects. The way these defects rearrange sets the shape of the particle and its direction of movement.
For postprocessing purposes the centre of the i-th particle at time t n is computed as r n i = Bi rψ n dr/ Bi ψ n dr, with B i a small circle around the maximum of the i-th density peak. The radius of B i is related to d. The i-th particle velocity follows as v n i = (r n i − r n−1 i )/(t n − t n−1 ), and the mean particle velocity magnitude is the average over all v n i , computed [20,28] we define for every time t n the translational order parameter φ n T and the rotational order is the unit i-th particle velocity vector andr n i = r n i / r n i the unit i-th particle position vector at time t n . In case of collective translation or collective rotation, we get φ T,i ≈ 1 or |φ R,i | ≈ 1, respectively. However, also collective orientation in synchronously spinning particles leads to φ T,i ≈ 1. To distinguish translational and orientational order φ O measures synchronously changing orientation. The frequency of the oscillation in φ O (t) determines the collective angular spinning velocity.

Single particle
We first consider the situation of one particle. It is placed in the centre of the domain, and we consider the effect of v 0 . Figure 2 shows the particle velocity, the eccentricity and the asymmetry of the defect arrangement as a function of v 0 . The eccentricity is defined as are the minimal and maximal distances between the centre of mass and the 0.15 level set of ψ for particle i at time t n , respectively. The 0.01 level set is considered as the particle boundary. The asymmetry of the defect arrangement is computed as the deviation from the centre of mass with respect to length and angle, as a n where d n i,1 and d n i,2 are the positions of the two +1/2 defects for particle i at time t n . Various approaches exist to determine defects in nematic liquid crystals, see [46] for a comparison of various methods. We here consider them as degenerate points of Q for which Q 11 = Q 12 = 0. This allows an easy detection of the position of a defect. For a nematic liquid crystal in 2D two types of topological defects predominate +1/2 and −1/2. Considering the sign of δ = ∂Q11 ∂x ∂Q12 ∂y − ∂Q11 ∂y ∂Q12 ∂x allows to distinguish between them. Due to the setting within a particle and the specified vertical anchoring only +1/2 defects occur in the considered parameter regime. Figure 2(left) shows three regimes: (a) resting, characterized by a zero velocity, the cell shape deforms with increasing v 0 and the defect positions are symmetric, (b) circular or spinning, the velocity fluctuates, which has an effect on the eccentricity and the asymmetry of the defect positions, and (c) translation, with increasing velocity, constant shape and symmetric arrangement of defects. The nonequilibrium steady state of the circular or spinning regime is shown in Fig. 2(right) for different v 0 . The oscillations underpin the correlation between velocity, eccentricity and defect asymmetry. While they are strongest for v 0 = 1.5, they decrease for v 0 = 1.75 and are almost gone for v 0 = 2.0, in accordance with the error bars in Fig. 2(left).
To further highlight the connection between particle velocity, eccentricity and asymmetry of the defect positions Fig. 3(left) shows the particle shape together with the principle eigenvector of the largest eigenvalue of Q (director field) and the defect positions for various v 0 . As the defects can also be located at the 0.01 level set, the nematic liquid crystal, which is forced to decay to zero in regions with ψ ≤ 0, is also shown in the vicinity of the particle. The defects deform the director field and the deformed director field is responsible of the symmetry breaking. For the circular or spinning regime the shape deformation is asymmetric with respect to the direction of movement and the defect asymmetry increases with v 0 . The direction (up or down) depends on the splitting of the +1 defect into two +1/2 defects and the resulting shape deformation. For the translation regime the shape is symmetric with respect to the direction of motion and also the defect arrangement is almost symmetric. With increasing v 0 the defects are located closer to the symmetry axis and the velocity of movement, which only slightly deviates from the symmetry axis, increases. Figure 3(right) shows a typical circular path together with the corresponding director field and the velocity in the shown time instances. Due to the small radius of the circulation, which is almost independent on the strength of activity v 0 , we denote this motion as spinning in the following.
Having spinning and translation as possible motility modes within one model, determined by the strength of activity, offers the possibility to switch between both modes. If these switches occur randomly in time and the spinning mode randomizes the translation direction, the behaviour is reminiscent to run-and-tumble particles and, however, here controlled by regulating the activity.

Collective behaviour
The behaviour in the resting and translation regimes essentially coincides with that of the polar active phase field crystal model [28]. This also remains true for the emerging collective behaviour in unconfined and confined geometries, see Appendix A. We thus only concentrate on the spinning regime in more detail. First, we characterize the behaviour of interacting spinning particles in unconfined and confined geometries for an intermediate packing fraction of 0.57. To compute the packing fraction we consider the 0.01 level set of ψ to determine the area of the particles as A N = I {ψ>0.15} dr. The area of one particle A = A N /N ≈ 0.9d 2 with d = 4π/ √ 3 the lattice distance in the phase field crystal model. This essentially motivates to consider the 0.01 level set. The packing fraction results as A N /|Ω| or A N /|Ω c |.

Synchronization in unconfined and confined geometries
We first consider 120 particles in the square domain Ω with periodic boundary conditions. The self-spinning particles form crystalline structures with local triangular order, with dislocations and regions with no particles, which dynamically rearrange. The particles are self-spinning, and due to local interactions some particles also move to positions further away than the spinning radius. This is consistent for all considered selfpropulsion strength v 0 . However, only for v 0 = 2.0 the translational order parameter φ T ≈ 1, which indicates translational order or in the current context synchronized spinning. This is confirmed by the angular order parameter φ O , which oscillates with fixed periodicity, see Fig. 4a, b. The behaviour in the circular confinement Ω c is similar, see Fig. 4c, d. Also in this setting the translational order parameter φ T ≈ 1 for v 0 = 2.0 and the angular order parameter φ O oscillate with fixed periodicity. This nonequilibrium steady state is reached much faster than in the unconfined geometry. One could conclude that in this setting confinement helps to synchronize the particles. Deviations from synchronized spinning in the reached nonequilibrium steady state are only found at the edge. This corresponds with regions with crystalline defects or no particles. In the centre a crystal with perfect triangular lattice and synchronously spinning particles emerges.
In contrast to the translational regime considered in Appendix A with translational and rotational motion as nonequilibrium steady states, here both settings behave similar. In unconfined and confined geometries the initially independently spinning particles undergo a transition to a nonequilibrium steady state of positional and orientational order, a synchronized spinning crystal. The simulations only confirm this for v 0 = 2.0. If this state is also reached at later times for the other self-propulsion strength v 0 remains open.

Varying packing fraction and emerging edge currents
All previous simulations consider a packing fraction of 0.57. We now vary this in the circular confinement and observe different behaviour, see Fig. 5. For a smaller packing fraction of 0.48, at least within the considered time neither synchronized spinning nor crystal formation can be observed. Instead only small crystalline patches form and dynamically rearrange. Due to the available space particle interactions lead to local positional rearrangements. Similar to the situation in unconfined geometries particles can move to positions further away than the spinning radius. Figure 5b shows the coarse-grained trajectories of the particles (without the spinning component), and Fig. 5c  of 0.48 underpin the described behaviour. The bond number gives an indication of crystalline order and is computed for particle j as b n 6j = ( k∈Nj e 6iθ n jk )/N j , with N j the nearest neighbours of particle j within a specified radius related to d and θ n jk the angle between r n k − r n j and the x-axis. The considered averaged bond numberb 6j accounts for the average over various times t n .b 6j = 0 considers the situation of an isolated particle andb 6j = 1 that of a perfect triangular lattice, a particle with six neighbours. The nearest neighbours The situation changes for increasing packing fraction. For 0.66 and 0.7 the dominating situation of a crystal with triangular lattice and synchronously spinning particles is destroyed. The edge currents increase and propagate towards the centre, see Fig. 5b. While for 0.66 a triangular lattice still exists at least over some time span before it gets rearranged, the averaged bond order for 0.7 has even less indication of such stable crystalline order, see Fig. 5c. The coarse-grained particle trajectories, see Fig. 5b, show a transition towards a global vortex. The snapshots in Fig. 5a further indicate that the particles no longer spin synchronously. The Supplementary Movies further confirm this behaviour. The coarse-grained particle movements (without the spinning component) are further analysed in Fig. 6, which confirms the above discussion. The kymographs show the orthoradial and radial components of the coarse-grained velocity averaged over all particles which are located at a distance R from the centre of the domain Ω c . While there is almost no movement in orthoradial direction, the slight edge currents for packing fraction 0.57 in the radial component and their increased strength and extension towards the centre for packing fractions 0.66 and 0.7 are clearly visible. The direction of the emerging vertex rotation depends on the shape of the particles at the boundary as a result of their spontaneous symmetry breaking. The majority decides on the emerging direction at the edge, which persists towards the centre.

Discussion
The proposed minimal model allows to explore different dynamical regimes by varying one activity parameter only. The direct coupling of the self-propulsion strength v 0 with the internal nematic structure and the deformability of the particle lead to slightly deformed resting states if v 0 is below some threshold. It induces within a certain parameter range chirality, which leads to circular or spinning motion. If it is above some threshold, a symmetric shape and translational motion emerge. All regimes have been investigated in unconfined and confined geometries. While the translational regime is more or less identical with the behaviour of the active polar phase field crystal model [28] and the observed rotational behaviour in circular confinements reminiscent of various experiments, e.g. on highly concentrated bacterial suspension, which self-organize into a single stable vortex [47], collective behaviour of self-circulating or self-spinning particles is much less explored. Selfspinning particles are computationally considered in circular confinements [48]. While fundamental issues differ, e.g. our particles are deformable, our spinning radius is significantly larger and our spinning velocity significantly lower, the emerging behaviour is similar. The competition between circular confinement, selfpropulsion and steric interactions can lead to the emergence of edge flows and rotations. Within a wider perspective, similar edge flows and rotations have also been observed in chiral fluids [49]. It is shown that in systems of synchronously spinning particles both parity (or mirror) symmetry and time-reversal symmetry are broken. Hydrodynamic theories with additional terms to account for these broken symmetries, e.g. rotational viscosity, tend to force the fluid as a whole to rotate with the angular velocity of the spinning particles. However, the motion of the fluid in the bulk is suppressed by friction. As a result, the fluid moves mostly at the boundary and the penetration depth of the vorticity of the fluid from the boundary into the bulk is controlled by the shear viscosity. These results are similar to the edge currents in our simulations and their propagation towards the centre with increasing packing fraction. These similarities with other systems range from collective rotation of chiral molecules of a liquid crystal [50], to sperm cells [51,52], colloidal and millimetre scale magnetes [53,54] and rotating robots [55]. An interesting biological example is provided by Chlamydomonas reinhardtii (C. reinhardtii). This micron-sized unicellular algae is able to self-propel to perform translational motion, but also has the ability to self-rotate [56]. Rotation is used to sense the direction of light to optimize efficiency of phototaxis [57].
The proposed microscopic field theoretical model can be extended towards various directions, e.g. hydrodynamic interactions. This is already considered within the phase field crystal model for passive systems, e.g. [58][59][60][61]. Other possibilities consider multicomponent systems [29,39].

A Collective behaviour in translational/rotational regime
If the self-propulsion strength is above some threshold, a symmetric arrangement of defects and a symmetric shape of the particles is enforced. As a result, the particles behave as in the polar active phase field crystal model [28]. This behaviour leads to collective translation in unconfined geometries and collective rotation in circular confinements, see Fig. 7. It shows simulations for different values of v0. Collective behaviour is only reached in the considered simulation time for the largest values, ϕT ≈ 1 for v0 = 4.0 and |ϕR| ≈ 1 for v0 = 3, 5 and 4.0. This behaviour is in qualitative agreement with results in [28] and corresponding large-scale simulations of multiphase-field models [62].