Viscoelastic confinement induces periodic flow reversals in active nematics

We use linear stability analysis and hybrid lattice Boltzmann simulations to study the dynamical behaviour of an active nematic confined in a channel made of viscoelastic material. We find that the quiescent, ordered active nematic is unstable above a critical activity. The transition is to a steady flow state for high elasticity of the channel surroundings. However, below a threshold elastic modulus, the system produces spontaneous oscillations with periodic flow reversals. We provide a phase diagram that highlights the region where time-periodic oscillations are observed and explain how they are produced by the interplay of activity and viscoelasticity. Our results suggest new experiments to study the role of viscoelastic confinement in the spatio-temporal organization and control of active matter.

We use linear stability analysis and hybrid lattice Boltzmann simulations to study the dynamical behaviour of an active nematic confined in a channel made of viscoelastic material.We find that the quiescent, ordered active nematic is unstable above a critical activity.The transition is to a steady flow state for high elasticity of the channel surroundings.However, below a threshold elastic modulus, the system produces spontaneous oscillations with periodic flow reversals.We provide a phase diagram that highlights the region where time-periodic oscillations are observed and explain how they are produced by the interplay of activity and viscoelasticity.Our results suggest new experiments to study the role of viscoelastic confinement in the spatio-temporal organization and control of active matter.
Living systems across scales exhibit collective motion, and thus spatio-temporal patterns, vividly manifested as, for instance, motility-induced phase separation [1], spontaneous flow transitions [2][3][4][5], and turbulence at low Reynolds number [6][7][8].Not only biochemical and genetic cues but mechanical interactions of the system with its surroundings are important in dictating such emergent dynamics.Adding to this complexity, biological environments are often endowed with viscoelastic properties, for example, biofilms where bacterial cells colonize in a polymeric matrix [9], migration of cells through extracellular matrix [10][11][12][13], notably the phenomenon of durotaxis [14], and change in the swimming behaviour of microorganisms due to the presence of polymers in biofluids [15,16].In a different context, traction force microscopy has become an indispensable tool to probe force fields in cellular structures.The technique assumes a one-way mechanical interaction of cells with an elastic substrate [17,18].Therefore, clarifying the interplay of the viscoelasticity of a confining medium and activity of the living system is crucial from understanding measurements in mechanobiology to biological events such as wound healing [19], morphogenesis [20], and cancer invasion [21].Besides, identifying universal pathways of pattern formation is a central goal of active matter research.
It is well known that active nematics, a versatile model fluid for active matter, confined in a rigid channel displays a transition-mathematically analogous to the Fredericks transition in passive liquid crystals [22]-from quiescence to a flow state when the activity is increased beyond a threshold value [2,[23][24][25][26].Further increase in activity induces a cascade of dynamical transitions resulting in oscillatory flows [4,27], dancing topological defects [27,28], and active turbulence [29][30][31][32][33]. Thus, channelconfined active nematics have become a paradigm for understanding the dynamical behaviour active systems [34].Therefore we investigate the interaction between activity and viscoelasticity by analyzing an active nematic flowing in a soft channel.
Previous studies that address the role of viscoelasticity in living systems considered either active particles within a viscoelastic fluid [35][36][37][38][39][40] or active matter in contact with a viscoelastic environment [41][42][43].In the former case, oscillating vortices and drag reduction effects are seen to arise due to the presence of polymers [38,39].
In the latter, less studied case, numerical simulations demonstrate that temporal pulses in activity drive reversal of spontaneous flows [42].In this letter, we demonstrate analytically and numerically that, above a critical activity, viscoelastic confinement produces spontaneous, oscillatory flow states of an active nematic that switches flow directions periodically.The direction-reversing oscillatory flows exist only in 'soft' channels, and they disappear when the elastic modulus of the confinement increases above a critical value.Building on our findings, we explain the origin of oscillations as the interplay of activity and viscoelasticity, demonstrate the generality of the phenomenon and discuss the consequences.We consider a two-dimensional channel of width 2L and infinite length which contains the active nematic.The borders of the channel which span a width of (β − 1)L on either side are made up of viscoelastic material (see Fig. 1).Let x and y denote the directions parallel and perpendicular to the channel length, with y = 0 the centerline of the channel.The relevant hydrodynamic variables are Q and v representing the orientational order and velocity field in the active nematic respectively and u the displacement field in the viscoelastic layers.
Active nematics may develop orientational order either due to the elongated shape of the constituents, [44,45] or as an emergent feature of deformability of particles, such as cells [46] or due to activity itself [47].The nematic order is measured using an orientational order parameter Q = 2q(nn − I/2), where n = (cos(θ), sin(θ)) is the director field, θ ∈ (−π/2, π/2) is the angle that the nematogens form with the positive-x direction, q is the magnitude of the nematic order and I is the identity tensor.The nematic order parameter tensor evolves according to [48] where is the strain rate tensor, and Ω = ((∇v) − (∇v))/2 is the vorticity tensor.The flow aligning parameter λ is determined by the shape of the nematogens.In Eq. ( 1), γ is the rotational viscosity and H = −δF/δQ is the molecular field which drives the system to the minimum of the free energy with energy density F = Here, K is the elastic constant, and A and C are material parameters, chosen so that the system is in the nematic phase at equilibrium.
The velocity field v obeys the incompressible Navier-Stokes equations [6,49]: where the total stress tensor σ is given by the sum of (i) the viscous stress σ viscous = 2η 1 E, where η 1 is the viscosity of the active nematic, (ii) the elastic stress , where P 1 is the bulk pressure, and (iii) the active stress σ active = −ζQ.
The dynamics of the incompressible viscoelastic layers is described by the displacement field u from the equilibrium position, that evolves according to [50,51] where ρ 2 is the gel density and P 2 is the bulk pressure in the viscoelastic layers.The stress tensor τ is model dependent and we consider two simple yet powerful constitutive relations, namely the to capture the rheological response of the viscoelastic layers that confine the active nematic.In the above, D/Dt is the upper convected derivative [51], E and η 2 are the elastic modulus and viscosity respectively.A Maxwell (Kelvin-Voigt) material is composed of a spring and a dashpot connected in series (parallel).It behaves as an elastic solid at short (long) times and as a viscous liquid at long (short) times, with a single crossover timescale η 2 /E.Eqs. ( 1)-( 3) govern the dynamics of the system and we solve them (i) analytically as a linear stability problem and (ii) numerically using a hybrid lattice Boltzmann method [52].We assume translational invariance in the x-direction, so that v y = 0; v x = v x (y) and u y = 0; u x = u x (y).The viscoelastic material is in contact with a no-slip wall at y = ±βL.At the interface between the active nematics and the viscoelastic layer, we impose noslip conditions, v x (±L) = ∂ t u x (±L), and continuity of the stress tensor σ xy (±L) = τ xy (±L).For simplicity, we consider strong planar anchoring of the director field at the interface, i.e., θ(±L) = 0.
To investigate the interplay of activity and viscoelasticity, we perform linear analysis to calculate the stability of a small perturbation around the static nematic state with (v x , u x , θ, q) = (0, 0, 0, q 0 ) where q 0 = −A/(2C).
For simplicity, we first consider a purely elastic material bounding the nematic fluid, corresponding to the limit η 2 → ∞ (η 2 → 0) for the Maxwell (Kelvin-Voigt) model.In the limit of large elastic modulus E → ∞, the boundaries at y = ±L are rigid and we recover the classical result of Voituriez et al. [2]: increasing the activity beyond a critical value ζ c wall , the nematically ordered state is unstable and spontaneous flows develop driven by the distortions in the director field.The critical activity is calculated from Eq. ( 4), In the opposite limit E = 0, corresponding to a freestanding film of active nematic, an analogous transition to a steady flow is observed at activity ζ c free = ζ c wall /4.The critical activity of the system at intermediate values of E, obtained from Eq. ( 4), is summarized in Fig. 2 (red line).The critical activity ζ c = ζ c free at E = 0, and increases with increase in the elastic modulus E, until a threshold elastic modulus E = E c .Beyond E c the critical activity "freezes" to ζ c = ζ c wall , that corresponding to a rigid wall.E c can be determined analytically from Eq. (4), see [52].
Interestingly the transition mechanism at ζ c , at which the ordered nematic state becomes unstable, is different for E < E c and E > E c .We find that, for E < E c the route to instability is via a Hopf bifurcation where the complex conjugate eigenvalues ω cross the imaginary axis with a finite imaginary part at ζ = ζ c (Fig. 2b).Consequently, the ensuing instability is oscillatory and the active nematic transitions from a quiescent to an oscillating state where the flow direction is reversed periodically.On the other hand, for E > E c , the instability becomes stationary (Im(ω) = 0) and no oscillations are observed.The numerical simulations show that the oscillations are replaced by steady flow at sufficiently high activity (see Fig. 2).
The oscillatory state can be understood by following the temporal evolution of a system which is at its critical activity ζ = ζ c , and with 0 < E < E c and |ζ c free | < |ζ| < |ζ c wall | (such as a point marked ' ' in Fig. 2(a)).At time t = 0 (see Fig. 3a), the elastic layer is not deformed (u x = 0), and the stress at the active nematic-elastic interface (y = ±L) vanishes.This condition corresponds to a free standing active nematic film (no resistance from the elastic layer), which will have a critical activity ζ c free .Since the activity of the system exceeds this critical value, |ζ| > |ζ c free |, spontaneous flow develops in the active film.The velocity profile ṽx (y) is an odd function of y [52] similar to that of a shear flow.These flows, in turn, drive the deformation of the elastic confinement.Eventually, the elastic response of the channel wall slows down the flow and the deformation rate at the active-elastic interface vanishes.In this configuration, the effect of elastic confinement is the same as that of a rigid wall and the critical activity for the active nematic is ζ c wall .However, since |ζ| < |ζ c wall |, the active forcing is not sufficient to sustain the flows and they die out (Fig. 3b).The elastic energy stored in the elastic medium pushes the flow in the opposite direction, leading to a flow reversal (Fig. 3c).Hence, the oscillations arise because the activity is too high to remain in the quiescent state (ζ > ζ c free ) but too low to sustain the flow (ζ < ζ c wall ).The period T of the oscillations is set by the elasticity, viscosity and L/(γ −1 K), the relaxation timescale of the director field.The period T close to the critical point ζ = ζ c can be obtained analytically [52] from Eq. ( 4) and is shown in Fig. 2c as a function of E. For E → 0, the activity ζ is only slightly larger than ζ c free required to initially start a flow, leading to a slowdown of the dynamics.Similarly, when E → E c , the activity ζ is only marginally below ζ c wall and the flow-reversal mechanism again slows down significantly.Indeed the time period diverges in the limiting cases: Hence, the crossover from oscillatory to steady flow at the two limiting cases, E > 0 to E = 0 and E < E c to E = E c occurs smoothly via an infinite-period bifurcation.The period T has a minimum at E = E * , reminiscent of the phenomenon of resonance and the elastic modulus can be optimally tuned to increase the frequency of oscillatory motion.
To gain further insight into the oscillatory modes of the instability, we next plot the trajectory of the system in a phase space spanned by the displacement of the elastic layer (u x (y = L, t)) and the velocity of the active nematic (v x (y = L/4, t)) as shown in Fig. 3e.The exact shape of the curve depends on the choice of parameters, but note that the phase space trajectory encloses a finite area indicating the phase lag in the the velocity field of active nematic and the displacement field of elastic confinement.Interestingly, the phase space trajectory does not coincide with the time-reversed trajectory (u x (L, −t), −v x (L/4, −t)), manifestly breaking the timereversal symmetry and showing the non-equilibrium nature of the active-dissipative system under consideration.While non-reciprocal oscillatory motion, the sine qua non for self-propulsion (the scallop theorem), is abundant in life at low Reynolds number [53,54], the current analysis demonstrates that the mechanical coupling of activity and elasticity automatically generates such nonreciprocal motion in active systems.
Having established that the genesis of oscillations is the elasticity of the confining channel we can analyze more complex constitutive relations.For the Maxwell model, on time scales smaller than η 2 /E the viscoelastic confinement behaves as an elastic solid and the coupling between activity and elasticity still leads to oscillations as illustrated in Fig. 4 .The behaviour at small E can be understood in a similar fashion.In this limit, the viscoelastic timescale η 2 /E is large compared to the period of the oscillations T ∼ 1/ √ E, and the Maxwell material behaves as an elastic solid exhibiting an η r independent behaviour of ζ c .In particular, ζ c ≈ ζ c free at small elasticity E.
Opposite behaviours are observed when the channel confinement is the Kelvin-Voigt material.For E > E c the instability is still stationary but the confining material now behaves as an elastic solid and wall .This results in the threshold elastic modulus E c being independent of the viscosity η 2 as illustrated in Fig. 4(b).On the other hand, for E → 0, the viscoelastic timescale η 2 /E is large compared to the period of the oscillations and the Kelvin-Voigt material behaves as a viscous fluid.Hence, the critical activity ζ c for small E strongly depends on the viscosity ratio η r .To summarize, choice of different constitutive models of the channel confinement leads to quantitative differences but does not change the physics of the oscillations.
Our results highlight a novel pathway to spatiotemporal pattern formation in active matter.It is indeed remarkable to note that the time periodic, oscillatory flows arise even at constant activity.Our predictions can be tested experimentally, by confining cell layers [5] or microtubule-based active fluids [25] in channels with soft walls.Moreover, traction force microscopy provides a potential platform to study the role of an active-elastic boundary [17,18].In addition to extracting work from active materials, the nonreciprocal dynamics that arise from the interplay of activity and viscoelasticity might also be utilised to make self propelling objects.
Solving the gel equation ( 12) and imposing the no-slip boundary condition ũx (±βL) = 0 we find for y > L (y < −L).From Eq. ( 11), we obtain where we have defined We first consider the even solution where we have imposed the boundary condition θ(±L) = 0 and defined Using Eq. ( 10), we find Imposing the boundary conditions at y = ±L, we find the following condition for ω γ Considering the odd solution, we find and the condition For the range of parameters considered in the paper, we find that the even solution (corresponding to no net flow in the channel) is dominant, i.e., it becomes unstable at lower values of the activity.Hence, in the main text, we only focus on the even mode.The odd solutions may be favored by introducing weak anchoring.

Asymptotic behaviors
In this section, we extract the asymptotic behavior of the solution of Eq. ( 23).For simplicity, we set q 0 = 1 and we consider the case ρ 1 = ρ 2 = 0.In this limit, the condition in Eq. ( 23) becomes where We first consider the limit of small E. For E = 0 (corresponding to a free surface), the critical value of the activity can be computed analytically and reads We set E = , ζ = ζ c free + a 1 , and ω = a 2 √ .We then expand Eq. ( 27) in powers of , yielding Setting the coefficients to zero, we get and As expected, the growth rate ω is purely imaginary.
To investigate the asymptotic behavior of the system close to the transition, we set Eq. ( 36) is transcendental and must be solved numerically to determine ζ c visc .Considering higher order expansions, one can find expressions for a 1 and a 2 .In the limit of a purely elastic medium (η 2 → ∞) we find Kelvin-Voigt model In this section, we perform the linear stability analysis in the case of the Kelvin-Voigt model.The constitutive relation reads Hence, the gel displacement u x (y, t) evolves according to Assuming u x (y, t) = e ωt ũx (y), we find where we have defined Following the same derivation as for the Maxwell model, we find two instabilities, corresponding to the even and odd solutions for θ.The condition for the growth rate of the even solution reads For the odd solution, we find As for the Maxwell model, for the range of parameters considered in the paper, we find that the even solution is dominant.

FIG. 1 .
FIG. 1. Schematic representation of the system: an active nematic layer of width 2L is confined between two viscoelastic layers, each of width (β −1)L.Thus the bounding rigid plates are separated by a distance 2βL.The active nematic is a dense suspension of elements that generate active stress.The viscoelastic layers are shown as made up of Maxwell elements.

FIG. 2 .
FIG. 2. a):Phase diagram in the (ζ, E) plane, illustrating the states of an active nematic when confined in a soft channel.The continuous red line is the critical activity ζ c , obtained from linear stability analysis, at which the nematic state becomes unstable, driving flows.The symbols are obtained from hybrid lattice Boltzmann simulations with η1 = 10/3, η2 = ∞ (elastic limit), γ = 10, ρ1 = 20, ρ2 = 0, K = 0.1, q0 = 0.25, L = 10, λ = 0, β = 2.For E < Ec ≈ 0.00326, the instability leads to periodic oscillations.For E > Ec, a steady flow of active nematic is obtained.b): The growth rate ω in the complex plane, for E = 0.002 < Ec (blue lines) and E = 0.0042 > Ec (orange line).The arrows indicate the direction of increasing |ζ|.c): Time period of oscillations T as a function of E/Ec.T diverges when E → 0 and E → Ec.

FIG. 3 .
FIG. 3. Temporal evolution of the states of the system with the same parameter values as in Fig. 2 at ζ = ζ c (E).The panels a), b), c), and d) show the hydrodynamic fields in the oscillatory phase for t = 0, t = T /4, t = T /2, and t = 3T /4.Panel e) shows the phase space trajectory in the (ux(L, t), vx(L/4, t)) plane for different values of E. The dashed lines display the time-reversed trajectories (ux(L, T − t), −vx(L/4, T − t)), showing time-irreversibility.