Stability analysis of flow of active extensile fibers in confined domains

The study of active matter in confinement is very relevant to understanding many biological systems as well as in guiding potential applications. In contrast with passive systems, whose dynamics is the direct result of external agents (e.g., a pump in the Poiseuille flow), active matter systems are able to use stored energy or extract it from external sources at small (or at individual) length scales and collectively enhance it at the macroscopic scale, often in cooperation with the environment (e.g., bacteria swimming in liquid crystal media), and convert it into work.^1,2 Typical examples of these systems also include cytoplasm networks and cancer cells.

The aim of our work is to provide theoretical support to recent experiments of shear flow of extensile fibers in confined channels of varying width in order to understand complex phenomenology of active flow in a controlled setting.³ Earlier theoretical work on extensile fibers found that planar confined active nematics that are aligned laterally undergo an instability, leading to a spontaneous laminar flow when the channel width reaches a certain threshold value that depends on the strength of the activity. In an analogous work, Voituriez et al.⁴ analyzed the stability of an aligned state of an active polar gel. Their results were later confirmed in experiments with spindle-shaped cells.⁵ 

The active system studied in earlier work³ that motivates our analysis is composed of microtubules concentrated in bundles, tens of microns long, by the addition of the depleting agent polyethylene glycol (PEG). Kinesin motors within each bundle bridge neighboring microtubules and walk toward their positive ends causing internal sliding of filaments of opposite polarity. The aqueous suspension of active filaments is placed in a custom-made pool 5 mm in diameter and covered with a layer of silicon oil. About 30 min later, the aqueous/oil interface is coated with a layer of the active filaments, self-organized as a two-dimensional active nematic liquid crystal. The confinement of the system is achieved by placing micro-printed polymer ridges of $L1:=1.5$ mm in length, $100μ$m in depth, with a width range between $L3:=30$–$300$ $μ$m. We point out the small aspect ratio of the flow confining domain. Because it seems clear from the experiment that the long dimension along the channel does not play any role, our analysis below considers an invariant system of fixed length $0<L≤L1$ within the flow domain.

The main control parameters of the experiment are the channel width and the concentration of ATP, which determines the activity of the system. The experiments identify two distinguished regimes: (i) a defect-free shear flow, observed in channels of width $L3<80μ$m and (ii) a dancing disclination regime for $L3>90μ$m, with pairs of defects of fractional degree $±12$ that form and annihilate on the boundary. The events of defect annihilation drive a flow in the channel, as reported in an earlier article.⁶ Our work focuses on the small channel width, defect-free, shear flow aligning regime.

The Ericksen–Leslie equations of incompressible liquid crystal flow and their numerical simulations have been used in many studies of active systems, especially in the turbulent regimes.⁷ The variables of the model include the velocity field $v$⁠, the pressure $p$ associated with the incompressibility constraint, and the unit nematic director $n$⁠, representing the local average of angular alignment of the fibers in the flow.

The dimensionless form of the equations involves, in addition to the aspect ratio $ℓ$ of the channel, three main parameter groups. Two of them, the Reynolds $Re$ and Ericksen $Er$ numbers, are defined for passive flows. The activity parameter $A$ quantifies the activity level in the system. Although $Re$ is significantly smaller than other nondimensional parameters, it influences the stability behavior of the active system in subtle ways. The $Er$ represents the ratio of the viscous to the elastic torques and may take very large values.

In this work, we focus on the plane shear flow of an active liquid crystal in the flow aligning regime. The domain of the flow is taken as the patterned channel domain of the experiment. We find solutions with linear velocity profiles and constant angles of alignment that agree with experimental observations of active filaments in a channel when $A>0$⁠.³ The ratio of the extensional effect of a shear flow vs the rotating one is represented by the positive shear aligning parameter $λ$⁠. Aligning regimes, characterized by the dominance of the extensional effect, correspond to $λ>1$ (for tumbling regimes $λ<1$⁠). Indeed, the dominance of the extensional effect causes the liquid crystal to align at a certain angle with respect to the direction of the velocity. We find that the magnitude of the velocity gradient is proportional to the activity parameter, and the angle of alignment solely depends on the ratio of two relevant anisotropic viscosity coefficients, as in the case of passive shear flow. The velocity profile quantitatively matches the profile found in the experiments.³ 

