Dancing disclinations in confined active nematics

TL;DR

This paper demonstrates how confining active nematic fluids in channels induces organized disclination patterns and vortex lattices, enabling control over chaotic flows and paving the way for active microfluidic devices.

Contribution

It reveals a method to achieve ordered disclination arrangements and flow structures in confined active nematics, advancing the design of active topological microfluidic systems.

Findings

Regular motion of disclinations observed in channels

Formation of dynamic vortex lattice

Disclination pairs exchange partners, creating an ordered state

Abstract

The spontaneous emergence of collective flows is a generic property of active fluids and often leads to chaotic flow patterns characterised by swirls, jets, and topological disclinations in their orientation field. However, the ability to achieve structured flows and ordered disclinations is of particular importance in the design and control of active systems. By confining an active nematic fluid within a channel, we find a regular motion of disclinations, in conjunction with a well defined and dynamic vortex lattice. As pairs of moving disclinations travel through the channel, they continually exchange partners producing a dynamic ordered state, reminiscent of Ceilidh dancing. We anticipate that this biomimetic ability to self-assemble organised topological disclinations and dynamically structured flow fields in engineered geometries will pave the road towards establishing new active topological microfluidic devices.

Equations

$$\left({\partial }_t+\mathbf{u}\cdot \nabla \right)\mathbf{Q}-\mathbf{S}$$

$$\mathbf{S}$$

$$\mathbf{H}$$

$$f$$

$$\nabla \cdot \mathbf{u}$$

$$\left({\partial }_t+\mathbf{u}\cdot \nabla \right)\mathbf{u}$$

$${\mathbf{\Pi }}^{\text{elastic}}$$

$$-\nabla \mathbf{Q}:\frac{\delta f}{\delta \nabla \mathbf{Q}}+\mathbf{Q}\cdot \mathbf{H}-\mathbf{H}\cdot \mathbf{Q},$$

$${\mathbf{\Pi }}^{\text{act}}$$

Full text

I Introduction

Bacterial suspensions Wensink et al. (2012), cellular monolayers Vedula et al. (2012) and sub-cellular filament/motor-protein mixtures Sumino et al. (2012); Sanchez et al. (2012) are continuously driven far from equilibrium through intrinsic energy injection by their biological constituent elements. This microscopic energy input can result in the spontaneous emergence of large-scale collective behaviour, including flocking Toner and Tu (1998); Cavagna et al. (2015); Attanasi et al. (2015), unidirectional flows Voituriez et al. (2005), meso-scale turbulence Dombrowski et al. (2004); Wensink et al. (2012); Creppy et al. (2015), topological defect pair production Giomi et al. (2013); Thampi et al. (2014a); Giomi et al. (2014), and phase separation Putzig and Baskaran (2014). Moreover there is a rapidly growing list of biological systems, including anterior-posterior establishment by active cytoplasmic streaming in oocytes Khuc Trong et al. (2015), pupal wing morphogenesis Etournay et al. (2015), wound healing Poujade et al. (2007); Lee and Wolgemuth (2011), and cancer invasion Wolf et al. (2007); Janet et al. (2012), that has been identified as leveraging collective flows for active self-organisation. Recent works have shown that nematic disclinations in growing bacterial colonies can determine morphological changes Doostmohammadi et al. (2016a) and play a crucial role in the dynamic structure of cell groups such as in ameoboid cells Kemkemer et al. (2000); Duclos et al. (2014) and spindle-shaped cells Duclos et al. (2016), and the emergence of 3D mounds in layers of stem cells Kawaguchi et al. (2016).

These examples suggest the possibility of harnessing active driving by manufacturing active microfluidic devices. In a few engineered instances, spontaneous collective motion has been controlled through geometrical constraints. Ratchet systems have been constructed to sort cells against entropy Galajda et al. (2007); Kaehr and Shear (2009); Mahmud et al. (2009), while individual asymmetric gears Hiratsuka et al. (2006); Angelani et al. (2009); Di Leonardo et al. (2010); Sokolov et al. (2010) and arrays of smooth rotors Thampi et al. (2016) have been designed to have persistent rotation when submersed in bacterial baths. Long-range static orientational ordering of nematic disclinations has recently been probed in experiments on filament/motor protein mixtures DeCamp et al. (2015) and theoretical works have predicted positional disclination ordering in crystal-like configurations due to hydrodynamic screening in over-damped active nematics Doostmohammadi et al. (2016b); Putzig et al. (2016); Oza and Dunkel (2016); Srivastava et al. (2016). Nevertheless, the ability to control disclination dynamics in active matter is currently lacking and is a necessary prerequisite to the control of active topological microfluidics.

