Bifurcation analysis of spontaneous flows in active nematic fluids

Continuum models of active nematic gels have proved successful to describe a number of biological systems consisting of a population of rodlike motile subunits in a fluid environment. One of the prominent features of active systems is their ability to sustain, above a critical threshold of the active parameter, an autonomous collective motion that results in a spontaneous flow of particles. In this paper we show that in a simple channel geometry, the characteristics of this spontaneous motion are largely independent of the model that is used to describe the dynamics of the active system, but are dictated by material symmetry. The natural symmetry for active nematics in a channel is found to be described by the Klein four-group K4≃Z2×Z2 . We perform a Lyapunov–Schmidt reduction and an equivariant bifurcation analysis to show that the K 4-equivariance of the problem generically results in two pitchfork bifurcations with four solution branches.


Introduction
An active system is a collection of self-driven units that are able to continuously exchange energy with the environment and create a collective motion. It is a recent physical paradigm used to represent the collective dynamics of many biological systems such as filaments of the cytoskeleton [15,16], dense bacterial suspensions [20,22,23], and, on a more macroscopic scale, also flocking of birds or fishes [1].
Since these systems are intrinsically out of equilibrium, it is possible to observe many interesting instabilities that lead to very complex patterns of collective motion and self-organized structures [3-5, 9, 10, 14], with features not observed in passive systems, such as internally generated flow patterns, large-scale collective motion, active turbulence, and sustained oscillations.
A distinguishing feature of active system, and the most prominent example of collective motion, is the spontaneous flow of particles induced by activity in a simple shallow channel geometry. Specifically, when the active parameter crosses a suitable threshold the active subunits spontaneously flow along the channel in the absence of any external field.
However, even for this simple example, different models apparently predict different flow profiles. For instance, the analyses in [3,4,15,19] predict only one bifurcation mode at the transition. In particular, when an active gel is confined between two channel walls and the velocity vanishes on both walls, no net flow of particles is observed, but only a shear flow with the formation of bands flowing in opposite directions. By contrast, other analytical and numerical studies [11,12,18,21] on similar equations show that both spontaneous flow in a single direction and banding are possible at the same critical value of the activity. Therefore, the flow-transition seems to have two equivalent modes, a condition which we refer to as 'twofold degeneracy' of the bifurcation.
In the present paper we show that this degeneracy is a generic feature for every K 4 -equivariant system. That is, any active nematic model with symmetric boundary conditions, regardless of the assumptions made to derive it and the specific equations obtained, must have this degeneracy in the bifurcation modes. The bifurcation diagram is essentially determined by the material symmetries of the system. The paper is organised as follows. We present the problem in section 2. For the sake of concreteness we will adopt the model described in [18], but any other set of equations, derived from similar models and possessing the same equivariant properties, could be used. In section 3, we show that the active instability is robust against changes of the model parameters. To have a more precise picture of the bifurcation diagram beyond linear analysis, we perform a Lyapunov-Schmidt reduction of our equations in section 4. In section 5 we find that the natural symmetries of an active nematic system are represented by the Klein group K 4 . By using an equivariant bifurcation analysis, we then show that the obtained bifurcation diagram is a generic feature of any system possessing K 4 invariance. Finally, we draw the conclusions in section 6.