Several works on active liquid crystals found in the literature use the Beris–Edwards model based on the evolution of the order tensor $Q$ rather than the director of the Ericksen–Leslie theory.⁸ $Q$ represents a symmetric, traceless second order tensor whose two independent eigenvectors are the director fields of the theory. They reduce to a single eigenvector $n$ in the uniaxial case, with $Q$ admitting the representation $Q=s(n⊗n−13I),$ where the scalar $s∈(−12,1)$ corresponds to the single, independent eigenvalue. Our work would correspond to a constant value $s=s0>0,$ the positiveness agreeing with dealing with rod-like liquid crystals. A recent article⁹ presents the effects of a pressure gradient on the spontaneous flow of an active nematic liquid crystal in a channel, subject to planar anchoring and no-slip conditions on the boundaries of the channel. The work is also based on the Ericksen–Leslie theory with active stress, studies the director angle at the asymptotic limit of high activity, and numerically solves the nonlinear equations to determine the effect of the pressure on the equilibrium solutions. A great bulk of the work on active systems involves numerical studies of the high activity, chaotic regime, with the presence of defects.^6,7,10–15

We carry out a normal mode analysis of the shear flow, with $ω$ representing the spatial frequency of the perturbation and $Im(c)$ the corresponding growth rate. We find that a perturbation with frequency $ω$ is unstable for $Im(c)>0$⁠. The threshold $Im(c)=0$ represents neutral stability. Since the perturbing flows that we study are planar, we formulate the corresponding stability problem in terms of the stream function $Ψ$⁠. We choose boundary conditions so that the boundary values of the base flow are not altered. The zero boundary conditions on perturbations suggest that a discretization using Chebyshev polynomials is appropriate. We apply the Chebyshev–QZ-algorithm to solve the generalized eigenvalue problem resulting from the linearization of the governing system about the base shear flow. We emphasize the secondary role of the director field boundary conditions. Indeed, one important difference between active and passive liquid crystals is that whereas the latter can be aligned by surface anchoring, this is not the case for active flows, known to respond to alignment by flow only.

The stability plots in the $Aω$-plane indicate threshold values of $A>0$ at which the system becomes unstable. The critical value decreases with increasing values of $Er$ and $ℓ$ and with the range of unstable frequencies also increasing. The growth rate profile calculated with respect to $A$ shows a parabolic shape at low frequencies, becoming linear at a threshold frequency value. The profile of the growth rate $Im(c)$ with respect to $ω$ starts being positive for all values of $A$⁠, becoming $0$ at a threshold value of $ω$ that appears to be independent of $A$⁠. With further increase of $ω$⁠, $Im(c)$ reaches a positive maximum, subsequently tending to zero as $ω$ increases, a feature that is also independent of $A$⁠. The instability profiles obtained from the analysis qualitatively follow the experimentally reported trends,³ in that the range of unstable frequencies increases with activity and channel width. Our work also provides quantitative data on the latter, but it cannot be directly verified from the available experimental information.

This paper is organized as follows. Section II presents the formulation of the Ericksen–Leslie system for active liquid crystals, including the assumptions that lead to aligning flows and a description of the boundary conditions. The scaling of the problem and non-dimensionalization of the equations are presented in Sec. III. The shear flow geometry and calculation of the corresponding steady states are presented in Sec. IV. The framework for the stability analysis is developed in Sec. V. The numerical method to analyze the stability of shear flows is presented in Sec. VI. In Sec. VII, we present the stability results. In Sec. VIII, we present the conclusions.

The fully developed equations for the Ericksen–Leslie and Ericksen models are presented in the supplementary material. A summary of the plots that resulted from the numerical simulations is also given there.