Many dense active systems are well modelled as active nematics, and confining such an active nematic is known to result in spontaneous symmetry breaking leading to self-sustained, unidirectional laminar flows at moderate activities Voituriez et al. (2005); Woodhouse and Goldstein (2012); Ravnik and Yeomans (2013); Edwards and Yeomans (2009). It is further known that such unidirectional active flows transition to oscillatory flows at sufficiently high activity Giomi et al. (2012); Ramaswamy and Jülicher (2016) and to meso-scale turbulence at higher activities yet Wensink et al. (2012); Dunkel et al. (2013); Thampi et al. (2014a, b); Giomi (2015). Linear stability analysis has been performed on the formation of oscillations in the laminar flow with moderate activity Giomi et al. (2012) but increased activity causes non-linear effects to quickly dominate. It has recently been shown Doostmohammadi et al. (2016c) that moderate activities in confining (quasi-1D) geometries can produce a system-spanning lattice of flow vortices intermediate to oscillating flows and meso-scale turbulence. The vortex lattice has been observed in simulations of sheared active fluids Słomka and Dunkel (2016), as well as experimentally in quasi-1D suspensions of bacteria Wioland et al. (2016) and associated simulations Wioland et al. (2016); Theillard et al. (2017). Similarly, spatially ordered textures and disclination lattices have previously been noted in confined and flowing passive liquid crystals Grecov and Rey (2003); Gupta et al. (2005); Rey (2010).

Here, we confine an active nematic within a two dimensional micro-channel. The results of our simulations show that the intermediate activity state is characterized by dancing disclination dynamics, where the stationary vortex lattice is accompanied by a dynamically ordered state of disclinations (Fig. 1(a); Movie 1). The disclination dynamics accompanying the vortex lattice have not hitherto been reported. These dancing disclinations are positively charged, long lived, and continually navigate through the channel. As they move past each other on the vortex lattice, the positive disclinations form short-lived pairs, leading us to dub this state topological Ceilidh dynamics, in rough analogy to the traditional Gaelic dance in which participants form two inward-facing parallel lines and pairs of dancers continually exchange partners as they circulate between the lines. Our results for quasi-1D active nematics show that, in this simple geometry, ordered dancing-disclination dynamics emerge as a system-spanning intermediate regime between unidirectional flow and meso-scale turbulence.

We next treat the Ceilidh dynamics state as a spatially ordered configuration between the unidirectional flow state and meso-scale turbulence. As a system-spanning ordered state, the Ceilidh dynamics pattern can itself possess irregularities, which we term lattice defects. A subset of these lattice defects are found to drive a well-defined net flux that increases incrementally with the quantized number of drift-lattice defects. The dynamically ordered topological and steady-flow structures of Ceilidh dynamics within a channel represent an ideal system for studying the emergence of ordered dynamical structures in active matter.

II Methods

Many living fluids can be modelled as continuous active nematic liquid crystals because they generally consist of a dense suspension of many (continuum) shape-anisotropic (nematic) self-propelled (active) particles moving through a fluid medium (force-free) Marchetti et al. (2013). We utilize numerical simulations of active nematohydrodynamics to solve for the density, velocity (with its associated vorticity $\boldsymbol{\omega}\left(\boldsymbol{r},t\right)$), and orientation tensor $\boldsymbol{Q}(\boldsymbol{r},t)$ (with associated director $\boldsymbol{n}\left(\boldsymbol{r},t\right)$ and scalar order $q\left(\boldsymbol{r},t\right)$) fields Marenduzzo et al. (2007). We assume that the nematic fluid has a single nematic elastic constant $K$ and that the activity coefficient $\zeta>0$, representing extensile active systems such as pusher-type bacteria Lauga and Powers (2009). The 2D active nematic is confined between two parallel no-slip channel walls separated by a distance $h$, producing a quasi-1D system. Strong anchoring at the boundaries sets the director perpendicular to the confining walls.

II.1 Governing equations

Active nematohydrodynamics have been extensively applied to intrinsically out-of-equilibrium biological systems composed of rod-like constituents, including bacterial suspensions Volfson et al. (2008), filament/motor protein mixtures Sanchez et al. (2012); Giomi et al. (2013); Thampi et al. (2013), and cellular monolayers Bittig et al. (2008); Duclos et al. (2014); Doostmohammadi et al. (2015). The continuum fields that must be solved are the total density $\rho$, the velocity $\boldsymbol{u}$, and the orientation $\boldsymbol{Q}$ fields. The orientational order of the fluid is described by the nematic tensor $\boldsymbol{Q}=\frac{3q}{2}(\boldsymbol{n}\boldsymbol{n}-\boldsymbol{I}/3)$, where $q$ denotes the magnitude of the orientational order, $\boldsymbol{n}$ is the director, and $\boldsymbol{I}$ the identity tensor de Gennes and Prost (1995).

