Open Access Article
This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

Chiral stresses in nematic cell monolayers

Ludwig A. Hoffmann a, Koen Schakenraad ab, Roeland M. H. Merks bc and Luca Giomi *a
aInstituut-Lorentz, Leiden University, P.O. Box 9506, 2300 RA Leiden, The Netherlands. E-mail: giomi@lorentz.leidenuniv.nl
bMathematical Institute, Leiden University, P.O. Box 9512, 2300 RA Leiden, The Netherlands
cInstitute of Biology, Leiden University, P.O. Box 9505, 2300 RA Leiden, The Netherlands

Received 13th September 2019 , Accepted 28th November 2019

First published on 2nd December 2019


Abstract

Recent experiments on monolayers of spindle-like cells plated on adhesive stripe-shaped domains have provided a convincing demonstration that certain types of collective phenomena in epithelia are well described by active nematic hydrodynamics. While recovering some of the hallmark predictions of this framework, however, these experiments have also revealed a number of unexpected features that could be ascribed to the existence of chirality over length scales larger than the typical size of a cell. In this article we elaborate on the microscopic origin of chiral stresses in nematic cell monolayers and investigate how chirality affects the motion of topological defects, as well as the collective motion in stripe-shaped domains. We find that chirality introduces a characteristic asymmetry in the collective cellular flow, from which the ratio between chiral and non-chiral active stresses can be inferred by particle-image-velocimetry measurements. Furthermore, we find that chirality changes the nature of the spontaneous flow transition under confinement and that, for specific anchoring conditions, the latter has the structure of an imperfect pitchfork bifurcation.


I. Introduction

Multicellular systems of prokaryotes, such as suspensions of planktonic bacteria, have historically played a pivotal role in the development of the hydrodynamics of active fluids, since the early work of Batchelor on the stress distribution in suspensions of microswimmers.1 By contrast, multicellular systems of eukaryotes have entered only recently into the realm of active hydrodynamics, following a number of inspiring experimental works on epithelial and mesenchymal cell layers and tissues (see e.g.ref. 2–8). Among these, monolayers of clotured spindle-like cells, such as NIH 3T3 mouse embryo fibroblasts,4 murine neural progenitor cells (NPCs),6 human bronchial epithelial cells (HBEC),7 Retinal Pigment Epithelial (RPE1) cells8 and C2C12 mouse myoblasts8 represent an especially promising class of model systems, because of their connection with active nematic liquid crystals (see e.g.ref. 9).

First identified as a broken symmetry in certain types of cell cultures,10 and later exploited to decipher their static11 and dynamical properties,2,4–8 nematic order has surged as one of the central themes in collective cell dynamics. In layers of spindle-like cells, where the local orientation can be unambiguously identified, nematic order is marked by the presence of ±1/2 disclinations,4,6i.e. point-like singularities about which the average cellular orientation rotates by ±π. Consistently with the predictions of active nematic hydrodynamics,12,13 these defects self-propel and pairwise annihilate until cell crowding freezes the system into a jammed configuration. Before dynamical arrest, the collective motion of the cells gives rise to a decaying turbulent flow at low Reynolds number, whose statistics, spatial organization and spectral structure are in exceptional agreement with the hydrodynamic picture14 (but see ref. 15 for an alternative theoretical picture based on glassy dynamics).

Another remarkable demonstration of active hydrodynamic behavior in eukaryotic cell layers has been recently reported by Duclos et al., upon confining spindle-like RPE1 and C2C12 cells within adhesive stripe-shaped domains.8 Depending on the width of the stripe the system was found either in a stationary state, with the cells parallel to the longitudinal direction of the confining region, or in a collectively flowing state characterized by a spontaneous tilt of the cells toward the center of the stripe. The latter picture, often referred to as spontaneous flow transition, had been anticipated for over a decade by Voituriez et al.16 and represents one of the hallmarks of active liquid crystals. While confirming this seminal prediction, however, Duclos et al. have also highlighted a number of unexpected features that could be ascribed to the existence of chirality over length scales larger than the typical size of a cell.

The notion of chirality is not new in active matter and has been theoretically explored well before the interest around collective cell dynamics in eukaryotes had started to blossom. Fürthauer et al., for instance, demonstrated that microscopic torque dipoles, such as those arising from rotating molecular motors or flagella, give rise to antisymmetric stresses and angular momentum fluxes, which, in turn, drive rotating flows and other chiral patterns on the large scale.17,18 More recently, Banerjee et al. showed that rotational motion at the microscopic scale further enriches the spectrum of hydrodynamical behaviors or chiral active fluids by giving rise to non-dissipative “odd” viscosity,19 analogous to that found in quantum Hall fluids.20 Whereas undoubtedly interesting and relevant for a broad class of biological and synthetic systems, these mechanisms appear however unsuited to account for the chirality observed in the experiments by Duclos et al., because of the manifest lack of rotational motion at the scale of individual cells.

In this article, we show that macroscopic chirality can arise in nematic cell layers as a consequence of a misalignment between the cell's local orientation and active forces, even in the absence of microscopic rotational motion (Section II). Collectively, this gives rise to a chiral and yet symmetric stress tensor that, in two dimensions, complies with the symmetries of the nematic phase. Next, we explore the effect of such a chiral stress on the active flow generated by ±1/2 disclinations and identify a characteristic signature of chirality from which the ratio between chiral and non-chiral active stresses can be experimentally estimated (Section III). Finally, following Duclos et al.,8 we investigate the hydrodynamic stability of a chiral nematic cell monolayer confined on adhesive stripes and subject to various boundary conditions and classify all possible scenarios arising from the interplay between the geometry of the confining region, the extensile/contractile stresses and chirality (Section IV).

II. Chiral stresses

Let us consider a two-dimensional volume element and let n be the nematic director representing the average direction of the enclosed cells (Fig. 1a). The most generic form of the stress tensor associated with such a volume element can be expressed in the basis of the nematic director n and its normal vector n, namely:
 