Similar to its passive counterpart, an active liquid crystal is assumed to be a viscous anisotropic and incompressible fluid with activity sources drawn from internal mechanisms or from the environment. Let $Ω⊂R3$ be an open domain occupied by the liquid crystal with the smooth boundary $∂Ω$⁠. The Ericksen–Leslie equations of balance of linear and angular momentum and the incompressibility and unit director constraints for the velocity field $v$⁠, pressure $p$⁠, and director field $n$ in $Ω$ and at time $t>0$ are^16–19 

with $ρ>0$ denoting the constant mass density and $γ1$ and $γ2$ being combinations of Leslie coefficients defined in (10). We point out that since the system is strongly dissipative, rotational inertia has been neglected in Eq. (2). Moreover, such an equation results from taking the cross product of the original equation of balance of angular momentum by the vector $n$⁠. This has the advantage of explicitly suppressing the Lagrange multiplier associated with the unit director field constraint. The function $WOF$ denotes the Oseen–Frank energy of the liquid crystal, quadratic in the gradients of $n$⁠,

with $k1$⁠, $k2$⁠, $k3>0$⁠, $k2>|k4|$⁠, and $2k1≥k2+k4$ denote the Frank elastic constants. The total energy is

The Cauchy stress tensor $σ$ is the sum of the elastic, viscous $σ^$⁠, and active $σA$ components, respectively,

where

Here, the superimposed dot denotes the material time derivative, that is, $f˙(t,x)=∂f∂t+(v⋅∇)f.$ The Leslie coefficients $αi$⁠, $1≤i≤6$⁠, represent the anisotropic viscosities of the liquid crystal. In particular, $α4$ corresponds to the isotropic or Newtonian viscosity. The parameter $a$ in (8) quantifies the activity of the system, with $a=0$ corresponding to the standard Ericksen–Leslie system for passive liquid crystals.

The active part (8) of the stress tensor accounts for the non-conservative forces generated by the individual fibers that are assumed to be dipolar. Their expressions were obtained from the symmetry of the flow field that they generate, with $a>0$ corresponding to the extensile regime and $a<0$ to the contractile one^1,2,10,11 as illustrated in Fig. 1. In the terminology of swimmers, extensile particles are known as pushers and contractile ones as pullers.

The rate of dissipation function is quadratic on the time-rate quantities $N$ and $A$⁠. It takes the form

The second law of thermodynamics in the form of the Clausius–Duhem inequality reduces to the positivity of the rate of dissipation function, $Δ≥0$⁠. Necessary and sufficient conditions for the latter result in the well-known inequalities,¹⁶ 

Parodi’s relation, a consequence of Onsager’s reciprocal relations in the microscale description of liquid crystals, is an additional assumption of the theory,

This condition renders the rate of dissipation function of a potential for the viscous stress; that is, $σ^=∂Δ∂∇v$⁠. We consider a class of liquid crystals able to align under the effect of flow of a small velocity gradient. This requires that

where $λ$ is known as the flow alignment parameter. It represents the ratio between the extensional and rotational effects of the shear flow, with the former dominating in the case $λ>1$ and, therefore, the director aligns along the flow direction. The tumbling regime corresponds to $λ<1$⁠, with a prevailing rotational couple that prevents $n$ from choosing an aligning direction.^11,20

The behavior of active liquid crystals on the domain boundary may be significantly different from that of their passive counterparts. Whereas actin fibers may anchor and stick to the boundary, some active liquid crystals, such as microtubules, do not align by anchoring on the boundary and can only become oriented by flow. There is no evidence of non-slip behavior in experiments where microtubules are found to slide along the walls. Guided by experiments, slip-free boundary conditions were used in numerical simulations in Ref. 3. Hence, we require

where $ν$ and $τ$ denote the outer unit normal and tangent vectors to the boundary, respectively. For the director field, we impose zero Neumann boundary conditions,