The transport equation for the nematic tensor field is Beris and Edwards (1994)

[TABLE]

where $\Gamma$ is a rotational diffusivity and $\boldsymbol{S}$ is the co-rotational advection term that accounts for the impact of the strain rate $\boldsymbol{E}=(\boldsymbol{\nabla}\boldsymbol{u}^{T}+\boldsymbol{\nabla}\boldsymbol{u})/2$ and vorticity $\boldsymbol{\Omega}=(\boldsymbol{\nabla}\boldsymbol{u}^{T}-\boldsymbol{\nabla}\boldsymbol{u})/2$ on the director field. The co-rotational advection has the form

[TABLE]

where $\boldsymbol{\mathcal{Q}}\equiv\boldsymbol{Q}+\boldsymbol{I}/3$ and the alignment parameter $\lambda$ controls the degree of coupling between the orientation field and velocity gradients and determines the objective time derivative of orientation, with $\lambda=\pm 1$ corresponding to an upper and lower convected derivative respectively. The alignment parameter determines whether the nematogens align or tumble in a shear flowThampi et al. (2014b). The relaxation of the orientational order is controlled by the free energy \mathcal{F}$$=\int fdV through the molecular field,

[TABLE]

The free energy density has two components, which are the Landau-de Gennes bulk free energy density and the elastic free energy density due to spatial inhomogeneities in the order parameter de Gennes and Prost (1995):

[TABLE]

assuming a single elastic Frank coefficient $K$.

The local density and velocity field obey the incompressible Navier-Stokes equations

[TABLE]

where $\boldsymbol{\Pi}$ is the generalized stress tensor that includes both nematic and active contributions, in addition to the viscous stress $\boldsymbol{\Pi}^{\textmd{visc}}=2\eta\boldsymbol{E}$. The stress due to elastic contributions arising from nematic ordering within the liquid crystal is

[TABLE]

which includes the isotropic pressure $P$ Beris and Edwards (1994). The active stress accounts for changes in the flow field caused by continual energy injection at the microscopic scale. Activity generates flows for nonzero gradients of $\boldsymbol{Q}$ and takes the form Aditi Simha and Ramaswamy (2002)

[TABLE]

The activity parameter $\zeta$ determines the strength of the active flows with positive and negative values denoting extensile and contractile fluids, respectively.

II.2 Numerical implementation

The equations of active nematohydrodynamics (Eq. 1-5) are solved using a hybrid lattice Boltzmann and finite difference method Marenduzzo et al. (2007); Fielding et al. (2011); Thampi et al. (2014b), with the discrete space and time steps defining the simulation units. Parameters can be mapped to physical units that correspond to microtubule-kinesin bundles Sanchez et al. (2012), dense bacterial suspensions Wensink et al. (2012), or other active nematics of interest when the material properties are well known Cates et al. (2008); Thampi et al. (2013). Simulations presented here use the parameters $\Gamma=0.34$, $K=0.04$, $\eta=2/3$, and $\rho=1$, except where explicitly stated otherwise.

The alignment parameter is taken to be $\lambda=0.3$, which is in the flow tumbling regime. The active fluid is extensile with $\zeta>0$, though Ceilidh dynamics were also observed for $\zeta<0$ and $\lambda<0$. Unless otherwise stated, $\zeta=0.04$. The Landau-de Gennes bulk free energy parameters are chosen to be $\mathcal{A}=0$, $\mathcal{B}=-0.3$ and $\mathcal{C}=−0.3$ in lattice Boltzmann units. The bulk energy density scale is therefore identified from these choices Marenduzzo et al. (2007) to be $\mathcal{E}_{G}\simeq 0.1$.

All simulations are performed in a 2D channel of height $h=25$, unless otherwise stated. The walls enforce a no-slip condition on the velocity field through bounce-back boundary conditions and strong homeotropic anchoring of the orientation field. The channel runs in the $\hat{x}$-direction with periodic boundary conditions. The channel length is $L=250$ (Fig. 2-4 and Movie 1-2), $L=300$ (Movie 3), $L=750$ (Movie 4), or $L=1000$ (Fig. 5-9). Multiple channel lengths ($L=\left\{100,250,300,500,750,1000\right\}$) are utilized in Fig. 9(a). The fluid is initialized with a zero velocity and with the director everywhere perpendicular to the walls. Initial simulations showed steady state dynamics are reached in $\sim 10000$ time steps so an initial warmup of $50000$ is utilized throughout.

III Results