σa = σnn + σnn + τ(nn + nn).(1)
Here σ and σ represent the stresses experienced by the volume element in the direction parallel and perpendicular to n, whereas τ is the shear stress. By construction, σa is invariant under the transformation n → −n, thus is consistent with the symmetry of the nematic phase, but is evidently chiral for any non-zero τ value, as the stress distribution depicted in Fig. 1a does not coincide with its mirror image. Furthermore, taking nn = Inn, with I the identity tensor, allows one to cast the active stress as the sum of an isotropic contribution, equivalent to an active pressure, a deviatoric term, common to all active nematic liquid crystals regardless of their chirality,21,22 and a chiral term, namely:
 
image file: c9sm01851d-t1.tif(2)
where:
 
image file: c9sm01851d-t2.tif(3)
Microscopically, the chiral stress τ might originate from the fact that the force exerted by an individual cell is tilted with respect to the cell orientation. To illustrate this concept let us consider an individual cell whose major and minor axes are parallel to the unit vectors νc and νc (Fig. 1b). Following Lau and Lubensky,23 one can express the stress tensor as σa = Σcdcδ(rrc), where the index c runs over all the cells in the system, rc is the position of the c-th cell and
 
image file: c9sm01851d-t3.tif(4)
is a force–dipole tensor. In eqn (4), Rc is the distance between the cell's surface and center of mass, −dAfc the force exerted by the cell's area element dA on the surrounding medium and the integral is calculated over the surface Σc of the c-th cell. A complete derivation can be found in ref. 23. For an effectively two-dimensional system, such as the one considered here, dA = d[small script l]h, with h the thickness of the cell along the z-direction, here assumed to be uniform throughout the sample, and [small script l] the one-dimensional arc-length.

image file: c9sm01851d-f1.tif
Fig. 1 Schematic representation of stresses and forces in cellular nematic monolayers. (a) A volume element whose faces are conventionally oriented in the direction of the nematic director n and its normal vector n. In the most general setting, the volume element is subject to three independent stresses: the two normal stresses σ and σ and the shear stress τ. (b) The chiral stress τ arises at the microscopic level when the force exerted by an individual cell is neither parallel nor perpendicular to the cell axis.

Now, the magnitude of the active stresses Pa, α, and τ depend exclusively upon the distribution of the forces exerted by the cells along their contour. The simplest approximation of the force density field fc consists then of a dipole of the form:

 
fc = Fcδ(Rcaνc) − Fcδ(Rc + aνc),(5)
where a is the cell's major semi-axis (Fig. 1b and ref. 22). To make progress, we express the cellular forces in the {νc,νc} basis, namely Fc = Fνc + Fνc, with F and F the longitudinal and transverse components of the force exerted by the c-th cell, here assumed for simplicity uniform throughout the system. Replacing this in eqn (5) and coarse-graining σa over the length scale of a volume element Ω(r) centered at r (Fig. 1a) yields:
 
σa = [2Fνcνc〉 + Fνcνc + νcνc〉],(6)
where 〈⋯〉 is the spatial average within Ω(r) and ρ is the local cell number density. Finally, calculating the averages and comparing with eqn (2) readily gives:
 
Pa = −aρF, α = 2aρSF, τ = aρSF,(7)
where S = 2〈|νc·n|2〉 − 1 is the local nematic order parameter in two dimensions. As expected, in the absence of nematic order (i.e. S = 0), the active forces exerted by the cells result exclusively in an effective pressure, while both the deviatoric and chiral stress vanishes identically. When nematic order is not uniform throughout the system, hence S varies in space, the active stress tensor of eqn (2) is more conveniently expressed in terms of tensor order parameter Qij = S(ninjδij/2), namely:
 
σaij = −Paδij + α0Qij − 2τ0εikQkj,(8)
where α0 = 2aρF and τ0 = aρF are constants independent on the nematic order parameter S and εij is the two-dimensional Levi-Civita tensor (i.e. εxx = εyy = 0 and εxy = −εyx = 1).

Some comments are in order. Although eqn (5) is only a rudimental approximation of the force field generated by an irregularly-shaped cell, considering a more involved force distribution does not change the qualitative picture with respect to the emergence of chiral stresses, as long as this is asymmetric with respect to the cell's longitudinal direction. To illustrate this concept, we discuss in Appendix A the case of a quadrupolar force distribution. Whereas the exact origin of this asymmetry is beyond the scope of the present article, the biophysical literature is not scarce of examples where chirality can be detected at the single-cell level. For instance, various mammalian cells, when plated on micropatterns, can break the left-right symmetry by suitably positioning their internal organelles with respect to the cell body.24 Analogously, chirality can emerge at the scale of the entire cell from the self-organization of the actin cytoskeleton.25,26 The broken symmetry can furthermore propagate over the mesoscopic scale and bias the cell's collective migratory motion.27

III. Defect motion

The motion of ±1/2 disclinations has become a hallmark of active nematic liquid crystals. As it was theoretically predicted12,13 and experimentally verified in both microtubules-kinesin28–31 and actin-myosin suspensions,32 ±1/2 disclinations in active nematics drive a flow, whose structure and speed is strictly related to the geometry of the defect. For +1/2 disclinations, in particular, this flow has a Stokeslet structure that results in a propulsion of the defect in the direction of its symmetry axes p or −p depending on whether the system is contractile or extensile (Fig. 2 and ref. 33). Exceptionally, the same mechanism has been experimentally identified in various types of cell cultures, including spindle-shaped NIH 3T3 mouse embryo fibroblasts,4 murine neural progenitor cells (NPCs),6 Madin Darby canine kidney cells (MDCKs)5 and human bronchial epithelial cell (HBEC).7 Unlike in suspensions of rod-like cytoskeletal filaments, however, the direction of motion of +1/2 disclinations is not necessarily parallel to ±p and also the reconstructed flow around the defects is less symmetric than that found in cytoskeletal suspensions. Whereas it is not unlikely for such a feature to originate from statistical errors, here we demonstrate that chiral active stresses affect the dynamics of +1/2 defects precisely by rotating their direction of motion with respect to p. In fact, the angle between p and the direction of motion of the defects could be used to indirectly measure the relative magnitude of the chiral stress τ compared to the deviatoric stress α. This phenomenon shares some similarities with recent results by Maitra and Lenz about the dynamics of +1/2 disclinations in rotating active nematics.34
image file: c9sm01851d-f2.tif
Fig. 2 Pressure (a and b) and velocity (c and d) fields in proximity of ±1/2 defects obtained from the analytical solutions of eqn (9) for an extensile chiral active nematic with α = 4τ < 0. The configuration of the director is indicated by white lines for the case of +1/2 (a) and −1/2 (b) defects. For +1/2 defects, the velocity field at the core, thus the direction of motion of the defect, is tilted by an angle θtilt = arctan(1/2) ≈ 27° with respect to the defect polarity direction p = [x with combining circumflex]. For both chiral and achiral +1/2 defects, the pressure is anisotropic and larger toward the direction of motion of the defects.