When $Ω$ is a semi-infinite stripe domain of width $2L3$⁠, suppose that the flow and director field are restricted to the $xz$-plane and denote by $ϕ$ the angle between $n$ and the $z$-axis. Then, the conditions (13)–(15) become

respectively.

Next, we formulate the governing equations in terms of dimensionless time and space variables. For this, we choose the positive quantities (⁠$0<L≤L1$⁠), $L1$ representing the length of the channel and $V$ to be the characteristic length and the characteristic velocity, respectively. Moreover, we choose the characteristic viscosity $η$ to be equal to the isotropic viscosity coefficient $α4$ and let $p0>0$ represent the typical pressure. The resulting dimensionless variables are

The quantity $T>0$ represents the typical time scale. We subsequently divide the governing equations, term by term, by the expression $VηL2$⁠. The resulting equations involve well-known dimensionless parameter groups: the Reynolds number $Re$⁠, the Ericksen number $Er$⁠, and the activity parameter $A$⁠, as well as the aspect ratio of the channel, $ℓ$⁠. The quantity $Er$ represents the ratio of the viscous to the elastic torques. Furthermore, we multiply the equation of balance of angular momentum by $L2k12$⁠, which also brings $Er$ into the expression.

We also need the following nondimensional parameters:

where $α~i$ and $γ~i$ are the dimensionless viscosity coefficients and $k~$ is the dimensionless elastic modulus. The list of the dimensionless parameters groups of the model is

In summary, the list of the model parameters is

Likewise, ${v~,n,p~}$ and the Lagrange multiplier, maintaining the unit director constraint, are the unknown fields of the Ericksen–Leslie model.

Since our work deals with the active liquid crystal confined in a channel, in later sections, we will recast the governing equations in terms of two space dimensions.

We suppress the superimposed bar notation and assume that all the variables are already dimensionless. We look for solutions such that the velocity field is unidirectional $v=(U(z,t),0,0)$ and the director angle $ϕ=ϕ(z,t)$⁠. From the component form of dimensionless equations (A1)–(A3) in the supplementary material for velocity and stress, we get the reduced one-dimensional problem^18,21 in (24)–(26),

where

The function $g(ϕ)$ is the dimensionless form of the rate of dissipation $Δ$ of the flow,

Next, we eliminate the pressure $p$ from the equations of balance of linear momentum. Taking the partial derivative of Eq. (25) with respect to $x$⁠, we get

Likewise, taking the derivative with respect to $z$ in Eq. (24) and applying to it the previous result, we get

Integrating once with respect to $z$⁠, we get

where $c1(t)$ is arbitrary. The governing system reduces now to Eqs. (26) and (32). Subsequent use of Eq. (25) determines the pressure.

Let us now look for one-dimensional fields $v=(U(z),0,0)$⁠, $ϕ=ϕ(z)$⁠, and $s=s(z)$⁠. They satisfy the system of equations,

where $c1$ and $c2$ are arbitrary constants obtained in setting $Ut≡0$ in Eq. (32) and further integrating it with respect to $z$⁠. Let us now determine the constants $c1$ and $c2$⁠. We point out that combining Eq. (24) in the steady state case with Eq. (33) yields

We shall take $c1=0$⁠, which corresponds to the assumption that there is no applied external pressure gradient driving the flow. Furthermore, requiring that the velocity gradient vanishes when the activity is equal to zero implies that $c2=0$⁠. The resulting expression of the velocity gradient is then

Rewriting the previous equation in terms of the original dimensional variables, we get

where now the “prime” notation indicates a derivative with respect to the dimensional transverse variable, $z$⁠, and $a$ denotes the activity parameter in (8). This indicates that the shear rate depends on the activity only but not on the width of the channel in agreement with experiments.³ 

Substitution of Eq. (35) into (34) gives

Next, we look for solutions of Eq. (37) such that $ϕ$ and consequently $U′(z)$ are constant. The latter corresponds to the observed states with the linear velocity profile.The angle of orientation of the director is

From the properties of $γ1$ and $γ2$ for rodlike nematics, we observe that $cos2ϕ<0$⁠; therefore,

