The Taylor–Couette motor: spontaneous flows of active polar fluids between two coaxial cylinders

We study the dynamics of active polar fluids in a Taylor–Couette geometry where the fluid is confined between two rotating coaxial cylinders. This system can spontaneously generate flow fields and thereby set the two cylinders into relative rotation either by spontaneous symmetry breaking or via asymmetric boundary conditions on the polarization field at the cylinder surfaces. In the presence of an externally applied torque, the system can act as a rotatory motor and perform mechanical work. The relation between the relative angular velocity of the cylinders and the externally applied torque exhibits rich behaviors such as dynamic instabilities and the coexistence of multiple stable steady states for certain ranges of parameter values and boundary conditions.


Introduction
Living matter is internally driven far from thermodynamic equilibrium by active molecular processes such as the action of motor molecules. These motors consume a chemical fuel, adenosinetriphosphate (ATP), to generate directed movements and mechanical work by interacting with structurally polar filaments. In living cells, large numbers of motors and filaments collectively generate dynamic processes such as cell division, cell locomotion and organelle transport in cells. The network of filaments and motors is called the cytoskeleton and provides a key example of an active complex fluid with unconventional material properties. Other examples are suspensions of microswimmers or collections of many cells, for example, in developing tissues.
Active fluids often exhibit long-range nematic or polar order. For example, in crawling cells, cytoskeletal filaments orient typically in the direction of motion. Similarly, the direction of motion of microswimmers defines a polarity vector. Macroscopic polar order emerges when many filaments align or swimmers form swarms [1]. The interplay of actively generated flow fields with the orientation of polar order can render ordered states unstable [2,3]. Instabilities of active polar fluids can have interesting physical consequences and might also be responsible for observed 'low Reynolds number turbulence' in bacterial suspensions [4]. Furthermore, the polarity dynamics of active fluids could be responsible for physiologically relevant structures such as contractile rings in dividing cells [5][6][7], flows during cellular wound healing [8] or the reorientation of planar cell polarity in tissues [9].
In the hydrodynamic limit, i.e. on large length and long time scales, the generic properties of visco-elastic active polar gels are described by generic equations. Hydrodynamic modes are the result of conservation laws and of broken continuous symmetries [10]. Hydrodynamic equations can be obtained systematically by first identifying conjugate pairs of thermodynamic fluxes and forces. Constitutive material relations can be expressed by writing all linear coupling terms permitted by symmetry, respecting Onsager reciprocity relations and the signature 3 with respect to time-reversal [11][12][13]. Note that similar approaches have also been used to describe bird flocks [14], swarms of hydrodynamically interacting swimmers [2], active nematic fluids [15,16] or active solids [17][18][19]. Furthermore, an extension to multi-component active fluids has been developed [20].
A large number of instabilities of externally driven flows are known, some of which have been analyzed in remarkable detail [21]. Among the best-studied instabilities are those of a fluid flowing in the interstice between two rotating coaxial cylinders, known as the Taylor-Couette system. In this context, complex fluids such as ferrofluids in the presence of a magnetic field or visco-elastic fluids have also been considered. For the latter, notably, finite-amplitude instabilities of linearly stable states have been found even at low Reynolds numbers if elastic stresses decay sufficiently slowly [22,23].
Here, we study the flow of an active polar fluid in the Taylor-Couette geometry to highlight some of the interesting rheological properties of active fluids in a setup that might be experimentally accessible. This requires an analysis of the angular momentum flux, which is related to the externally applied torque. In particular, we discuss the role of the Ericksen stress in the torque balance and extend the discussion of the so-called magic spiral [24,25] to active fluids. Our main result is that an active polar gel in a Taylor-Couette system becomes a rotatory motor which can generate spontaneous rotation against external torques. We distinguish between symmetric motors where motion occurs by spontaneous symmetry breaking and asymmetric motors in which symmetry is broken by boundary conditions. Our analysis is also relevant to biological processes, notably for active flows in the cytoskeleton. Moreover, it might serve as a starting point for the design of biomimetic rotational motors. This paper is organized as follows. In section 2, we present a derivation of the hydrodynamic equations of active polar fluids, paying particular attention to the conservation of angular momentum. The analysis of the Taylor-Couette motor is presented in section 3. The paper concludes with a discussion of our results in section 4.