An analytical approximation of the flow driven by the active stresses in the surrounding of a disclination can be obtained by solving the incompressible Stokes equation with a body force resulting from the active stress associated with an isolated ±1/2 defect. Namely:

 
η2v + fa± = ∇P, ∇·v = 0,(9)
where η is the shear viscosity of the tissue, here assumed isotropic for simplicity, and fa± = ∇·σa is the body force arising from spatial variations of the active stress of eqn (2) in the presence of a ±1/2 defect. This approach, introduced in ref. 13, does not account for feedback of the flow on the orientation of the director and hence the structure of the defect. Nevertheless it provides a simple and faithful approximation of the defective flow as well as an estimate of the defect velocity. A generic solution of eqn (9) can be expressed in the form v = v± + v0, where v0 is a solution of the homogeneous Stokes equation and is dictated by the boundary conditions and v± is given by:
 
image file: c9sm01851d-t4.tif(10)
where G is the two-dimensional Oseen tensor given by:
 
image file: c9sm01851d-t5.tif(11)
with image file: c9sm01851d-t6.tif a constant set by the boundary conditions. Analogously, taking the divergence of eqn (9) allows to express the pressure field as the solution of the following Poisson equation:
 
2P± − ∇·fa± = 0.(12)
Now, in the presence of a disclination of strength s = ±1/2 located at the origin of the (x,y)-plane and oriented in the direction p = (cos[thin space (1/6-em)]ψ,sin[thin space (1/6-em)]ψ), the nematic director n = (cos[thin space (1/6-em)]θ,sin[thin space (1/6-em)]θ) has local orientation θ = + (1 − s)ψ, with ϕ = arctan(y/x) the polar angle.33 The body force fa± is then readily found in the form:
 
image file: c9sm01851d-t7.tif(13)
where p = (−py,px) and r is the distance form the defect core. The effect of chirality is most dramatic for +1/2 disclinations. In achiral active nematics, τ = 0 and the active force fa+ is purely longitudinal. For non-zero τ, fa+ acquires a transverse component resulting in a tilt in the direction of motion of the defect by an angle
 
image file: c9sm01851d-t8.tif(14)
with respect to the orientation p (Fig. 2). As anticipated, eqn (14) can in principle be used in combination with experimental reconstruction of defect trajectories in order to estimate the relative magnitude of the chiral and deviatoric stresses in nematic cell monolayers.

To calculate the flow velocity in proximity of a ±1/2 defect, we set, without loss of generality, ψ = 0, thus p = [x with combining circumflex] and p = ŷ and we assume the defect at the center of a circular domain of radius R. The exact velocity of the flow at the boundary of such a domain, hence the homogeneous solution v0, is not relevant for the purpose of this discussion. In practice, this will be determined by the chemistry of the substrate and the possible presence of other topological defects in the same region.4 Under these assumptions, and carrying out algebraic manipulations as those in ref. 13, the flow velocity caused by a ±1/2 defect can be found from eqn (10), (11) and (13) in the form:

 
image file: c9sm01851d-t9.tif(15a)
 
image file: c9sm01851d-t10.tif(15b)
These flows are illustrated in Fig. 2c and d for a specific choice of the angle θtilt. The corresponding pressure field is readily found from eqn (12) with ∇·fa+ = −(α[thin space (1/6-em)]cos[thin space (1/6-em)]ϕ + 2τ[thin space (1/6-em)]sin[thin space (1/6-em)]ϕ)/(2r2) and ∇·fa = 3(α[thin space (1/6-em)]cos[thin space (1/6-em)]3ϕ + 2τ[thin space (1/6-em)]sin[thin space (1/6-em)]3ϕ)/(2r2). This yields:
 
image file: c9sm01851d-t11.tif(16a)
 
image file: c9sm01851d-t12.tif(16b)
where P(0)± are harmonic functions depending on the boundary conditions. Interestingly, the pressure field given by eqn (16) is independent on the distance from the defect core, but varies with the angle ϕ and, depending on the sign of the stresses α and τ, is maximal or minimal at specific directions relatively to the polarity vector p. For instance, for extensile systems (i.e. α < 0) with negligible chirality (i.e. τ ≈ 0), eqn (16a) predicts a pressure maximum in the −p direction, thus toward the “head” of the comet-like structure characteristic of +1/2 defects. Despite our analysis revolving around incompressible systems, we can expect this trend to persist in the presence of small spatial variations of the density field ρ. In this case, and under the assumption that ρP, we expect a higher density of extensile (contractile) cells in the front (back) of +1/2 defects, consistent with the experimental observation of Kawaguchi et al.6 Finally, the lack of spatial dependence in eqn (16) originates for the specific parametrization of the orientation field θ in proximity of the defects, i.e. θ(ϕ) = ±ϕ/2. In a more realistic setting, θ(r,ϕ) = ±ϕ/2 + θfar(r,ϕ), where θfar is a non-singular function that determines the far-field configuration of the nematic director and vanishes at the defect core. Accounting for the far-field configuration of the director results into a dependence of the pressure on the distance from the defect core.

In summary, the presence of a symmetric chiral active stress, such as that embodied by the parameter τ, affects the flow generated by ±1/2 disclinations by stretching and rotating the velocity field in the surrounding of the defects (Fig. 2c and d). Most prominently, this results in a tilt in the direction of motion of +1/2 defects: i.e.vself = v+(r = 0) = R/(4η)(αp + 2τp). Thus +1/2 defects self-propel at an angle θtilt with respect to their orientation p [see eqn (14)]. Such an angle could in principle be measured in experiments on two-dimensional cell cultures, thus providing a direct measurement of the relative magnitude of the chiral stress. The same behavior has been reported in the case of actively rotating +1/2 defects, in the limit of vanishing angular velocity.34