That is, the angle $0≤π2−ϕ≤π4$ between the director field and the horizontal direction is smaller than $π4$ rad, $π2>π2−ϕ>π4.$ The velocity field is given by

where the constant of integration has been chosen to give the odd profile, illustrated in Fig. 2. Rewriting the previous equation in terms of the original dimensional variables as in (36), we find that the velocity profile does not depend on the channel width. We observe that the solutions show very good agreement with the experimental results. Indeed, in the shear flow regimes, the activity parameter does not directly influence the flow alignment, but it does increase the velocity gradient, i.e., the shear rate.

Note that Eq. (38) shows that the director angle does not depend on the activity parameter $A$ whose dependence enters the expression (39) of the velocity gradient. Both properties are found to be in full agreement with experiments.³ 

From henceforth, we study the stability of the solution $(U(z),ϕ0)$ under perturbations of the form

where $U$ in (39) and $ϕ0$ in (38) are the shear base flow solution and $ϕ0$ is constant. Substituting these expressions into the full two-dimensional system of governing equations (A1)–(A3) in the supplementary material, we obtain the linear system (B1)–(B3) in the supplementary material for the fields $(v1(t,x,z),v3(t,x,z),ϕ1(t,x,z))$⁠.

Moreover, we propose the following exponential expressions of the unknown fields, consistent with those used in the normal mode stability analysis,

where $Ψ$ denotes the stream function of the flow. Also, for the direction angle, we assume that

Here, $ω$ and $c$ are dimensionless complex numbers corresponding to the spatial frequency and to the speed or growth of the perturbation, respectively. Specifically, separating $Ψ$ and $Φ$ into their real and imaginary parts, it follows that $Re(ω)$ represents the spatial oscillatory part of the perturbation and $Im(c)$ corresponds to its time growth rate. From now on, we will restrict ourselves to the case when

Substituting the expressions (40)–(43) into the linear governing equations (B1)–(B3) in the supplementary material, we obtain a linear system (C3) and (C4) in the supplementary material for the new variables. The former is a fourth order linear ordinary differential equation for $ψ$ and the latter is a second order equation for $Φ$⁠.

A relevant quantity in the analysis of flow that also helps to quantify the transition to turbulence is the vorticity vector: $W=∇×v$⁠. Since we are dealing with two-dimensional flow, it reduces to a scalar,

The governing equations for the perturbations are given by (C3) and (C4) in the supplementary material, where their derivation is carried out.

We require that the perturbations do not change the boundary values of the velocity and alignment of the basic solutions. That is, we impose the conditions

The Chebyshev polynomials^22,23 in the variable $z$ form an orthogonal family on the $z$-interval $[−1,1]$⁠. We apply the Chebyshev-QZ algorithm^24,25 to solve the generalized eigenvalue problem (C3) and (C4) in the supplementary material and (45)–(47). We first approximate stream function $ψ(z)$ and angle function $Φ(z)$ by the truncated Chebyshev expansions,

where $Tn(z)$ denotes the nth order Chebyshev polynomial of the first-kind. Our goal is to determine the coefficients $an$⁠, $bn$⁠, and the eigenvalue $c$⁠. To keep the unknown coefficients evenly distributed over the two functions, we choose different orders of truncations for $ψ(z)$ and $Φ(z)$⁠. We collocate the Galerkin truncation at the extrema of the Chebyshev polynomial,

Thus, evaluating the governing equations at these extrema points, we obtain $2N−6$ linear algebraic equations. Six additional relations are provided by the corresponding boundary conditions (45)–(47), which complete the system.

The substitution of the truncated expansions (48) into the boundary conditions yields rows of zeros, which produce a spurious eigenvalue,²⁴ and therefore, we eliminate them. The full system of equations becomes the algebraic eigenvalue problem,