Hydrodynamic equations for an active polar fluid
We discuss the hydrodynamic equations of an isothermal active polar fluid. Following the logic outlined in [11,12,24], we express the conservation laws and introduce the relevant hydrodynamic variables. Based on the conjugate pairs of fluxes and forces in the entropy production, we write the dynamic equations of the system. In the Taylor-Couette geometry discussed here, angular momentum fluxes play an important role and are therefore specifically discussed.

Conservation laws
In a fluid, the conserved quantities are mass, momentum and angular momentum. Conservation of mass implies that where ρ denotes mass density. This equation defines the center of mass velocity v α . The conservation of momentum is given by where g = ρv is the momentum density. Here, the stress tensor σ tot denotes the total momentum flux and f ext denotes the density of external bulk forces. Greek indices take the values 1, . . . , d, where d is the spatial dimension. Angular momentum conservation is expressed by where αβ = r α g β − r β g α is the angular momentum density and τ ext αβ is the density of the external bulk torques. The total angular momentum flux M tot αβγ has contributions r α σ tot βγ − r β σ tot αγ resulting from momentum fluxes. In complex fluids, however, there exist in general additional contributions M αβγ to the angular momentum flux that are related to molecular anisotropies: Here, we are interested in polar fluids with a local polar anisotropy that is described by the polarization vector p. This polarization field is defined as the average of microscopic polarity in small volume elements of the system. For a polar fluid, the expression for the intrinsic angular momentum flux M is given in appendix A. Using the force and torque balance equations (2) and (3), we see that the antisymmetric part of the total stress in general does not vanish and obeys Here σ tot,a αβ = (σ tot αβ − σ tot βα )/2, see also [44] 5 . We only consider situations where the total external force and torque applied to the system is zero and where external bulk forces and torques vanish, i.e. f ext α = 0 and τ ext αβ = 0. In this case, external forces σ tot αγ dS γ and torques can still act locally on the boundary ∂ V at a surface element dS γ pointing outward. Vanishing external force and torque imply ∂ V σ tot αγ dS γ = 0 and Note that even if the total flux M tot αβγ vanishes at the boundary, the contributions due to M αβγ and σ tot αβ can be non-zero according to equation (4).

Free energy and hydrostatic stress
We consider a fluid that is locally in thermal equilibrium, but globally evolves in time. Local equilibrium implies that the free energy density f is well defined. We write the free energy density as

5
The first term is the kinetic energy of center-of-mass motion. The second term is the free energy density in the local rest frame. Here, n = ρ/m is the particle density and m denotes the molecular mass. We define the field h α = −δ F/δp α conjugate to the polarization field, and the chemical potential µ = ∂ f 0 /∂n, conjugate to the particle density. The total free energy is given by F = d d x f . The hydrostatic or Ericksen stress, can be obtained by considering variations of the free energy as shown in appendix B and is a generalization for anisotropic fluids of the hydrostatic pressure. Exploiting the translation invariance of the system, we obtain the Gibbs-Duhem relation, which links the intensive thermodynamic variables, The hydrostatic stress tensor has an antisymmetric component, This equation is a consequence of the rotational invariance of f 0 . In equilibrium, when h α = 0 and σ tot αβ = σ e αβ , equation (11) becomes a special case of equation (5). Note that even though M αβγ and σ tot αβ do not necessarily vanish at equilibrium, the total angular momentum flux M tot αβγ must be zero if external forces and torques are absent. This implies that Consequently, the torques due to polarization field and due to the Ericksen stress cancel each other. This fact is illustrated by the so-called magic spiral [24,25].