IV. Spontaneous flow on adhesive stripes

In this section we revisit a classic problem of the hydrodynamic stability and spontaneous flows of active nematics in a quasi-one-dimensional channel (Fig. 3). First discussed in a seminal paper by Voituriez et al.,16 later elaborated by many others35–37 and recently observed in experiments on cell monolayers8 and suspensions of microtubules and kinesin,38 this phenomenon consists of a continuous transition between a stationary and uniformly oriented configuration to a state characterized by a spontaneous distortion of the nematic director coupled to an internally driven shear flow. The transition, in many aspects similar to the Fréedericksz transition in passive liquid crystals,39 results as a consequence of two different mechanisms. First, a distortion of the nematic director drives a shear flow as illustrated by eqn (9) for the time-independent case; second, the nematic director rotates in a shear flow. As a consequence, when the hydrodynamic torque driven by the active stresses outweighs the elastic restoring torque, the uniformly oriented configuration becomes unstable to splay or bending deformations, depending on the sign of the active stress α and other material parameters. Roughly speaking, this occurs when the active length scale image file: c9sm01851d-t13.tif, defined by the ratio between the passive and active torques, with K the Frank elastic constant, becomes comparable to the width L of the channel.14
image file: c9sm01851d-f3.tif
Fig. 3 Schematic representation of the channel of infinite length in x-direction and width L in y-direction.

In the following, we extend and generalize the theoretical analysis by Duclos et al.8 by considering various experimentally relevant scenario in terms of boundary anchoring and flow. The hydrodynamic equations governing the dynamics of an incompressible (i.e. ∇·v = 0) active nematic liquid crystal are given by:16,36,37

 
image file: c9sm01851d-t14.tif(17a)
 
image file: c9sm01851d-t15.tif(17b)
where D/Dt = ∂t + v·∇ is the material derivative, uij = (∂ivj + ∂jvi)/2 and ωij = (∂ivj − ∂jvi)/2 are respectively the strain-rate and vorticity tensor and h = −δFn is the so-called molecular field, governing the relaxational dynamics of the nematic director, with F the Frank free energy. In one-elastic-constant approximation, the latter is simply given by image file: c9sm01851d-t16.tif and h = K2n. The material parameters λ and γ are respectively the flow-alignment parameter, controlling the tendency of the nematic director to rotate in a shear flow, and the rotational viscosity of the nematic fluid. In eqn (17b), the pressure P includes the active pressure Pa given in eqn (2) and (7), but, due to incompressibility, is not an independent field and σe is the elastic stress resulting from the departure of the director configuration from the ground state of the Frank free energy. Although it does affect the onset of the spontaneous flow instability, the latter is often unimportant for the phenomenology of active nematics and will be neglected here for simplicity.

We solve eqn (17) in a rectangular channel which is infinitely long in x-direction and has width L in y-direction (Fig. 3). Because the channel is invariant for translations along the longitudinal direction, we assume for simplicity n and v to be independent of x, both while stationary and spontaneously flowing. Furthermore, incompressibility and mass conservation demand the y-component of the velocity to vanish identically, i.e. vy = 0. Thus taking n = (cos[thin space (1/6-em)]θ,sin[thin space (1/6-em)]θ) and v = (vx,0), the eqn (17) reduce to:

 
image file: c9sm01851d-t17.tif(18a)
 
image file: c9sm01851d-t18.tif(18b)
whereas from the y-component of eqn (17b) we obtain the following expression for the pressure field:
 
image file: c9sm01851d-t19.tif(19)
with P0 a constant. Next, we look for stationary solutions of eqn (18) subject to different boundary conditions in terms of orientation of the nematic director at the boundary and whether or not the cells are allowed to slide along the channel walls, thus: θ = θ(y), vx = vx(y). In order for the fluid to be stationary, the shear stress must be uniform across the channel. Thus, from eqn (18b):
 
image file: c9sm01851d-t20.tif(20)
Solving eqn (20) with respect to ∂yvx and substituting this into eqn (18b) yields a single homogeneous equation for θ, namely:
 
image file: c9sm01851d-t21.tif(21)
Before analyzing specific cases in detail, we consider the generic situation in which the nematic director is anchored to the channel walls at an arbitrary angle θ0 ≠ arccos(1/λ)/2, thus:
 
θ(0) = θ(L) = θ0.(22)
In the following, we will separately analyze the scenarios in which the cells are stationary at the boundary of the channel (Section IV.A) and when, on the other hand, they are able to slide while keeping their orientation fixed (Section IV.B).

Our analysis is complemented by numerical solutions of eqn (18) with various boundary conditions. For this purpose we rescale time by the viscous time scale τν = ρL2/η, length by the channel width L and stress by the viscous stress scale σν = ρL2/τν2, i.e. tt/τν, yy/L, σσ/σν. All the other quantities in Fig. 4 and 5 are rescaled accordingly.


image file: c9sm01851d-f4.tif
Fig. 4 Bifurcation diagram of the spontaneous flow transition obtained from numerical (dots) and analytical (lines) solutions of eqn (18) for K/γ = L = 1, λ = −0.5 and various τ values. (a) No-slip boundary conditions and parallel anchoring (see Section IV.A.1). The chiral stress τ does not influence the critical α value, but weakly affect the post-transitional configuration of the nematic director. (b) Stress-free boundary conditions and special boundary anchoring θ(0) = θ(L) = −θtilt/2 (see Section IV.B.3). The chiral active stresses, embodied by the parameter τ, explicitly break the clock-counterclockwise symmetry of the lowest free-energy configuration rendering the pitchfork bifurcation “imperfect”. In this case, only one of the two branches of the bifurcation diagram is connected to the trivial solution, which may then be the only one observed experimentally.

image file: c9sm01851d-f5.tif
Fig. 5 Numerical solutions of eqn (18) in the dimensionless quantities defined above and for K/γ = 1, λ = −0.5 as well as different values of τ and α for three sets of boundary conditions. (a) and (b) are numerical solutions for the eqn (18) with parallel anchoring as well as no-slip boundary conditions [panel (a) and analyzed analytically in Section IV.A.1] and stress-free boundary conditions [panel (b) and Section IV.B.1], respectively. (c) Numerical solutions for stress-free walls and the stationary solution θ0 = −θtilt/2 = −arctan(2τ/α)/2 imposed at the boundary, see Section IV.B.3.

