In this article, we study shear flow of active extensile filaments confined in a narrow channel. They behave as nematic liquid crystals that we assumed are governed by the Ericksen–Leslie equations of balance of linear and angular momentum. The addition of an activity source term in the Leslie stress captures the role of the biofuel prompting the dynamics. The dimensionless form of the governing system includes the Ericksen, activity, and Reynolds numbers together with the aspect ratio of the channel as the main driving parameters affecting the stability of the system. The active system that guides our analysis is composed of microtubules concentrated in bundles, hundreds of microns long, placed in a narrow channel domain, of aspect ratios in the range between and dimensionless units, which are able to align due to the combination of adenosine triphosphate-supplied energy and confinement effects. Specifically, this work aims at studying the role of confinement on the behavior of active matter. It is experimentally observed that, at an appropriately low activity and channel width, the active flow is laminar, with the linear velocity profile and the angle of alignment analogous to those in passive shear, developing defects and becoming chaotic, at a large activity and a channel aspect ratio. The present work addresses the laminar regime, where defect formation does not play a role. We perform a normal mode stability analysis of the base shear flow. A comprehensive description of the stability properties is obtained in terms of the driving parameters of the system. Our main finding, in addition to the geometry and magnitude of the flow profiles, and also consistent with the experimental observations, is that the transition to instability of the uniformly aligned shear flow occurs at a threshold value of the activity parameter, with the transition also being affected by the channel aspect ratio. The role of the parameters on the vorticity and angular profiles of the perturbing flow is also analyzed and found to agree with the experimentally observed transition to turbulent regimes. A spectral method based on Chebyshev polynomials is used to solve the generalized eigenvalue problems arising in the stability analysis.
This article presents a comprehensive stability analysis of the aligning shear flow of active fibers driven by adenosine triphosphate (ATP)-supplied energy. The study addresses extensile fibers confined in narrow channels. The material and the experimental settings provide a model environment to study confinement of in vivo biological systems that may prove useful in medical research. The setup of the problem and the results of the analysis closely follow experiments of microtubules as reported in recent articles. In our study, to make predictions for experiments, the activity level, the ratio of the viscous to the elastic torques, and the channel geometry are quantitatively characterized as the mechanisms that trigger instability, and their role on the vorticity profiles of the perturbing flow is analyzed. A spectral method is applied to solve the generalized eigenvalue problems arising in the stability analyses.
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.3 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.4 analyzed the stability of an aligned state of an active polar gel. Their results were later confirmed in experiments with spindle-shaped cells.5
The active system studied in earlier work3 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 mm in length, m in depth, with a width range between – 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 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 m and (ii) a dancing disclination regime for m, with pairs of defects of fractional degree that form and annihilate on the boundary. The events of defect annihilation drive a flow in the channel, as reported in an earlier article.6 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.7 The variables of the model include the velocity field , the pressure associated with the incompressibility constraint, and the unit nematic director , 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 and Ericksen numbers, are defined for passive flows. The activity parameter quantifies the activity level in the system. Although is significantly smaller than other nondimensional parameters, it influences the stability behavior of the active system in subtle ways. The 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 .3 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 (for tumbling regimes ). 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.3
Several works on active liquid crystals found in the literature use the Beris–Edwards model based on the evolution of the order tensor rather than the director of the Ericksen–Leslie theory.8 represents a symmetric, traceless second order tensor whose two independent eigenvectors are the director fields of the theory. They reduce to a single eigenvector in the uniaxial case, with admitting the representation where the scalar corresponds to the single, independent eigenvalue. Our work would correspond to a constant value the positiveness agreeing with dealing with rod-like liquid crystals. A recent article9 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 the corresponding growth rate. We find that a perturbation with frequency is unstable for . The threshold 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 -plane indicate threshold values of at which the system becomes unstable. The critical value decreases with increasing values of and and with the range of unstable frequencies also increasing. The growth rate profile calculated with respect to shows a parabolic shape at low frequencies, becoming linear at a threshold frequency value. The profile of the growth rate with respect to starts being positive for all values of , becoming at a threshold value of that appears to be independent of . With further increase of , reaches a positive maximum, subsequently tending to zero as increases, a feature that is also independent of . The instability profiles obtained from the analysis qualitatively follow the experimentally reported trends,3 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.
II. ACTIVE LIQUID CRYSTALS: THE ERICKSEN–LESLIE EQUATIONS
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 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 , pressure , and director field in and at time are16–19
with denoting the constant mass density and and 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 . This has the advantage of explicitly suppressing the Lagrange multiplier associated with the unit director field constraint. The function denotes the Oseen–Frank energy of the liquid crystal, quadratic in the gradients of ,
with , , , , and denote the Frank elastic constants. The total energy is
The Cauchy stress tensor is the sum of the elastic, viscous , and active components, respectively,
Here, the superimposed dot denotes the material time derivative, that is, The Leslie coefficients , , represent the anisotropic viscosities of the liquid crystal. In particular, corresponds to the isotropic or Newtonian viscosity. The parameter in (8) quantifies the activity of the system, with 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 corresponding to the extensile regime and to the contractile one1,2,10,11 as illustrated in Fig. 1. In the terminology of swimmers, extensile particles are known as pushers and contractile ones as pullers.
A. The Leslie coefficients and the dissipativity of the system
The rate of dissipation function is quadratic on the time-rate quantities and . 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, . Necessary and sufficient conditions for the latter result in the well-known inequalities,16
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, . 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 and, therefore, the director aligns along the flow direction. The tumbling regime corresponds to , with a prevailing rotational couple that prevents from choosing an aligning direction.11,20
B. Boundary conditions
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 , suppose that the flow and director field are restricted to the -plane and denote by the angle between and the -axis. Then, the conditions (13)–(15) become
In the calculations for the stability analysis, we impose Dirichlet boundary conditions on the perturbation angle rather than the Neumann condition (18). This has been done to simplify the stability analysis, but one would naturally expect different results for other boundary conditions.
III. SCALING AND NON-DIMENSIONALIZATION
Next, we formulate the governing equations in terms of dimensionless time and space variables. For this, we choose the positive quantities (), representing the length of the channel and to be the characteristic length and the characteristic velocity, respectively. Moreover, we choose the characteristic viscosity to be equal to the isotropic viscosity coefficient and let represent the typical pressure. The resulting dimensionless variables are
The quantity represents the typical time scale. We subsequently divide the governing equations, term by term, by the expression . The resulting equations involve well-known dimensionless parameter groups: the Reynolds number , the Ericksen number , and the activity parameter , as well as the aspect ratio of the channel, . The quantity represents the ratio of the viscous to the elastic torques. Furthermore, we multiply the equation of balance of angular momentum by , which also brings into the expression.
We also need the following nondimensional parameters:
where and are the dimensionless viscosity coefficients and 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, 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.
IV. SHEAR FLOW OF THE ERICKSEN–LESLIE MODEL
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 and the director angle . From the component form of dimensionless equations (A1)–(A3) in the supplementary material for velocity and stress, we get the reduced one-dimensional problem18,21 in (24)–(26),
The function is the dimensionless form of the rate of dissipation of the flow,
Next, we eliminate the pressure from the equations of balance of linear momentum. Taking the partial derivative of Eq. (25) with respect to , we get
Likewise, taking the derivative with respect to in Eq. (24) and applying to it the previous result, we get
Integrating once with respect to , we get
A. Steady state flow
Let us now look for one-dimensional fields , , and . They satisfy the system of equations,
where and are arbitrary constants obtained in setting in Eq. (32) and further integrating it with respect to . Let us now determine the constants and . We point out that combining Eq. (24) in the steady state case with Eq. (33) yields
We shall take , 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 . 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, , and 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.3
Next, we look for solutions of Eq. (37) such that and consequently 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 and for rodlike nematics, we observe that ; therefore,
That is, the angle between the director field and the horizontal direction is smaller than rad, 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.
V. STABILITY ANALYSIS
From henceforth, we study the stability of the solution under perturbations of the form
where in (39) and in (38) are the shear base flow solution and 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 .
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 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 represents the spatial oscillatory part of the perturbation and 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: . 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
VI. THE NUMERICAL METHOD
The Chebyshev polynomials22,23 in the variable form an orthogonal family on the -interval . We apply the Chebyshev-QZ algorithm24,25 to solve the generalized eigenvalue problem (C3) and (C4) in the supplementary material and (45)–(47). We first approximate stream function and angle function by the truncated Chebyshev expansions,
where denotes the nth order Chebyshev polynomial of the first-kind. Our goal is to determine the coefficients , , and the eigenvalue . To keep the unknown coefficients evenly distributed over the two functions, we choose different orders of truncations for and . We collocate the Galerkin truncation at the extrema of the Chebyshev polynomial,
Thus, evaluating the governing equations at these extrema points, we obtain 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,24 and therefore, we eliminate them. The full system of equations becomes the algebraic eigenvalue problem,
where and , , , and denote 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
The numerical study yields plots of with respect to , for different values of , , and , 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 are fixed and seek how the dimensionless parameters , and 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.
VII. STABILITY ANALYSIS
In this section, we perform a stability analysis of the basic shear flow solution for extensile fibers, that is, in the case . 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 , , or , prompts unstable behavior of the shear flow regime.
In particular, Fig. 3 shows that, for fixed and , the critical value of the activity number , above which the system becomes unstable, decreases with , with dimensionless units for in Fig. 3(a), for in Fig. 3(b), and for in Fig. 3(c).
Figures 4(a)–4(c) show the analogous behavior but with respect to increasing . The increase of , while keeping the other parameter fixed, also promotes unstable behavior. For , , the threshold values of above which the flow becomes unstable are found to be for for , and for .
Figure 5 shows that increasing 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 , even for and but with
Figure 6 shows the regions of stability in the -plane, with large Reynolds number but small Ericksen number . For these results, the width parameter 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 and . Figure 7(b) shows the results for and . is fixed for both plots.
Figure 8 shows the growth rate given by with respect to , for the wavelengths , , and . The remaining parameters are chosen as , , and . We find that the system is stable  for small , with a stability threshold between (for ) and (for ) dimensionless units. For small frequencies, the growth rate shows a parabolic profile, reaching a maximum whose -location increases with frequency, with a possible return to stability at higher For instance, for , the system changes from stable to unstable at and then regaining stability at . The profile becomes linear at .
Figure 9 shows the graphs of the growth rate with respect to the frequency, at three different activity values, , 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 .
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 microtubules29 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.30
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 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 , 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 () 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.