Entropy production and constitutive equations
Starting from an expression for the rate of entropy production, we identify thermodynamic forces and the conjugate thermodynamic fluxes and express constitutive relations between these fluxes and forces. The rate of entropy production˙ for an isothermal system [11] can be written as Here T denotes the temperature, the dots indicate time-derivatives and J F denotes changes in free energy due to fluxes at the boundary. Moreover, W denotes the work exerted by external forces and torques. Using the conservation laws, we can evaluateḞ and find the entropy production rate, Here, we have assumed that the system is driven by a chemical reaction of a fuel which is provided everywhere at constant concentration by a reservoir. In the biological context this fuel is ATP, which is hydrolyzed to adenosinediphosphate (ADP) and inorganic phosphate P. The chemical free energy transduced per fuel molecule µ = µ ATP − µ ADP − µ P is the difference of the chemical potentials of ATP and its hydrolysis products. The reaction rate describing the number of fuel molecules hydrolyzed per unit time is denoted by r . Therefore, the rate of entropy production due to ATP hydrolysis is r µ/T . We separate σ tot − σ e into a symmetric part and an antisymmetric part where equations (5) and (11) have been used. Note that equation (16) differs from the expression given for σ a αβ in [13] by an overall sign. Equation (14) then becomes Here u αβ = 1 2 (∂ α v β + ∂ β v α ) are the components of the strain rate tensor, while the convected co-rotational derivative of the polarization vector is given by where ω αβ = 1 2 (∂ α v β − ∂ β v α ) is the vorticity of the velocity field. From equation (17) we thus identify the thermodynamic fluxes σ αβ , Dp α /Dt and r as well as the respective corresponding thermodynamic forces u αβ , h α and µ. Expanding the fluxes in terms of the forces, respecting the symmetries of the system and Onsager's reciprocity principle, we obtain the generic constitutive equations of the active polar fluid [13]: The viscosities η and γ denote the resistance of the fluid to internal rearrangements. The coefficients ν 1 ,ν 1 capture the coupling of the polar distortion field to flows. Finally, the coefficients denoted by ζ ,ζ , ζ describe the active stresses. Positive active stresses are called expansive, while negative active stresses are called contractile. The values of these phenomenological constants either need to be measured or can be derived by coarse-graining microscopic theories [26][27][28][29][30].
The equations of motion of the system are given by the continuity equations (1) and (2) and by equation (20) for the time evolution of p. Using equations (11) and (5) and the Gibbs-Duhem relation (10), the force balance becomes

The Taylor-Couette motor
We now discuss the flows of an active polar fluid in the Taylor-Couette geometry at low Reynolds number. The fluid is confined to the space between two impermeable concentric cylinders of radii R + and R − , as shown in figure 1(a). We will discuss situations where the two cylinders are either stationary or rotate relative to each other at a rate ω. At low Reynolds number we can choose to keep the inner cylinder fixed without loss of generality. The rotation rate of the outer cylinder then equals ω. We only consider cases that are invariant with respect to rotations around and translations along the cylinders' long axes. Furthermore, the fluid is incompressible, which together with translational and rotational invariance imposes v r = 0. We only consider cases where the vector p is confined to the r -θ -plane. Furthermore, since the magnitude of p is not a hydrodynamic variable, we only consider the orientation dynamics of the polarity vector. The orientational dynamics is captured by choosing p 2 = 1. Therefore, the polarity p can be expressed in terms of the polarization angle ψ such that p r = cos ψ and p θ = sin ψ, see figure 1(b).
One can apply torques + and − per unit axial length on the outer and the inner cylinder surfaces, respectively, such that the total torque on the surface vanishes, + + − = 0. Using equation (6), which is independent of r for R − r R + because of angular momentum conservation. Using equations (4) and (16), this implies that