For activities where a lattice of flow vortices is stable, the disclinations exist as permanent pairs of oppositely charged $\pm 1/2$ topological defects de Gennes and Prost (1995) (Fig. 1(a); Movie 1). Each $-1/2$ disclination lingers at a point in the vicinity of the channel wall, while positive $+1/2$ disclinations continually dance; navigating through the channel and following approximately sinusoidal trajectories (Fig. 1(b)). Half of the positive disclinations travel in one direction along the channel and the other half travel along the mirrored trajectory, in the ideal case (Fig. 1(c)). As they move past each other on the vortex lattice, the positive disclinations form ephemeral pairs of dancers that continually exchange partners and periodically skew the vortices as they circulate between the $-1/2$ disclinations (Fig. 1(b); Movie 1).

III.0.1 Mid-channel vortex lattice and $+1/2$ disclinations

Plus half disclinations in extensile active nematics self-propel themselves along their comet-like tail Giomi et al. (2014) and move along the borders of counter-rotating vortices (Fig. 1(b); Movie 1). Although elastic interactions with the near-wall $-1/2$ disclinations attract the $+1/2$ disclinations towards the bounding planes, the vorticity field and the motility of the $+1/2$ disclinations keep them from leaving the mid-region. A disclination cannot switch paths, since this would require overcoming both the deformation free energy barrier of crossing the region between the two paths and also of reorientation.

Increasing the active energy input in the system, increases the speed of the $+1/2$ disclinations, which is seen to increase with the activity as $v_{+1/2}\sim\zeta$ (Fig. 2). The no-slip boundary conditions viscously dissipate this energy and we see that increasing the wall separation $h$ also increases the disclination velocity linearly (Fig. 2). These dependencies are in agreement with the theoretically predicted characteristic self-motility of $+1/2$ disclinations Giomi et al. (2014) of $v_{+1/2}\sim h\zeta/\eta$ where the channel height $h$ is taken to be the characteristic system size. This agreement indicates that nematic elasticity does not dominate the dynamics. In fact, increasing the Frank coefficient from $K=0.035$ to $0.07$ produces no observable effect (Fig. 2). For $K$ outside of this region, we do not see topological Ceilidh dynamics for the simulated channel height and activity. For $K<0.035$ values we instead observe meso-scale turbulence, while for $K>0.07$ we find oscillating flows — the dynamical transitions to and from Ceilidh dynamics, and the phase diagram of dynamical steady flow regimes will be discussed in more detail in § III.1 and § III.2.

III.0.2 Near-wall region and $-1/2$ disclinations

Unlike $+1/2$ disclinations, $-1/2$ disclinations in active nematics are not self-motile Giomi et al. (2014) and their dynamics are primarily set by elastic and hydrodynamic interactions. Line tension pins the $-1/2$ disclinations to the confining walls Denniston (1996); however, the $-1/2$ disclinations can never reside exactly on the surfaces because the director field is imposed by strong homeotropic anchoring boundary conditions, which ensure that the field cannot vary substantially in the near-wall region. It is noteworthy that the dancing dynamics in the mid-channel is independent of the anchoring conditions on the walls.

Although the $-1/2$ disclinations do not transverse the channel, they move in-plane slightly about their pinning point and in a way that is directly correlated to the oscillatory trajectories of the $+1/2$ disclinations (Movie 1). The in-plane oscillations of the $-1/2$ disclinations parallel to the channel walls have the same constant period as the lateral component of the $+1/2$ disclinations (Fig. 3(a)), although the amplitude is much smaller. Interestingly, the oscillations of the non-motile $-1/2$ disclinations are generally slightly out-of-phase with the self-propelled $+1/2$ disclinations. This phase difference results from the nematic elasticity $K$ (Fig. 3(b)) but the fact that it is small again suggests that, although present, nematic elasticity is not a significant factor in these simulations of topological Ceilidh dynamics.

III.1 Dynamical transitions to and from Ceilidh dynamics

We treat the flows (unidirectional, oscillating, Ceilidh dynamics and meso-scale turbulence) as ordered states that arise from a series of dynamical transitions as a function of dimensionless activity (Movie 2). The marked difference between laminar flows (Movie 2(a)) and the intermediate Ceilidh dance state (Movie 2(c)) is the formation and structure of the vortex lattice. As the activity is increased, we first observe oscillating but unidirectional laminar flow, with temporary traveling vorticity patterns that intermittently dissipate then re-form (Movie 2(a)). For the vortex lattice to form, the size of the vortices must be commensurate with the channel height (Fig. 4(inset)). When the height of the channel is increased then the number of vortices per unit length decreases (Fig. 4). Because vortices must occur in counter-rotating pairs, we find that increasing the channel height (for constant channel length) eventually results in a decrease of two in the number of vortices (Fig. 4).