where $x=(a3,…,aN,b2,…,bN−2)T∈C2N−5$ and $AR$⁠, $AI$⁠, $BR$⁠, and $BI$ denote $(2N−5)×(2N−5)$ real matrices. Using the QZ-algorithm of MATLAB, we obtain the eigenvalues and corresponding eigenvectors. Details of the steps leading to the system (50) are given by equations (E1)–(E3) in the supplementary material.

We end this section by listing the values of the Leslie viscosity coefficients used in the simulations. For extensile liquid crystals, we take

The data list to be used in the simulations of contractile liquid crystals is

These provide simple values that still maintain the anisotropy of the viscosity, satisfy the positivity of the rate of dissipation function, and also represent the aligning regime in each class. Moreover, we take the constant order parameter as $s0=1.0.$

The numerical study yields plots of $ω$ with respect to $A$⁠, for different values of $Er$⁠, $ℓ$⁠, and $Re$⁠, showing regions of the spatial frequency domain for which the corresponding perturbation is either stable or unstable, that is, whether it decays or grows in time. The study also yields plots of the streamlines, vorticity, and the director angle of the perturbation fields. The tangent vector field of the former corresponds to the velocity field of the system. We focus on quantitatively understanding the role of the parameter values in determining the instability behavior. Specifically, we assume that the Leslie coefficients $αi$ are fixed and seek how the dimensionless parameters $Er,A,ℓ2$⁠, and $Re$ affect the stability of the shear flow.

Our analysis follows along the lines of many previous investigations of the physical mechanisms that cause either instability or stability in terms of eigenmodes of the linearized system.^26–28 To characterize the stability of the shear flow steady state, we check the growth rate of the dominant unstable eigenmode of the perturbation that affects the system. For unstable systems, the largest value of the linear growth rate and the corresponding wavenumber excite the system and modify the basic state in some essential fashion. On the other hand, stable perturbation modes also modify the system but decay in time.

In this section, we perform a stability analysis of the basic shear flow solution for extensile fibers, that is, in the case $A>0$⁠. For this, we solve the system (50) following the method presented in Sec. VI. The result of the analysis is presented in the graphs of Figs. 3–9. Figures 10 and 11 show the profiles of the perturbing fields. We summarize the results of the computations as follows.

• The main general trend is that increasing, either one of the quantities $Er$⁠, $ℓ$⁠, or $A$⁠, prompts unstable behavior of the shear flow regime.

• In particular, Fig. 3 shows that, for fixed $Er=1000$ and $Re=5$⁠, the critical value of the activity number $A$⁠, above which the system becomes unstable, decreases with $ℓ$⁠, with $Ac≈350$ dimensionless units for $ℓ=0.1$ in Fig. 3(a), $Ac≈100$ for $ℓ=0.2$ in Fig. 3(b), and $Ac≈20$ for $ℓ=0.4$ in Fig. 3(c).

• Figures 4(a)–4(c) show the analogous behavior but with respect to increasing $Er$⁠. The increase of $Er$⁠, while keeping the other parameter fixed, also promotes unstable behavior. For $Re=0.1$⁠, $ℓ=0.2$⁠, the threshold values of $A$ above which the flow becomes unstable are found to be $A≈990$ for $Er=100,$ $A≈200$ for $Er=500$⁠, and $A≈100$ for $Er=1000$⁠.

• Figure 5 shows that increasing $Re$ may not decrease the critical activity number threshold leading to instability, but it does increase the unstable frequency range as shown in Fig. 5(b). Figure 5(c) reiterates the role of the channel width in promoting instability, showing the instability threshold of $A≈100$⁠, even for $Re=5$ and $Er=250$ but with $ℓ=0.4.$

• Figure 6 shows the regions of stability in the $Aω$-plane, with large Reynolds number $Re=100$ but small Ericksen number $Er=5$⁠. For these results, the width parameter $ℓ=0.2$ is fixed.

• Figure 7 illustrates the effect of the channel width in the $ℓΩ$-plane when activity is fixed for the system. Figure 7(a) shows the regions of stability for $Re=5$ and $A=10$⁠. Figure 7(b) shows the results for $Re=0.1$ and $A=100$⁠. $Er=1000$ is fixed for both plots.