Equation of motion and boundary conditions
We now present the equations of motion governing the system in polar coordinates. The only non-vanishing component of the strain rate tensor is u θr . From equation (19) we obtain where h = h · p and h ⊥ = h r p θ − h θ p r . Projecting equation (20) onto the directions parallel and perpendicular to the polarization vector yields The parallel component of the distortion field h = p · h acts as a Lagrange multiplier to impose the constraint p 2 = 1. The perpendicular component h ⊥ is the torque exerted by the polarization field and is given by the expression h ⊥ = −δ F/δψ. The force balance equation (22) becomes where we have neglected inertial terms. Integration of equation (28) yields where the integration constant − is the torque applied to the inner cylinder as can be seen by comparing equations (29) and (24). Combining equations (25)- (27) and (29) we obtain the equation of motion for the polarization field: whereζ = 1 2 (ζ + γ ν 1 λ 1 ). For a given polarity field ψ(r ), we can determine the flow field v θ by integrating equation (25).
The boundary conditions that we consider are the following. We prescribe the polarization angle at both cylinder surfaces, ψ(R − ) = ψ(R + ) = ψ 0 . We impose no-slip boundary conditions for the velocity of the fluid at both cylinders. At the fixed inner cylinder, this implies v θ (R − ) = 0, at the outer cylinder v θ (R + ) = ω R + . With u θr = (r/2)∂ r (v θ /r ), the rotation rate ω reads At the outer cylinder we have the choice of two ensembles: (i) prescribed torque + = − − or (ii) prescribed rotation rate ω. We will discuss both cases.

Equilibrium steady states
In equilibrium, all thermodynamic fluxes and forces vanish. In particular, h ⊥ = 0 and µ = 0. In addition, equilibrium requires that − = 0 and that no movements occur, v θ = 0. The polarization field ψ of an equilibrium steady state is obtained by solving the equation h ⊥ = 0, with the specified boundary conditions. Choosing the standard Landau-de Gennes form for the free energy of the polarization field, we have where primes indicate derivatives with respect to r . In this expression K is the splay elastic modulus and K + δ K is the bend elastic modulus. If the boundary conditions are chosen accordingly, equilibrium steady states with a constant ψ(r ) = ψ 0 exist for certain values of ψ 0 . These special solutions are 'asters' with ψ 0 = 0 or π and 'vortices' with ψ 0 = ±π/2.

Non-equilibrium steady states
We now discuss how the equilibrium situation is changed as µ becomes non-zero. We start by analyzing the stability of the symmetric aster and vortex states with boundary condition (i) and a prescribed torque − = 0.

No external torque at the outer surface.
We write ψ(r, t) = ψ 0 + δψ(r, t) with δψ = 0 at both boundaries. The time evolution described by equations (30) can be linearized around a steady state with ψ = ψ 0 . For aster solutions, i.e. ψ 0 = 0, we find that Using the separation ansatz δψ(r, t) = δψ(r )e st , we can solve for δψ(r ) where J n and Y n are Bessel functions of the first kind with n 2 = δ K /(K + δ K ) and Asters are stable if the real parts of all growth exponents s are negative. The possible values k i of k with i ∈ N are restricted by the boundary conditions δψ(R − ) = δψ(R + ) = 0. Explicitly, the values k i satisfy the equation We find that k 2 i and the corresponding values of s are real and we order these values such that k 2 1 < k 2 2 < · · · < k 2 i < · · · . By setting s = 0 in equation (36), we can thus determine a critical valueζ µ c for the activity, such that asters are unstable ifζ µ <ζ µ c . Explicitly, Following the same logic it is also possible to determine the stability of vortices. We find that the activity threshold for the stability of vortices is The values l i satisfy Y m (l i R + )J m (l i R − ) − J m (l i R + )Y m (l i R − ) = 0 with m 2 = −δ K /K and we order their values such that l 2 1 < l 2 2 < · · · < l 2 i < · · · . The limit R − → 0 corresponds to point defects that were studied in [13]. Note that our results differ slightly since the sign of the antisymmetric part of the deviatory stress was chosen incorrectly in [13].
The steady states of equation (30) beyond the instability can be obtained numerically by a shoot-and-match algorithm. In these states, the fluid flows spontaneously and the outer cylinder rotates. We present the bifurcation diagram of the Taylor  determined by spontaneous symmetry breaking. Therefore, two stable branches exist in the bifurcation diagram. We have explicitly determined the linear stability of these and all following states numerically. Moving further away from the bifurcation point, additional unstable steady states appear. These states correspond to l 2 2 and present an additional node for ψ(r ). Further unstable steady states corresponding to l 2 i with i > 2 exist for yet larger values of |ζ µ| and the number of nodes of the corresponding ψ(r ) is i − 1.