Model description and linear analysis
In order to study the spontaneous flow in a channel geometry, we consider the continuum equations, namely the balance of linear momentum and the director equation, as obtained in [18]. The model described in [18] is based on the classical neo-hookean nematic elastomer free-energy density per unit mass where n is the director or preferred orientation of the molecules, L is the elastomer shape tensor representing the volume-preserving uniaxial stretch along the director, k is the Frank constant, ρµ is the shear modulus, ρ is the density and B e is the effective left-Cauchy-Green deformation tensor. The shape tensor is typically written as where a is the material parameter 3 that gives the amount of spontaneous elongation along n in an uniaxially ordered phase. The shape tensor is spherical, prolate, or oblate, respectively, for a = 1, a > 1 and a < 1.
A key ingredient of this model is that the stress in the material is allowed to relax, so that the effective left-Cauchy-Green deformation tensor B e has its own evolution equation where D is the dissipation tensor, containing information about the relaxation times and viscosity coefficients, B ∇ e := (B e ) . − (∇v) B e − B e (∇v) T is the codeformational derivative of the effective left-Cauchy-Green deformation tensor, and T a is the active tensor. We will take T a = − 1 2 ρµζ I, where ζ is a coefficient that measures the strength of the activity. When the relaxation times contained in D are much shorter than the characteristic times of deformation, the material effectively behaves as a fluid and we obtain a model for an active nematic fluid. We refer to [17,18] for a more detailed derivation of the model equations.
The existence of a spontaneous flow in a shallow channel is a characteristic behaviour of active nematics. We consider an active nematic liquid crystal, confined between two parallel plates at z = 0 and z = L (see figure 1). For simplicity, we consider a system with translational invariance along x, and we assume that the unknown fields are constrained to lie in the (x, z)plane and depend only on the transverse variable z. The nematic director field is determined by the angle θ(z) that n forms with the x-axis. We assume parallel boundary conditions so that θ(0) = θ(L) = 0. The velocity of the fluid, v(z), is along the x-direction (so that the only non-vanishing component is v x (z)) and satisfies the no-slip condition at the channel walls.
Within the active fluid approximation, the governing equations are given in equations (27a) and (27b) of [18]. We write them in dimensionless form and introduce the following nondimensional variables, where τ is the material relaxation time (contained in the definition of D, see [17]) that controls the time evolution of B e via equation (3). The non-dimensional equations of motion read as where the dimensionless ratio r := µL 2 /k = (L/L e ) 2 , compares the length of the channel, L, with an 'elastic length', L e , defined as L e = √ k/µ. In equation (5), there are two unknown fields, V(ξ) and q(ξ), and three dimensionless parameters, namely, a, r and the activity strength ζ. In an abstract setting we can write equation (5) as F(U, ν) = 0 where U = (V, q), ν = (a, r, ζ) and we introduce a smooth mapping F : X × R 3 → Y, between Banach spaces X, Y. We can We observe that F(0, ν) = 0 for any value of the parameters ν and we want to explore the possible bifurcations from the trivial solution. A necessary condition for a bifurcation to occur at ν = ν cr is that the linear operator L, defined as the Frechet derivative L := F U (0, ν cr ), is not invertible. This amounts to studying the nontrivial solutions of the linearized equations It is straightforward to see that the linear problem ( with n ∈ Z − {0}. For simplicity, we only consider the first buckling mode and take n = 1 When ζ = ζ cr , there exist non-trivial solutions to Lu = 0. In such a case, the null-space N = ker(L) is two-dimensional and is generated by the vectors The first of these solutions corresponds to a net flow of particles, while the second represents the banding of particles flowing in opposite directions (see figure 2).  (9). Red arrows represent the velocity field and blue lines identify the director field. Left and right pictures show, respectively, V 1 , q 1 , and V 2 , q 2 .

Robustness against parameter changes
Before analysing the bifurcation in more details, it is interesting to see whether the linear bifurcation analysis, and in particular the mode degeneracy, is robust to model parameter changes. We consider: (1) a different active tensor, (2) one more relaxation time, (3) a different boundary condition for the velocity, and (4) free anchoring conditions for the director 4 . We refer the reader to [13] for the details of the calculations, which, anyway, are very similar to those presented in section 2. For shortness, in the present section we simply report the main results.
More precisely, in case (1) we take T a = − 1 2 ρµζ L, as opposed to T a = − 1 2 ρµζ I which we have used in the previous section; while in case (2) we posit a dissipation tensor of the form D = (L −1 ⊗ L −1 )T, where T now contains two different relaxation times (see [13,17] for details). Introducing more relaxation times entails considering a larger class of possible viscosity coefficients. As a third example, in case (3), we change the no-slip boundary condition at the wall z = L, and impose a free boundary condition v ′ x (L) = 0. Finally, in case (4), we assume free anchoring for the director: q ′ (0) = q ′ (L) = 0.
We find that the qualitative features of the bifurcation diagram do not change in cases (1) and (2), only the value of the critical activity is rescaled (see table 1). Therefore, the linear analysis of the bifurcation and the degeneracy of the modes is unaffected by different active tensors or viscosity coefficients in its fundamental picture. It is worth remarking that these two cases maintain the full symmetry of the original problem.
By contrast, in case (3), the new boundary condition induces a symmetry change in the problem, the symmetry group is now a proper subgroup of the original one. As we shall see, this reduces the dimension of ker(L) and changes the nature of the bifurcation.  (1), (2) and (4). However, in case (4) the banding mode disappears and is replaced by a solution with v = 0 and an arbitrarily uniformly tilted director field. By contrast, the boundary condition change considered in (3) breaks the K 4 -symmetry of the problem and modifies the number of the bifurcation modes.
Case (4) deserves a separate discussion, as we still have two independent solutions, but now ker(L) is generated by u to be compared with (9). Hence, the general solution still comprises a spontaneous flow mode, but the banding mode is replaced by a solution with v = 0 and an arbitrarily uniformly tilted director field. However, as we shall see in section 5, the topology of the bifurcation diagram still coincides with that of cases (1) and (2). The following table 1 summarizes the main conclusions

Lyapunov-Schmidt reduction and bifurcation diagram
In section 3 we have seen that the qualitative features of the bifurcation, and in particular the degeneracy of the modes, are robust to parameter changes. In order to obtain a more precise picture of the bifurcation diagram we perform a nonlinear analysis. A standard method in this context is the Lyapunov-Schmidt reduction [6,7]. This technique allows us to reduce the problem to finite dimensional, and write a finite-dimensional bifurcation equation that is much easier to analyse than the original problem. It is possible to show that L is a Fredholm operator, so that we can split the Banach spaces X and Y into the direct products and we have implicitly used the L 2 -scalar product ·, · to define the orthogonal complements N ⊥ and R ⊥ . Every U ∈ X can be decomposed into the sum U = u + w, with u ∈ N and w ∈ N ⊥ . Likewise, our problem F(U, ν) = 0 can be split into an equivalent pair of equations, where Q is the orthogonal projector onto R and I is the identity. Equation (11a), allows us to derive w as a function of u and ν, while equation (11b), once we substitute w = w(u, ν), is a finite-dimensional equation that yields a complete, albeit local, picture of the bifurcation. Finally, we recall that, in our functional setting, it is often convenient to use the identity R ⊥ = ker(L † ), where L † is the adjoint operator, to calculate R and R ⊥ . In order to solve (11a), we look for solutions that are linear combinations of u 1 and u 2 , as given in (9), plus unknown functions belonging to N ⊥ with To find the projector Q onto R, we use R = ker(L † ) ⊥ . The kernel of L † has a similar structure to N, indeed it is still two-dimensional and its generating vectors, u * 1 and u * 2 , are found to be [13] Furthermore, it is possible to show by direct computation that u * 1 , u * 2 = 0. Hence equation (11a) reads with U as given in (12). This equation allows us to uniquely find ( w v (ξ), w q (ξ) ) for any given values of the amplitudes α, β and the material parameters. However, the explicit calculations can be very cumbersome so that it is usually only possible to find an asymptotic approximation of w when α, β 1. After some algebra with the aid of Mathematica, which we omit for brevity, we find the following form of the ( where O(4) stands for correction terms of at least degree four in α and β. The long expressions of the coefficients are reported in appendix.
With the aid of (16), the operator (I − Q)F(u + w(u, ν), ν) is now a finite-dimensional map N → R ⊥ = ker(L † ), so that it is possible to solve the bifurcation equation (11b) to determine how the mode amplitudes α, β depend on the parameters. Since u * 1 and u * 2 are orthogonal among themselves the bifurcation equation is equivalent to the two conditions After some algebra, we find the following bifurcation equations where λ is such that ζ = ζ cr (1 + λ). We see that α = 0, β = 0 is always a solution, as expected. The non trivial solutions are obtained from (18) by eliminating one of the variables and are described by the following two curves ( The bifurcation diagram (figure 3) shows that two coincident pitchfork bifurcations occur at λ = 0, i.e. ζ = ζ cr , producing four solution branches. Two of these branches (solid red curve), correspond to the spontaneous-flow mode u 1 , while the other two branches (solid blue curve) represent the banding mode u 2 . The nonlinear analysis shows that, although dim(ker L) = 2, the two bifurcation modes do not mix: only solutions with either α = 0 or β = 0 are observed, at least close to the bifurcation. Furthermore, it is clear from equation (19) that the material parameter r do not affect the qualitative features of the bifurcation diagram since it only rescales the amplitudes of the solutions.

Symmetries and equivariant analysis
In sections 3 and 4 we have seen that the nature of the bifurcation is not affected by a change in the model parameters. Furthermore, this result is in agreement with the analyses in [11,21] which were performed on different sets of equations for active nematic fluids. Therefore, it is natural to ask whether the bifurcation is only determined by the symmetry of the system. In this section we give a positive answer to this question. The symmetries of a system of ODEs are specified in terms of a group of transformations of the variables that preserves the structure of the equations. An equation F(U, ν) = 0 is said to be equivariant under the action of a group G if This notion is important because if F is G-equivariant and U is a solution, also gU is a solution for every g ∈ G, and G is a symmetry group for the problem. In our case the system of equation (5) is equivariant under the Klein four-group K 4 ∼ = Z 2 × Z 2 . Specifically, we obtain an equivalent set of equations under a transformation (ξ, V, q) → (ξ,V,q) in the following four cases Hence, the derivatives transform according to 3 : Mirror reflection about the line x = 0. It corresponds to (ξ, V, q) → (ξ, −V, π − q).
Hence, the derivatives transform according to 4 : Rotation of π about the origin. It corresponds to (ξ, V, q) → (1 − ξ, −V, π + q). Hence, the derivatives transform according to . This is the product of g 2 and g 3 .
Correspondingly, we can identify a group action on R 2 that describes how g 1 , g 2 , g 3 , g 4 transform the coordinates in the plane (α, β). The group element g i is then represented by a 2 × 2 matrix T i that can be used to analyse the bifurcation equations (18) (or (19)) and explore the symmetry of the solutions. The matrices T i , i = 1, 2, 3, 4, are reported in table 2.
The bifurcation analysis is carried over the (α, β)-plane. To this end, it is useful to introduce a smooth mapping f : R 2 × R → R 2 : (x, λ) → f(x, λ) so that equation (18) writes as f(x, λ) = 0 where x = (α, β), and we only keep the explicit dependence on the bifurcation parameter λ. The bifurcation occurs at λ = 0.
A key notion in this analysis is that of fixed-point space, associated with a subgroup Σ, and defined as the subspace of all the vectors that are fixed by any g ∈ Σ. In our specific case, we have Because f(·, λ) is K 4 -equivariant, it maps Fix(Σ) to Fix(Σ), and the bifurcation analysis can be restricted to Fix(Σ). The symmetry of a solution x is specified by its stabilizer, Stab(x), consisting of all the transformations that leave x fixed. In our case, for a vector x ∈ R 2 , we can write For our problem, we have the following fixed-point spaces and stabilizers: (i) All points in the (α, β)-plane are fixed by the identity T 1 . Therefore, Fix(g 1 ) = R 2 .
This theorem guarantees the existence of a solution branch for every one-dimensional fixed point space. As such, it characterises the unique smooth solution branch in terms of the axial stabiliser subgroups. Indeed, x ∈ Fix(Σ) and hence the solutions obtained have symmetries including Σ, so Σ ⊂ Stab(w). Since Σ is a stabilizer subgroup, Σ = Stab(w), and the symmetry of the solution x is actually Σ.
Our previous analysis shows that we only have two axial subgroups of K 4 , namely Σ 1 and Σ 3 , so that, by the branching lemma, there must be two independent symmetry-breaking branches, having symmetry Σ 1 and Σ 3 . This result is generic for any system of equations that possesses K 4 -equivariance. By contrast, when we consider a free boundary condition v ′ x (L) = 0, as in section 3, we loose the global K 4 -symmetry since our problem is no longer equivariant with respect to g 2 and g 4 . In this case dim(ker(L)) = 1 so that we observe only one bifurcation mode and no degeneracy.
Furthermore, it is possible to show that generically, for K 4 -symmetry, the only solutions expected at bifurcation are those supplied by axial subgroups (see [7], theorem 2.24, p 42), so that there are no solutions other than those predicted by the equivariant branching lemma.

Conclusions
In the present article, we studied the spontaneous flow arising in active nematics and the robustness of its qualitatively features over material parameter changes. In particular, we noticed that the two-fold mode degeneracy observed in some analytical and numerical studies [11,18,21] does not depend on the particular set of equations used to describe active matter and is largely unaffected by parameter changes.
A Lyapunov-Schmidt reduction and a subsequent equivariant bifurcation analysis allowed us to explore the nature of the bifurcation more carefully. Hence, we find that the bifurcation diagram generically comprises two coincident pitchfork bifurcations lying in orthogonal planes, producing four solution branches that intersect only at the origin. Furthermore, we showed that this bifurcation diagram is generic for any system possessing K 4 Z 2 × Z 2 symmetry. Therefore it should be observed in any theory of active nematic material, regardless of the model assumptions, provided it satisfies the natural K 4 symmetry of these systems. In physical terms, two modes of instability (usually spontaneous flow and banding) should be observed in active nematics at the same critical value of the active parameter. In some cases, the physical realisation of the two bifurcation modes is different, and band instability does not occur. However, regardless of the physical characteristics of the modes, the bifurcation diagram has the same topology, essentially induced by the symmetry of the problem.
Finally, we believe that analysing the active flow through the conceptual framework of symmetry could be a fruitful endeavour also in more complicated cases. For example, it is known that 2D and 3D active flows have different characteristics [2,24], but little is currently known about the behaviour of confined active material in three dimensions. The simplification resulting from explicitly considering the symmetries of the system could make the 3D problem more tractable and possibly help to understand the key differences between 2D and 3D coherent flow.

Data availability statement
No new data were created or analysed in this study.

Acknowledgment
This work was partially carried out under the auspices of the Italian GNFM (Gruppo Nazionale di Fisica Matematica) and is part of the MUR grant activities of "Dipartimento di Eccellenza 2023-2027".