A. No-slip boundary conditions

In this section, we consider the case in which the cells are unable to slide along the boundary of the channel, i.e. vx(0) = vx(L) = 0, which, in turn, experience a non-vanishing stress σxy ≠ 0 resulting from the cellular forces. In this case, θ(y) = θ0 is always a trivial solution of eqn (21) with the boundary condition given by eqn (22), and the monolayer admits a stationary and uniformly aligned configuration. Because of the internal stresses, however, such a uniform state can become unstable with respect to splay or bending deformations for sufficiently large active stresses or channel widths. To illustrate this point, we take θ(y) = θ0 + δθ(y), with δθ(0) = δθ(L) = 0, and linearize eqn (21) about δθ = 0. This yields:
 
y2δθ + q2δθ = 0,(23)
where we have introduced the constant:
 
image file: c9sm01851d-t22.tif(24)
The solution of eqn (23) is readily found to be:
 
δθ = C[thin space (1/6-em)]sin(qy), qL = nπ,(25)
with a constant C and n[Doublestruck Z] an arbitrary integer. By virtue of eqn (20), the corresponding velocity field is given by:
 
ηyvx = δθ(2τ[thin space (1/6-em)]sin[thin space (1/6-em)]2θ0α[thin space (1/6-em)]cos[thin space (1/6-em)]2θ0).(26)
whose solution with no-slip boundary conditions is given by:
 
image file: c9sm01851d-t23.tif(27)
with m[Doublestruck Z] another arbitrary integer. Upon comparing eqn (25) and (27) we find that the first mode to be excited is (n,m) = (2,1), thus the trivial solution θ(y) = θ0 becomes unstable when q = qc = 2π/L or, equivalently, when L = Lc = 2π/q. To be more specific, we consider, in the following, two practically relevant cases, where θ0 = 0 (parallel anchoring) and θ0 = π/2 (homeotropic anchoring).
1. Parallel anchoring. If the nematic director is parallel to the channel walls θ0 = 0, thus:
 
image file: c9sm01851d-t24.tif(28)
As in non-chiral active nematics, the instability is triggered by splay deformations (i.e. transverse to the nematic director) and uniquely depends on the non-chiral active stress α. Furthermore, as in non-chiral active nematics,13 such a splay instability affects flow-aligning systems (i.e. λ > 1) in the presence of extensile active stresses (i.e. α < 0), and flow-tumbling systems (i.e. λ < 1) in the presence of contractile active stresses (i.e. α > 0). A critical α value is readily found in the form: αc = 8π2ηK/[γL2(1 − λ)].

At the onset of the transition, the constant C can be calculated upon expanding eqn (21) up to the third order in δθ, Then, using the solution of the linearized equation yields a cubic equation in C. Solving the latter gives (see Appendix B):

image file: c9sm01851d-t25.tif

The spontaneous flow instability consists, therefore, of a standard pitchfork bifurcation whose relevant fields, θ and vx, scale like (LLc)1/2 at criticality, see Fig. 4a. Despite the fact that the chiral stress τ does not affect the instability of the stationary state, it leaves a clear signature on the post-transitional behavior of the flowing monolayers. This can be seen in Fig. 5a, showing numerical solutions of eqn (18) in the flowing state for various α and τ values. The most prominent effect of chirality, in this case, is evidently to render both the distortion of the nematic director and the associated flow asymmetric with respect to the channel centerline.

2. Homeotropic anchoring. The case in which the nematic director is perpendicularly anchored to the channel walls, θ0 = π/2, yields:
 
image file: c9sm01851d-t26.tif(29)
In this case the instability is triggered by bending deformations (i.e. parallel to the nematic director). In contrast with the scenario of Section IV.A.1, flow-aligning systems are unstable in the presence of extensile active stresses, whereas strongly flow-tumbling systems (i.e. λ < −1), are unstable in the presence of contractile active stresses. The critical α value is readily found in the form: αc = −8π2ηK/[γL2(1 + λ)].

To conclude this section, we observe that for both parallel and homeotropic anchoring, the spontaneous flow instability crucially relies on the flow-alignment behavior of the system, governed by the phenomenological parameter λ. As for molecular liquid crystals, where λ depends upon the molecules shape and interactions, we expect the flow-alignment parameter to be affected by the cellular shape, which in turn is not fixed, and by the passive and active processes underlying the cell–cell and cell–substrate interactions.

B. Stress-free boundary conditions

In this subsection we consider the scenario in which the cells are allowed to slide along the boundary, while keeping a fixed orientation θ0 with respect to the channel walls. As in the adhesive stripes used in ref. 8, the channel walls do not comprise a real physical barrier, but represent instead the interface between two regions of the substrate with different coating. As the walls now do not exert any force on the cells σxy(0) = σxy(L) = 0. Mechanical equilibrium [i.e.eqn (20)] thus implies σxy = 0 everywhere.

The most striking difference with respect to the case discussed in Section IV.A, as well as the most prominent consequence of the chiral stress τ, is that the stationary and uniformly aligned configuration (i.e. θ = θ0 and vx = 0) is not a trivial solution of eqn (21) for non-vanishing α and τ values, unless the chiral and non-chiral active stresses cancel each other identically. Before considering this latter case (see Section IV.B.3), we find approximated expressions for the local orientation θ and the velocity vx in the limits in which the active stresses are either very small or very large.

For very large active stresses, the active terms in eqn (21) overweight the elastic term on the left hand side. As a consequence, the equilibrium configuration of the nematic monolayer consists of a region in the bulk of the channel where the cells have uniform orientation θ = −θtilt/2 = −arctan(2τ/α)/2 and two boundary layers, whose size is roughly image file: c9sm01851d-t27.tif, with σ0 = (α/2)sin[thin space (1/6-em)]2θ0 + τ[thin space (1/6-em)]cos[thin space (1/6-em)]2θ0, where the director interpolates between the bulk and boundary orientation (see Fig. 5b). This phenomenon closely resembles flow-alignment in nematics (see e.g.ref. 40) with −θtilt/2 playing the role of the Leslie angle θL = arccos(1/λ)/2. Whereas passive flow-alignment, however, requires λ > 1 (e.g. flow-aligning nematics), such an active flow-alignment occurs at any finite value of α and τ, provided the elastic boundary layer is sufficiently small to have a clear distinction between bulk and boundary alignment.