3.3.2.
The stalled system. The emergence of spontaneous rotations indicates that the system can act as a rotatory motor. Its stall torque is determined by setting ω = 0. To investigate the stability of the corresponding steady state, we employ again equation (30), but now − depends on time. Its value is determined by simultaneously solving equation (25) for v θ and imposing v θ (R + ) = 0. For the stationary vortex solution with ψ 0 = ±π/2, we again find an instability at a critical value ζ µ c at which spontaneous flows occur, as is shown in the bifurcation diagram figure 3. As in the case − = 0, the direction of the flow depends on the initial conditions via spontaneous symmetry breaking. In contrast to the case considered above, the solutions with ω = 0 are characterized by non-vanishing torques − .

The relation between the rotation rate and applied torque
We determine the relation between applied torque − and rotation rate ω by numerically solving for the stationary solutions of equations (25) and (27). In the passive caseζ µ = 0, the rotation rate ω is always proportional to the torque ω ∝ − and vanishes for − = 0. Our results for finiteζ µ are displayed in figure 4. For a sufficiently smallζ µ, the rotation rate  ω increases monotonically as a function of − , see figure 4(a). Forζ µ =ζ µ c −20KR −2 + , a critical point appears. Beyond this point, three solutions coexist forζ µ −20KR −2 + , two of which are stable (solid lines) and one is unstable (dashed line); see figure 4(b). As a consequence, the rotation rate exhibits discontinuous changes as a function of torque and hysteresis occurs. Ifζ µ is decreased further, a more complex set of branches of unstable steady states emerges, see figures 4(c) and (d). In all cases there are two stable branches and a regime of hysteresis. Note that the stable branches always have positive slope, while the slope of unstable branches can have either sign. This relation between torque and rotation rate is very similar to the properties of the force-velocity relation of symmetric collections of molecular motors [31,32].

Variation of the boundary conditions of the polarization field
Changing the polarization angle ψ 0 at the boundaries leads in general to a non-symmetric motor, which has a spontaneous rotation rate even in the absence of a torque. In figure 5, we show the steady state rotation rate ω for − = 0 as a function of ψ 0 for different values of the active stress. This figure shows that the system rotates spontaneously except for ψ 0 = 0 (aster) and ψ 0 = π/2 (vortex). As the magnitude of the active stress is increased, unstable branches appear and the rotation rate exhibits discontinuous changes and hysteresis. The vortex becomes unstable for lower values of the active stress than the aster, see figure 5(b). Figure 5(c) reveals the subsequent emergence of additional unstable branches as well as the instability of the aster.