We characterize the transition between the oscillating state and the Ceilidh dance (vortex lattice) regimes via the enstrophy $\boldsymbol{\omega}\left(\boldsymbol{r},t\right)\cdot\boldsymbol{\omega}\left(\boldsymbol{r},t\right)$, and measure activity in terms of the dimensionless activity number $A=\sqrt{\zeta h^{2}/K}$, which is motivated in more detail in the next section. At low activities, (in the unidirectional flow state; Movie 2(a)), the distribution of the enstrophy averaged across the channel $\Epsilon\left(x\right)=\left\langle\boldsymbol{\omega}\left(\boldsymbol{r},t\right)\cdot\boldsymbol{\omega}\left(\boldsymbol{r},t\right)\right\rangle_{y}$ is sharply peaked around small values (not shown). At moderate activities, the system enters the oscillating-flow state and the channel height-averaged enstrophy distribution spreads to higher values (Fig. 5(a); blue line). In this oscillating-flow state, the distribution is seen to decay monotonically from a maximum likelihood at its smallest non-zero value. At higher activity, however, the peak value $\Epsilon_{\text{peak}}$ shifts from the distributions smallest non-zero value to larger values (Fig. 5(a); red line). This peak-shift is due to the appearance of the vortex lattice and is a good indicator of the state of the flow structure (Fig. 5(b)).

A further sharp rise in the position of the enstrophy maximum in Fig. 5(b) marks the onset to the transition between the Ceilidh dynamics and meso-scale turbulence (Movie 2(d)). At this point the first signals of localized turbulent patches are identified with short mean lifetimes. Simulations of the following transition to system-spanning, steady-state meso-scale turbulence show that the transition belongs to the directed percolation universality class Doostmohammadi et al. (2016c). Experiments in microfluidic circuits suggest that the transition may be accompanied by a preferential direction depending on bacteria orientation at the walls Wioland et al. (2016).

To further investigate the nature of the dynamical transitions between steady states, we consider the transformation from laminar flow to Ceilidh dynamics as a function of activity number by measuring the system-averaged enstrophy $\left\langle\Epsilon\right\rangle_{x}$. As the activity is increased incrementally from zero, the system moves slowly through the dynamical steady states. When the activity is then lowered from the vortex-lattice state the Ceilidh dance state remains stable to lower activities, indicating hysteresis in this dynamic and intrinsically far-from-equilibrium flowing system (Fig. 5(c)).

III.2 Dynamical phase diagram

By mapping out the parameter space for which topological Ceilidh dynamics arise, we construct a phase diagram of dynamical steady flow regimes (Fig. 6). The competition between the activity driving a flow field and nematic elasticity resisting the associated deformation leads to the characteristic activity-induced length scale $\ell_{\zeta}\sim\sqrt{K/\zeta}$, which corresponds to the length scale over which the active stresses are accommodated by orientational elasticity Ramaswamy (2010); Thampi et al. (2013). However, confinement truncates the allowable range of $\ell_{\zeta}$ and screens interactions on separations $>h$. Thus, there is a competition between the activity-induced length scale $\ell_{\zeta}$ and the screening length $h$, as described by the dimensionless activity number $A=h/\ell_{\zeta}=\sqrt{\zeta h^{2}/K}$. For suspensions of microtubule-kinesin bundles with $\ell_{\zeta}\sim 100\mu\text{m}$ Sanchez et al. (2012) confined within channels of height $100\mu\text{m}\lesssim h\lesssim 1\text{cm}$ DeCamp et al. (2015), we expect a correspondingly broad range of activity numbers $A\sim 1-10^{2}$.

However, the activity number is not the only control parameter that affects the order parameter $\left\langle\Epsilon\right\rangle_{x}$ since the $+1/2$ disclination speed increases linearly with both activity and channel height (Fig. 2) and, therefore, does not depend directly on $A$. In addition to the dimensionless activity number $A$, we must account for the self-motility speed of the $+1/2$ disclinations Giomi et al. (2014) $v_{+1/2}\sim h\zeta/\eta$. The magnitude of the self motility must be judged against the characteristic velocity scale of the nematic liquid crystal Hemingway et al. (2016) $\nu_{Q}=\ell_{Q}\mathcal{E}_{Q}/\eta$, where $\mathcal{E}_{Q}$ is the characteristic nematic energy density scale (\mathcal{E}_{Q}=$$0.1) and $\ell_{Q}=\sqrt{K/\mathcal{E}_{Q}}$ is the equilibrium nematic persistence length, which is comparable to the defect core radius. The competition between these two velocity scales defines a characteristic self-motility number V=v_{+1/2}/\nu_{Q}\sim$$\ell_{Q}h\zeta/K. Considering microtubule-kinesin bundle suspensions, $+1/2$ disclinations move with a self-motility DeCamp et al. (2015) of $v_{+1/2}\approx 8\mu\text{m}\cdot\text{s}^{-1}$. An order of magnitude estimate of $\nu_{Q}\sim 10\mu\text{m}\cdot\text{s}^{-1}$ is found by estimating that the nematic elastic constant is $K\sim 1\text{pN}$, that the dense suspension at an oil/water interface has a viscosity roughly an order of magnitude larger than water, and that the nematic persistence length is approximately the observed defect core size $\ell_{Q}\sim 10\mu\text{m}$ DeCamp et al. (2015). Consequently, $V\sim 1$ is expected in such systems.