For small α and τ values, we can postulate that the nematic director will depart only slightly from its orientation at the boundary. Thus, taking again θ(y) = θ0 + δθ(y) and linearizing eqn (21) around δθ = 0, we obtain:

 
y2δθ + q2θ + δθ0) = 0,(30)
with q2 given, as before, by eqn (24) and:
 
image file: c9sm01851d-t28.tif(31)
For a general anchoring angle θ0, a solution of eqn (30) with boundary conditions δθ(0) = δθ(L) = 0 is given by:
 
image file: c9sm01851d-t29.tif(32)
The associated velocity field can then be found from a direct integration of the linearized equation:
 
image file: c9sm01851d-t30.tif(33)
The lack of a boundary condition for eqn (33) can be compensated with a global constraint on the total momentum, namely image file: c9sm01851d-t31.tif.

In the following, we provide explicit approximated expression for the velocity field in the special cases where θ0 = 0 (parallel anchoring) and θ0 = π/2. Furthermore, we will investigate the stability of the trivial solution of eqn (21) obtained when θ0 is such that the chiral and non-chiral active stresses cancel each other identically.

1. Parallel anchoring. For θ0 = 0 eqn (28) and (31) yield:
 
image file: c9sm01851d-t32.tif(34)
The corresponding velocity field is then readily obtained by integrating eqn (33). This gives:
 
image file: c9sm01851d-t33.tif(35)
where the wave number q is that given in eqn (28). Numerical solutions for this case are displayed in Fig. 5b for various α and τ values.

We stress that, whereas the flowing configurations resulting from the instability of the stationary state are left-right and clock-counterclockwise symmetric (i.e. the cells are equally likely to flow toward the negative or positive x-direction and, correspondingly, to tilt clock- or counterclockwise, see Section IV.A.1), in this case the direction of the tilt as well as that of the flowing monolayer is set by the signs of the constants α and τ.

2. Homeotropic anchoring. For θ0 = π/2 from eqn (28) and (31) we find that the amplitude δθ0 is given, once again, by eqn (34). Thus, the expressions for δθ and vx are formally identical to those given in eqn (32) and (35), but with wave number q as given in eqn (29).
3. Stationary solution. To conclude this subsection, we consider a special situation where the orientation of the cells at the boundary is fixed, as before, but such that the chiral and non-chiral stresses cancel each other identically. Thus: θ0 = −θtilt/2 = −arctan(2τ/α)/2. In this case the orientation of the nematic director in the bulk of the channel, determined by the balance between the chiral and non-chiral active stress, is equal to that at the boundary. As a consequence, the boundary layer described in Section IV.B disappears and the system can achieve a stationary and uniformly aligned configuration. As those described in Section IV.A, however, the latter is unstable for sufficiently large active stresses or channel width.

Using the same algebraic manipulations adopted in Section IV.A, one can show that the perturbation δθ is again of the form given in eqn (25) with:

 
image file: c9sm01851d-t34.tif(36)
Analogously, the velocity is given by eqn (27), but with no constraint on the phase qL, because of the stress-free boundary conditions. As a consequence, the first mode to be excited is n = 1, thus the stationary state becomes unstable when q = qc = π/L or, equivalently, when L = Lc = π/q. Some numerical solution of eqn (18), in this regime, are shown in Fig. 5c.

Notably, the transition from stationary to flowing is, in this case, no longer left-right and clock-counterclockwise symmetric, as in the examples discussed in Section IV.A, for any τ ≠ 0. This is well illustrated by the bifurcation diagram of Fig. 4b, showing the departure in the director orientation from the boundary value at the center of the channel [i.e. θ(1/2) − θ0, with L = 1]. The dots have been obtained from a numerical integration of eqn (18), whereas the solid lines correspond to analytical solutions obtained by solving a third order equation for the constant C in eqn (25), as in Section IV.A.1. We find:

 
image file: c9sm01851d-t35.tif(37)
where image file: c9sm01851d-t36.tif and image file: c9sm01851d-t37.tif. For α > αc, the solution consists of two branches, of which only one is connected with the stationary solution θ(y) = θ0. Furthermore, the gap between the two branches increases monotonically with τ. If the instability is triggered upon applying a small random perturbation to the stationary state, this will always select the closest branch, thus the one connected to the trivial solution. As a consequence, a chiral cellular monolayer driven out of the stationary state by a small perturbation, will systematically tilt and flow in the same direction, which is in turn determined by the sign of the chiral active stress τ. The transition described above is known in bifurcation theory as a perturbed or imperfect pitchfork bifurcation and occurs when a standard pitchfork bifurcation, whose normal form is θ3μθ = 0, is biased by a small symmetry-breaking perturbation: i.e. θ3μθ + PL + PSθ2 = 0, where μ, PL, and PS are constant parameters. If PL = PS = 0, the equation is invariant under θ → −θ. Thus, for μ > 0, the trivial solution is unstable and the transition is supercritical, while for μ < 0, only the trivial solution is stable, and the bifurcation is subcritical. By contrast, for non-vanishing PL and PS, the equation is no longer invariant under θ → −θ and the bifurcation is no longer symmetric (see ref. 41 and 42 for an overview).

In our case, the role of the symmetry-breaking perturbation is played by the chiral stress τ. Thus, in the unperturbed scenario, τ = 0 and the stationary solution is θ = 0, with the critical α value being αc = 2π2ηK/[γL2(1 − λ)]. When τ ≠ 0, on the other hand, expanding eqn (21) around θ = 0 and using ∂y2θ = −q2θ one finds:

 
image file: c9sm01851d-t38.tif(38)
Evidently, this coincides with the normal form of a perturbed bifurcation for any finite τ value. For τ = 0, on the other hand, one recovers the normal form of the symmetric pitchfork bifurcation. Upon increasing τ, the bifurcation is shifted toward smaller α values, until, for τ = π2/ηK(γL2), αc = 0 and the system is never stationary.