Conclusions and discussion
In summary, we analyzed the dynamics of an active polar fluid confined in the Taylor-Couette geometry between two coaxial cylinders. We found that similar to active polar fluids confined to a straight channel [3], the active stress can induce spontaneous flows. The circular flow in the Taylor-Couette system drives a rotation of the outer cylinder relative to the inner cylinder even in the absence of an externally applied torque. If the polarization angle at the boundary is either parallel or perpendicular to the cylinder wall, relative motion between the two cylinders occurs by spontaneous symmetry breaking if the contractile active stress exceeds a critical value. Beyond this instability two steady states with opposite rotation rates coexist. For any other boundary condition imposed on the polarization, rotation in the absence of torque occurs for non-zero active stress. Here, we have restricted our analysis to states that are invariant with respect to rotations around and translations along the cylinder axis. Since we know that in the passive system rotationally symmetric states are stable for sufficiently low Reynolds numbers, we expect that the rotationally symmetric active steady states discussed above are stable for small enough active stresses. Additional instabilities to those discussed here could therefore occur at higher values of active stresses. Notably, active analogues of Taylor vortices, which break rotational and translational invariance, might exist. The classical Taylor-Couette system displays a large variety of flow patterns, including turbulent flow for very large Reynolds numbers. It will be interesting to see which of these structures can be generated in the low Reynolds number regime by active stresses. In particular, flow patterns that break the axial and rotational symmetries we imposed in this study should be analyzed. In view of 'low Reynolds number turbulence' observed in bacterial suspensions [4], one might also be interested in investigating the possibility of spatiotemporal chaos in the Taylor-Couette motor. Considering the finite relaxation time of elastic stresses exhibited by a number of active gels, the system might also be capable of spontaneously generating oscillations. Consequently, the classical system of two rotating coaxial cylinders is a promising geometry to study also experimentally the rich dynamic behavior of active gels.
If an external torque is applied between the two cylinders, the system behaves as a rotary motor that can perform mechanical work. In the special cases where the polarization is parallel or perpendicular to the boundary, this motor only works beyond the instability and can act in both senses of rotation. We find dynamic transitions between states of different rotation rates and hysteresis. These behaviors are very similar to those found for collections of symmetric coupled molecular motors where motion occurs via spontaneous symmetry breaking [31,32]. The behaviors of the Taylor-Couette motor for other choices of boundary conditions can be compared to those found for asymmetric collections of motors. Indeed, as already described for collective motors we find here dynamic transitions between states of different rotation rates at finite applied torques. Coupling to a rotational spring would thus lead to spontaneous oscillations [33][34][35]. Due to the possible appearance of additional unstable states, such oscillations may exist only in a finite range of active stresses.
The Taylor-Couette motor described in this paper could be realized experimentally using actin gels with myosin motors in small Taylor-Couette rheometers. We will now estimate the orders of magnitude of the quantities characterizing the system's dynamics. First, we note that the active stress together with the elastic constant K , which has the dimension of a force, see equation (32), introduce a length scale ξ = (K /ζ µ) 1/2 , see also [3]. The instability criterion derived above can be expressed in terms of this length scale: for ξ R + − R − spontaneous motion occurs in a symmetric system. The critical active stress beyond which spontaneous motion occurs is thus of the order of K /R 2 + for which ξ R + − R − . In passive systems, the value of K has been estimated using dimensional analysis, as K ∼ U/a where U and a are a characteristic interaction energy and length, respectively [24]. Using U = 10 pN nm and a = 1 nm implies a value of K of the order of several pN, which is consistent with estimates provided by Odijk for semi-flexible polymers [36]. In an actin gel, the stiffness K could result from the action of motor molecules and therefore be larger. For example, estimating K by the typical stress exerted by myosin mini-filaments on an actin filament in the network suggests values of several tens of pN. The active stress generated by the actin cytoskeleton has been estimated to be up to ζ µ 10 3 Pa [37], which would correspond to ξ 1 µm.
In a Taylor-Couette system with an outer radius R + = 5 mm and an inner radius R − = 1 mm, the critical active stress would be of the order of 10 −6 Pa and thus the system should start to rotate even at very small active stresses. The order of magnitude of the stall torque can be obtained from ∼ |ζ µ|R 2 + , see equation (24) and figure 3. Using again R + 5 mm and ζ µ 10 3 Pa, we estimate a stall torque per unit length 25 × 10 −3 N. To estimate the corresponding rotation rate ω ∼ u r θ , we note that it scales as ω ∼ |ζ µ|/η, see equations (31) and (25) and figure 2. The viscosity η of an actin gel can be estimated as η Eτ , where E is the gel elastic modulus and τ denotes the characteristic visco-elastic relaxation time. The elastic modulus of a passive actin gel scales as E k B T L p λ −4 , where L p is the persistence length of actin filaments and λ is the mesh size of the gel. For k B T 4 pN nm, L p 17 µm and λ 100 nm [38], we estimate that E 10 3 Pa, which is consistent with measurements for actin gels under various conditions [39][40][41]. Using τ 30 s [42], we obtain η 10 4 Pa s. For ζ µ 10 3 Pa, we estimate a rotation rate of about ω 0.1 s −1 .
The Taylor-Couette motor could be realized at small scales with radii of several µm using microfabrication techniques [43]. This leads to a stall torque of the order of several nN and since the rotation rate is rather insensitive to the system size again to rotation rates of the order of 0.1 s −1 . Note that when the system size approaches the mesh site of the acto-myosin gel, the hydrodynamic limit will no longer provide an accurate approximation. Note also that a symmetric system would only rotate if the characteristic length ξ is smaller than the system size. However, even in micrometer-sized systems, we expect many of the general features and the scaling relations to still hold. Thus, the Taylor-Couette motor described here is a promising candidate for the construction of artificial rotational micromotors.