By plotting flow structures as a function of activity number $A$ and self-motility number $V$, the flow regimes are well separated into distinct dynamical steady states (Fig. 6). The non-flowing quiescent state is shown at the lower left and meso-scale turbulence at the upper right, with unidirectional flow, oscillating flow and Ceilidh dynamics found in between. The oscillating-flow state forms one boundary with the Ceilidh dynamics and the transition does not depend strongly on activity number $A$ at moderate $V$ (Fig. 6). Likewise, the first measurable occurances of active puffs denoting the early onset of meso-scale from Ceilidh dynamics is seen to occur at an activity number $A$ that depends only weakly on $V$.

The region of the dynamical phase diagram (Fig. 6) that we can access in our simulations is limited by numerical stability and it would be interesting to reach further into the low $A$/high $V$ and high $A$/low $V$ regimes. In particular the form of the phase boundaries suggest that a direct transition from the oscillating flow state might occur at higher self-motility. While we have simulations in the low $V$/high $A$ region of the diagram, it is difficult to clearly identify Ceilidh dynamics prior to the transition to meso-scale turbulence. While our previous work has shown that the transition from the Ceilidh dance to meso-scale turbulence belongs to the directed-percolation universality class Doostmohammadi et al. (2016c), this is suggestive of a secondary route to turbulence that goes directly from oscillating flows to chaotic flows. Furthermore, it will be interesting for future studies to explore the large $V$ and moderate $A$ region of phase space to see if there exists additional routes directly from oscillating to disordered active flow states.

III.3 Geometries that eliminate periodic boundary conditions

The simulations presented thus far are all for a channel geometry with periodic boundary conditions. To test whether Ceilidh dynamics arise in experimentally realizable geometries, we consider two bounded geometries: annuli (Fig. 7(a-b)) and closed channels (Fig. 7(c)). For consistency with the previous results using periodic boundary conditions, we keep anchoring at the boundaries, initial conditions and physical parameters the same, and set the gap between parallel confining walls to $h=25$ in all cases.

Annuli

As expected, annular geometries of two concentric bounding circles with small curvature (Fig. 7(a)) correspond closely to the results for Ceilidh dynamics in a channel with periodic boundary conditions: Non-motile $-1/2$ disclinations are evenly spaced around the inner and outer circular walls, while self-propelled $+1/2$ disclinations travel azimuthally around the annular gap with half moving clockwise and half counter-clockwise. The vortex lattice is preserved Theillard et al. (2017) and remains relatively unperturbed from the results of the channel geometry with periodic boundary conditions. Further increasing the curvature does not cause the Ceilidh dynamics to become unstable, though the vortex lattice skews substantially (Fig. 7(b)) and becomes comparable to the vortex fields observed in the vicinity of moderately confined rotors Thampi et al. (2016).

Closed channels