• Figure 8 shows the growth rate given by $maxIm(c)$ with respect to $A$⁠, for the wavelengths $ω=50$⁠, $75$⁠, and $100$⁠. The remaining parameters are chosen as $Re=5$⁠, $Er=1000$⁠, and $ℓ=0.2$⁠. We find that the system is stable [$maxIm(c)<0$] for small $A$⁠, with a stability threshold between $A=100$ (for $ω=50$⁠) and $A=150$ (for $ω=100$⁠) dimensionless units. For small frequencies, the growth rate shows a parabolic profile, reaching a maximum whose $A$-location increases with frequency, with a possible return to stability at higher $A.$ For instance, for $ω=50$⁠, the system changes from stable to unstable at $A≈100$ and then regaining stability at $A≈500$⁠. The profile becomes linear at $ω=100$⁠.

• Figure 9 shows the graphs of the growth rate with respect to the frequency, at three different activity values, $A=150$⁠, 200, and 250. We observe the nonconvex shape of the graphs, for small frequencies, reaching a positive maximum. All the graphs tend to neutral stability, with increasing $ω$⁠, following a profile nearly independent of $A$⁠.

• Figure 10 illustrates the form of the perturbation fields that contribute to the shear flow. We present contour plots of the stream function, the angle, and the vorticity contour (Fig. 11).

Extensile fibers are known to form rodlike nematic liquid crystalline phases due to their elongated molecular shapes. The material that we study consists of self-propelled elongated fiber units formed by bundled microtubules²⁹ that are powered by adenosine triphosphate (ATP)-consuming kinesin 23. The confinement model of fibers provides an appropriate framework to study biological systems, such as actin, in the extracellular network as well as in-cell structures.³⁰ 

Our flow model is based on the Ericksen–Leslie equations of liquid crystals with an added activity stress term.^16,31 Guided by the active matter confining experiments reported in Ref. 3, we study laminar flow of active liquid crystals in a channel. Specifically, we quantitatively analyze the mechanisms of instability of well-aligned shear flows, with the experimentally observed linear velocity profile and characterize them in terms of relevant dimensionless parameters, the activity and Ericksen numbers, and the channel aspect ratio. We also find that the vorticity and speed of the perturbing flow increase as the activity number $A$ increases. A numerical method based on the discretization of the linear system by Chebyshev polynomials is used in the stability analysis, specifically in solving the underlying spectral problems. The analysis illustrates how geometrical confinement tends to control active flows, replacing bulk chaotic behavior with laminar flow configurations.

Our next goal is the measurement of angular profiles of the fibers by analyzing images obtained from experiments. This will provide information on the ratio of the viscosity coefficients $γ1γ2$⁠, quantitatively measure the relation between rotation and extension rates of the fibers, and provide some needed information on the Leslie coefficients of the material. Another goal is the analysis of higher activity and wider channel regimes, where defects cannot be neglected, focusing, in particular, on the dancing disclination state, where formation and subsequent pairwise annihilation of degree (⁠$±12$⁠) defects drive the flow.

See the supplementary material for detailed derivation of model equations and numerical implementations.

M. C. Calderer gratefully acknowledges the support of the National Science Foundation (NSF) through Grant Nos. DMS-1435372 and DMS-DMREF 1729589. D. Golovaty acknowledges funding from the NSF grant (No. DMS-1729538). L. Yao acknowledges support from the NSF grant (No. DMS-1852597). L. Zhao acknowledges support from the AWM-NSF mentor grant. J. Ignés-Mullol and F. Sagués acknowledge funding from the Ministerio de Economia, Industriay Competitividad, Spain (Project No. FIS2016-78507-C2-1-P, Agencia Estatal de Investigación-European Regional Development Fund).

Data sharing is not applicable to this article as no new data were created or analyzed in this study. The numerical algorithm is available in the supplementary material. The data that support the findings of this study are available from the corresponding author upon reasonable request.