V. Discussion and conclusions

In this article we have investigated how a chiral and yet symmetric stress tensor might arise microscopically in nematic cell monolayers and how such a chiral stress influences some of the hallmark phenomena of active nematics. In Section II we proposed a microscopic model where a misalignment of the active force dipole and the cell's orientation is demonstrated to lead to a macroscopic chiral active stress tensor of the form eqn (2). In Section III we showed how the presence of chiral active stresses tilts the flow around ±1/2 disclinations, thereby leading to a misalignment between the defect polarity and the direction of motion, by an angle θtilt = arctan(2τ/α), with τ and α the chiral and non-chiral active stress respectively. In Section IV we investigated the spontaneous flow transition in a quasi-one-dimensional channel for both no-slip and stress-free boundary conditions as well as for various types of anchoring. For no-slip boundaries (Section IV.A), we recovered the classic pitchfork bifurcation first discussed by Voituriez et al.16 In this case the chirality does not affect the transition itself, but does leave a signature on the post-transitional configurations of the nematic director and velocity field, in the form of asymmetry with respect to the channel centerline. In case of stress-free boundaries (Section IV.B), we found that chirality renders the stationary and uniformly aligned configuration incompatible with most of the anchoring conditions. As a consequence, the cellular monolayer is always in motion, for any non-vanishing chiral and non-chiral active stress. For very large active stresses, in particular, we found an active analog of flow-alignment in nematics, with the bulk orientation (analogous to the Leslie angle40) set by the ratio between chiral and non-chiral active stresses, i.e.θtilt/2. Finally, in the special case in which the nematic director is anchored at an angle −θtilt/2 at the channel walls, we found that the spontaneous flow transition becomes asymmetric, i.e. only one of the two branches of the pitchfork bifurcation is connected to the trivial solution, which may then be the only one observed experimentally. This latter result could potentially explain the experimental observations by Duclos et al.,8 who found that NIH 3T3 cells are more likely to tilt clockwise then counterclockwise once the spontaneous flow transition sets up. In ref. 43 and 44 it was found that at the interface between an active nematic phase, with negligible distortions of the director field, and an isotropic phase there is an active anchoring angle that is set by the activity. This discussion can easily be repeated for the chiral active stress tensor in eqn (2) with the result being that the angle is changed due to the presence of chirality. Furthermore, considering a confined nematic phase without isotropic phase but with distorted director field, the case considered in Section IV being one example, it is readily found that the active anchoring angle is equal to the angle θtilt introduced above. Thus, the active anchoring angle considered in ref. 43 and 44 (no distortions of director field but region with active nematic and isotropic phase) and the anchoring angle θtilt derived above (distorted director field but only a nematic phase) can be seen as being the two special cases of the more general case with a nematic–isotropic phase interface and non-negligible distortions of the director field.

Further experimental investigations into the influence of chirality would be interesting. In particular, the tilt of the flow around ±1/2 disclinations has, according to our knowledge, not yet been observed. Thus, measurements of the tilt angle and experimental investigations of the flow field are needed to compare the theory with real-life cell monolayers. Additionally, as mentioned, the tilt angle opens a possibility to determine the relative magnitude of the chiral stress directly by particle-image-velocimetry measurements. Furthermore, since the cells used in ref. 8 were only weakly chiral the effects of chirality were not as pronounced. Performing similar experiments with cells with stronger chirality and for different boundary conditions would enable further tests of the presented theory.

Conflicts of interest

There are no conflicts to declare.

Appendix A: stress owing to a force quadrupole

The derivation of the active stresses given in Section II can be straightforwardly generalized to account for a more complex force distribution. For illustrative purposes we consider here the case of a quadrupole consisting of two force dipoles applied at the ends of the cell in longitudinal and transverse directions. In this case the force density field fc is given by:
 
fc = F(a)cδ(Rcaνc) − F(a)cδ(Rc + aνc) + F(b)cδ(Rcbνc) − F(b)cδ(Rc + bνc),(A1)
where a and b are, respectively, the major and minor semiaxis. Taking F(i)c = F(i)n + F(i)n, with i = a, b, and coarse-graining over the scale of a volume element, we find: Pa = −(aF(a) + bF(b)), α = 2(aF(a)bF(b)) and τ = aF(a) + bF(b).

Appendix B: nonlinear expansion

Expanding eqn (21) about θ0 = 0 up to third order in δθ yields
image file: c9sm01851d-t39.tif
which can be written as
 
image file: c9sm01851d-t40.tif(B1)

The stress σxy is determined by the no-slip boundary condition and force balance, i.e.,

 
image file: c9sm01851d-t41.tif(B2)
An expansion up to third order in δθ = C[thin space (1/6-em)]sin(2πy/L) around θ = 0 yields the condition
 
image file: c9sm01851d-t42.tif(B3)
and thus σxy = τ − 2τC2[thin space (1/6-em)]sin2(2πy/L) = τ(1 − 2δθ2). This determines σxy and can be used in eqn (B2) to find
 
image file: c9sm01851d-t43.tif(B4)
Comparing with the normal form given in Section IV.B.3 we find that, for L > Lc, μ = 3(λ − 1)/(2(4λ − 1)) is positive for all λ < 1/4 and λ > 1 and vanishes for λ = 1. With δθ = C[thin space (1/6-em)]sin[thin space (1/6-em)]y/L solving for C at y = L/2 and expanding about LLc to leading order yields the trivial solution C = 0 and
 
image file: c9sm01851d-t44.tif(B5)
which has a singularity at λ = 1/4. As remarked on in the general discussion in Section IV.B.3 the system displays a supercritical bifurcation for μ > 0 but a subcritical bifurcation for μ < 0. Thus the character of the bifurcation is determined by the sign of μ and in the present case this sign changes at the singularity λ = 1/4 and at λ = 1, where the system transitions between flow-tumbling (λ < 1) and flow-aligning (λ > 1).

Acknowledgements

We are grateful to Mattia Serra for illuminating discussions about Section IV. This work was supported by funds from the Netherlands Organisation for Scientific Research (NWO/OCW), as part of the Frontiers of Nanoscience program (L. G.), the Vidi (L. G. and L. A. H.), and Vici (R. M. H. M.) schemes and the Leiden/Huygens fellowship (K. S.).