We next remove the periodicity altogether by considering a channel geometry that is closed by impermeable no-slip lateral walls (Fig. 7(c). We observe well-ordered Ceilidh dynamics in the centre of the channel, far from the lateral walls. Ceilidh dynamics emerge in the center of closed rectangular geometries with a sufficiently large aspect ratios ($\sim 10$). Near the lateral walls (within $\sim 4$ vortex sizes), a disorderly flow pattern occurs. In the immediate vicinity of the lateral walls disclination pairs are created, which feed $+1/2$ disclinations into the the ordered Ceilidh dancing in the centre of the channel. The resulting $-1/2$ disclination that arises from such pair-creation events drifts toward the adjacent lateral wall until it is annihilated by an outgoing $+1/2$ disclination.

III.4 Lattice defects

As a system-spanning ordered flow state, the Ceilidh dynamics state itself can possess various irregularities. We refer to these imperfections in the long-range structure as lattice defects, which should not be confused with the $\pm 1/2$ topological disclinations that occur in the continuous director field. These lattice defects arise when the system becomes trapped in a state in which a dancing pair of $+1/2$ disclinations is separated by other dancing-pairs. While the majority of irregularities are resolved as the steady state is approached, large elastic energy barriers can hinder rejoining of separated pairs and results in long-lived lattice defects. In channels with periodic boundary conditions, we do not observe additional disclinations (the number of topological disclination pairs is still equal to the number of vortices). Rather, there are only integer numbers of irregularities in the dynamic dance of the traveling $+1/2$ disclinations.

III.4.1 Broken-pair lattice defect

There is a finite probability that the distance between a pair of topological $+1/2$ disclinations does not reduce to $\lesssim h$ before another dancing-pair forms between them. The pair is then separated by a large energy barrier and we term this situation a broken-pair lattice defect (Fig. 8; Movie 3). In this broken-pair lattice defect, the number of disclinations moving leftward remains equal to the number moving rightward. Due to elastic repulsion between the like-signed disclinations and the confinement of the narrow channel, the separated disclinations cannot escape this imperfect state. However, it is important to recognize that no particular dancing-disclination pair permanently makes up the broken-pair lattice defect. Rather the Ceilidh dance ensures different dancing disclinations pass through the broken-pair lattice defect. In doing this, the disclinations shear the local vortex, causing it to become momentarily skewed and stronger (Fig. 8; Movie 3). Though the local vorticity field within the broken-pair lattice defect region is perturbed, the globally ordered topological dynamics of the rest of the system is not strongly affected.

III.4.2 Drift-lattice defect

In the broken-pair lattice defect, the number of disclinations on the leftward moving path remains equal to the number moving rightward. However, it is possible that upon the formation of the vortex lattice more disclinations reach one path than the other (Fig. 9(a); inset). As in the case of the broken-pair lattice defect, the energy cost of escaping this state is high and no disclination is observed to be able to change directions by switching dance-paths. In this case, two more disclinations move in one (spontaneously chosen) direction, and their activity drives a net flow $v_{\textmd{drift}}^{(1)}$ in that direction in addition to the zero-averaged vortex lattice flow. For this reason, we term this configuration a drift-lattice defect (Movie 4).

Although, increasing the activity is not seen to raise the probability of producing a drift-lattice defect in the Ceilidh state (Fig. 9(b)), increasing the length of the channel makes lattice defects more likely. For long enough channels, there is a non-zero probability that the number of lattice defects is greater than one. For each additional drift-lattice defect in the system there are two more $+1/2$ disclinations moving in the spontaneously chosen direction. In this way, the vortex lattice can be seen as a ground state and drift-lattice defects as excited states with quantization number $n$. The net flux of the system is seen to be cumulative, such that $v_{\textmd{drift}}^{(n)}=nv_{\textmd{drift}}^{(1)}$. The maximum number of drift excitations observed in our simulations is $n=3$ (Fig. 9(a); Movie 4(b)). We would expect that the probability of obtaining an ideally ordered system will decrease, that $n\geq 1$ will become more probable, and that a well defined number density of localized lattice defects will emerge as the length of the channel is increased.

The drift velocity $v_{\textmd{drift}}^{(1)}$ due to a single drift-lattice defect can be determined by estimating the net active force and viscous drag. The active force per unit volume $-\zeta\boldsymbol{\nabla}\cdot\boldsymbol{Q}$ is most significant near the $\pm 1/2$ disclinations, whose local deformation can be approximated as a solitary defect in bulk, $\boldsymbol{n}_{\pm 1/2}=\left[\cos\left(\pm\tfrac{1}{2}\phi\right),\sin\left(\pm\tfrac{1}{2}\phi\right)\right]$. By assuming that the ideal $+1/2$ disclinations are located at and oriented along the centerline of the channel and that the strong anchoring of the walls screen director deformations beyond a distance of $\sim h$, the net active force on the fluid due to a drift-lattice defect is found to be $F^{(1)}_{\text{act}}=c_{a}\eta v_{+1/2}$, where the numerical coefficient $c_{a}$ results from the precise integration area. For the case of $n$ drift lattice defects, the net force increases to $F^{(n)}_{\text{act}}=n\;F^{(1)}_{\text{act}}$. The drag that balances this active force is approximated by realizing that the train of vortices in the mid-channel region move en masse (Movie 4), while the flow in the near-wall region must viscously decrease from $\approx v_{\textmd{drift}}^{(n)}$ to zero at the no-slip boundary. Approximating the flow as simple shear and integrating over the near-wall regions gives the drag force to be $F_{\text{drag}}=-\left(c_{d}L/h\right)\eta v_{\text{drift}}^{(n)}$, where $c_{d}$ is a numerical drag coefficient whose details depend on the zero-averaged rotating vortex size and structure within the mid-channel region. These estimates predict a quantized drift velocity $v_{\text{drift}}^{(n)}\sim\left(\tfrac{h}{L}\right)nv_{+1/2}$ in agreement with Fig. 9(a).

IV Conclusions

It has been widely recognized that activity can drive collective motion in biological systems and the work presented here stresses the fact that activity can also produce highly ordered, but simultaneously dynamic, flow and topological fields. The topological Ceilidh dynamics that we report emerge as an intermediate state between spontaneous unidirectional flow and meso-scale turbulence for active nematics within a confining channel, when the characteristic vorticity length scale of meso-scale turbulence is comparable but smaller than the channel height.

While the time-varying topological disclination dynamics of the Ceilidh dance state were unexpected, the emergence of vortex lattices in active matter have been observed experimentally in motility assays of microtubles Sumino et al. (2012) and self-propelled particles Großmann et al. (2014); Jiang et al. (2017) with short range interactions and also numerically by hydrodynamic screening of activity-induced flows due to frictional damping Doostmohammadi et al. (2016b). Indeed, repeating vortex patterns have recently been reported in dense suspensions of bacteria confined within microfluidic channels Wioland et al. (2016), although the possibility of dynamically ordered disclination states was not discussed. In each of these examples, flows are screened. As such, the Ceilidh dance can be interpreted as arising due to the screening of the flow by the microfluidic channel walls.

Previous works have reported the emergence of long-range orientational (but not positional) order of active topological disclinations DeCamp et al. (2015); the current work advances this research direction towards dynamically ordered states and the effect of structural defects in the disclination patterning. Just as impurities, vacancies, dislocations and other defects play an important part in traditional condensed matter systems, the active-fluidic lattice defects discovered here may play critical but previously overlooked roles in precisely controlling the collective dynamics of biological systems, such as cellular monolayers or embryonic epithelial layers and will be relevant to the design of active microfluidic devices.

Acknowledgements.

Acknowledgments: This work was supported through funding from the ERC Advanced Grant 291234 MiCE and we acknowledge EMBO funding to TNS (ALTF181-2013). We thank Hugo Wioland for helpful discussions.

Appendix A Movie captions

Movie 1

Dancing of topological Ceilidh dynamics. (a) Vorticity field $\boldsymbol{\omega}\left(\boldsymbol{r},t\right)$ and streamlines of the highly ordered vortex lattice. The dynamically ordered state of disclinations within the channel is shown as $-1/2$ topological disclinations near the channel walls (magenta diamonds) and $+1/2$ disclinations transversing the mid-region (green dots). The system has zero net-topological charge and equal numbers of $+1/2$ disclinations travel to the left as to the right. (b) The scalar order field $q\left(\boldsymbol{r},t\right)$ and director field $\boldsymbol{n}\left(\boldsymbol{r},t\right)$ for the same simulation as Movie 1(a).

Movie 2

Vorticity field $\boldsymbol{\omega}\left(\boldsymbol{r},t\right)$ and and topological defects for the four types of flow states observed in simulations of active nematic fluids confined within two dimensional micro-channels. Magenta diamonds denote $-1/2$ topological disclinations and self-motile $+1/2$ disclinations are represented as green dots. (a) Unidirectional flow in a spontaneously chosen direction along the channel occurs at low activity. (b) Oscillating but unidirectional laminar flow, with temporary traveling vorticity patterns that intermittently dissipate then re-form. (c) Ceilidh dynamics with a well-defined vortex lattice at intermediate activity. (d) The chaotic flow state of meso-scale turbulence with the continual creation and annihilation of topological disclinations at the highest activity numbers.

Movie 3

Broken-pair lattice defect. Vorticity field $\boldsymbol{\omega}\left(\boldsymbol{r},t\right)$ and topological disclinations ($-1/2$ magenta diamonds; $+1/2$ green dots) for topological Ceilidh dynamics with a single broken-pair lattice defect, in which the number of disclinations moving leftward remains equal to the number moving rightward. The lattice defect locally skews the vorticity field about the two separated $+1/2$ disclinations.

Movie 4

Drift-lattice defects for a single defect (top) and for three defects. (a) Vorticity field $\boldsymbol{\omega}\left(\boldsymbol{r},t\right)$ and topological disclinations ($-1/2$ magenta diamonds; $+1/2$ green dots) for Ceilidh dynamics with a single drift-lattice defect ($n=1$), in which two more disclinations are moving leftward than rightward. A vertical dashed line is drawn across the channel to mark a static point against which the leftward drift of the vortex lattice can be identified. Alternatively, an arrow points to a single $-1/2$ defect so that the distance travelled can be assessed at all times. (b) The same system as Movie 4(a) but for an example simulation with three drift-lattice defects ($n=3$). Comparison with Movie 4(a) shows that the drift velocity is substantially faster than when only a single drift-lattice defect is present.