References

  1. G. K. Batchelor, J. Fluid Mech., 1970, 41, 545 CrossRef.
  2. G. Duclos, S. Garcia, H. G. Yevick and P. Silberzan, Soft Matter, 2014, 10, 2346 RSC.
  3. S. Garcia, E. Hannezo, J. Elgeti, J. F. Joanny, P. Silberzan and N. S. Gov, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 15314 CrossRef CAS PubMed.
  4. G. Duclos, C. Erlenkämper, J. F. Joanny and P. Silberzan, Nat. Phys., 2016, 13, 58 Search PubMed.
  5. T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212 CrossRef CAS PubMed.
  6. K. Kawaguchi, R. Kageyama and M. Sano, Nature, 2017, 545, 327 CrossRef CAS PubMed.
  7. C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos, L. Giomi and P. Silberzan, Phys. Rev. Lett., 2018, 120, 208101 CrossRef CAS PubMed.
  8. G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J.-F. Joanny, J. Prost and P. Silberzan, Nat. Phys., 2018, 14, 7 Search PubMed.
  9. A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 3246 Search PubMed.
  10. R. Kemkemer, D. Kling, D. Kaufmann and H. Gruler, Eur. Phys. J. E: Soft Matter Biol. Phys., 2000, 1, 215 Search PubMed.
  11. R. Kemkemer, V. Teichgräber, S. Schrank-Kaufmann, D. Kaufmann and H. Gruler, Eur. Phys. J. E: Soft Matter Biol. Phys., 2000, 3, 101 Search PubMed.
  12. L. Giomi, M. J. Bowick, X. Ma and M. C. Marchetti, Phys. Rev. Lett., 2013, 110, 228101 CrossRef PubMed.
  13. L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. C. Marchetti, Philos. Trans. R. Soc., A, 2014, 372, 20130365 CrossRef PubMed.
  14. L. Giomi, Phys. Rev. X, 2015, 5, 031003 Search PubMed.
  15. S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek and E. Bertin, 2019, arXiv:1901.04763.
  16. R. Voituriez, J.-F. Joanny and J. Prost, Europhys. Lett., 2005, 70, 404 CrossRef CAS.
  17. S. Fürthauer, M. Strempel, S. W. Grill and F. Jülicher, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 89 CrossRef.
  18. S. Fürthauer, M. Strempel, S. W. Grill and F. Jülicher, Phys. Rev. Lett., 2012, 110, 048103 CrossRef.
  19. D. Banerjee, A. Souslov, A. G. Abanov and V. Vitelli, Nat. Commun., 2017, 8, 1573 CrossRef.
  20. J. E. Avron, R. Seiler and P. G. Zograf, Phys. Rev. Lett., 1995, 75, 697 CrossRef CAS.
  21. T. J. Pedley and J. O. Kessler, Annu. Rev. Fluid Mech., 1992, 24, 313 CrossRef.
  22. R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101 CrossRef PubMed.
  23. A. W. C. Lau and T. C. Lubensky, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 011917 CrossRef CAS PubMed.
  24. L. Q. Wan, K. Ronaldson, M. Park, G. Taylor, Y. Zhang, J. M. Gimble and G. Vunjak-Novakovic, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 12295 CrossRef CAS PubMed.
  25. Y. H. Tee, V. Thiagarajan, R. F. Hariadi, K. L. Anderson, C. Page, N. Volkmann, D. Hanein, S. Sivaramakrishnan, M. M. Kozlov and A. D. Bershadsky, Nat. Cell Biol., 2015, 17, 445 CrossRef CAS PubMed.
  26. K. Schakenraad, J. Ernst, W. Pomp, E. H. J. Danen, R. M. H. Merks, T. Schmidt and L. Giomi, 2019, arXiv:1905.09805.
  27. K. E. Worley, D. Shieh and L. Q. Wan, Integr. Biol., 2015, 7, 580 CrossRef CAS PubMed.
  28. T. Sanchez, D. N. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 431 CrossRef CAS.
  29. F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic and A. R. Bausch, Science, 2014, 345, 1135 CrossRef CAS.
  30. P. Guillamat, J. Ignés-Mullol and F. Sagués, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 5498 CrossRef CAS.
  31. P. Guillamat, J. Ignés-Mullol and F. Sagués, Nat. Commun., 2017, 8, 564 CrossRef CAS.
  32. R. Zhang, N. Kumar, J. L. Ross, M. L. Gardel and J. J. de Pablo, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, E124 CrossRef CAS.
  33. A. J. Vromans and L. Giomi, Soft Matter, 2016, 12, 6490 RSC.
  34. A. Maitra and M. Lenz, Nat. Commun., 2019, 10, 920 CrossRef.
  35. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 031921 CrossRef CAS.
  36. L. Giomi, M. C. Marchetti and T. B. Liverpool, Phys. Rev. Lett., 2008, 101, 198101 CrossRef.
  37. S. A. Edwards and J. M. Yeomans, Europhys. Lett., 2009, 85, 18008 CrossRef.
  38. J. Hardoüin, R. Hughes, A. Doostmohammadi, J. Laurent, T. Lopez-Leon, J. M. Yeomans, J. Ignés-Mullol and F. Sagués, 2019, arXiv:1903.01787.
  39. V. Fréedericksz and A. Repiewa, Z. Phys., 1927, 42, 532 CrossRef.
  40. P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, Oxford University Press, Oxford, 1993 Search PubMed.
  41. M. Golubitsky and D. G. Schaeffer, Singularities and Groups in Bifurcation Theory, Springer-Verlag, New York, 1985, vol. I Search PubMed.
  42. P. M. Kitanov, W. F. Langford and A. R. Willms, Dynam. Cont. Dis. Ser. A, 2013, 20, 197 Search PubMed.
  43. M. L. Blow, S. P. Thampi and J. M. Yeomans, Phys. Rev. Lett., 2014, 113, 248303 CrossRef.
  44. A. Doostmohammadi, S. P. Thampi and J. M. Yeomans, Phys. Rev. Lett., 2016, 117, 048102 CrossRef.

This journal is © The Royal Society of Chemistry 2020