Existence and stability of kayaking orbits for nematic liquid crystals in simple shear flow

We use geometric methods of equivariant dynamical systems to address a long-standing open problem in the theory of nematic liquid crystals, namely a proof of the existence and asymptotic stability of kayaking periodic orbits for which the principal axis of orientation of the molecular field (the director) rotates around the vorticity axis in response to steady shear flow. With a small parameter attached to the symmetric part of the velocity gradient, the problem can be viewed as a symmetry-breaking bifurcation from an orbit of the rotation group~$\SO(3)$ that contains both logrolling (equilibrium) and tumbling (periodic rotation of the director within the plane of shear) regimes as well as a continuum of kayaking orbits. The results turn out to require expansion to second order in the perturbation parameter.


Introduction
Nematic liquid crystals, regarded as fluids in which the high aspect ratio, rigid, rod molecules require descriptive variables for orientation as well as position, are observed to exhibit a wide range of prolonged unsteady dynamical responses to steady shear flow. The mathematical study of these phenomena in principle involves the Navier-Stokes equations for fluid flow coupled with equations representing molecular alignment and nonlocal interactions between rod molecules, typically leading to PDE systems currently intractable to rigorous analysis on a global scale and resolved only through local analysis and/or numerical simulation. It becomes appropriate therefore to deal with simpler models as templates for capturing some of the dynamical regimes of interest and their responses to physical parameters. Stability and bifurcation behaviours that are robust for finite-dimensional dynamical systems, and that numerically reflect the same orbits of interest (specifically, kayaking orbits) in infinite-dimensional systems, provide a framework for extension of rigorous results to the infinite-dimensional systems.
Much of the work on dynamics of liquid crystals (and more generally, rigid large aspect ratio polymers) in fluid flow rests on models proposed by Hess [41] and Doi [19] that consider the evolution of the probability density on the 2-sphere (more accurately, projective space RP 2 ) representing unoriented directions of molecular alignment with the molecules regarded as rigid rods. Extensive theoretical and numerical investigations ( [6], [21], [22], [47], [54], [55], [62], [64], [65], [66] to cite only a few) of these and related nematic director or orientation tensor models in 2D or 3D reveal a wide range of periodic molecular dynamical regimes with evocative names [47] logrolling, tumbling, wagging and kayaking according to the behaviour (steady versus periodic) of the principal axis of molecular orientation (the nematic director) relative to the shear (flow velocity and velocity gradient) plane and vorticity axis (normal to the shear plane). Tumbling orbits, for which the principal axis of molecular orientation rotates periodically in the shear plane, are seen to be stable at low shear rates, but become unstable to out-of-plane perturbations and give way to kayaking orbits, for which the principal molecular axis is transverse to the shear plane, and rotates around the vorticity axis, reminiscent of the motion of the paddles propelling a kayak along the shear flow of a calm stream. The limiting case is logrolling, a stationary state where the principal axis of the rod ensemble collapses onto the vorticity axis, while wagging corresponds to oscillations (but not complete rotations) of the molecular orientation in the shear plane about some mean angle, although wagging regimes do not appear in our analysis. We note very recent experimental results [31] coupled with the high-resolution numerical results of the Doi-Hess kinetic theory [30] that provide overwhelming evidence that the kayaking orbit is responsible for the anomalous shear-thickening response of a high aspect ratio, rodlike, liquid crystal polymer with the acronym PBDT. The papers [27], [31] give extensive lists of literature references.
In the particular case of a steady shear flow and spatially homogeneous liquid crystal in a region in R 3 , the PDEs describing the evolution of orientational order can be simplified to an autonomous ODE in the setting of the widely-used Q-tensor model [17], [59], [72] for nematic liquid crystals. The assumption of spatial homogeneity of course rules out many important applications, to display technology for example, but nevertheless gives a worthwhile approximation in local domains of homogeneity (monodomains) away from boundaries and defects. In this setting the propensity of a molecule to align in any given direction in R 3 is represented by an order tensor Q belonging to the 5-dimensional space V of traceless symmetric 3 × 3 matrices, where { } t and tr denote transpose and trace respectively. The tensor Q is interpreted as the normalised second moment of a more general probability distribution on RP 2 . All such Qtensor models can be associated with a moment-closure approximation of the Smoluchowski equation for the full orientational distribution function [27]. The derivation of the equation yields technical problems concerning the approximation of higher-order moments, a topic of some discussion in the literature: see [23], [27], [46], [48] for example. In this context the dimensionless equation for the evolution of the orientational order takes the general form 1 dQ dt = F(Q, β ) := G(Q) + ω[W, Q ] + β L L L(Q)D (1.2) as an equation in V ∼ = R 5 ; here [W, Q ] = WQ − QW . On the right hand side of (1.2) the first term represents the molecular interactions in the absence of flow, derived for example from a Maier-Saupe interaction potential or Landau -de Gennes free energy: thus G is a frame-indifferent vector field in V . In the second term, W denotes the vorticity tensor, the anti-symmetric part of the (spatial homogeneous) velocity gradient, providing the rotational effect of the flow with constant coefficient ω. In the third term L L L(Q) is a linear transformation V → V applied to the rate-of-strain tensor D, the symmetric part of the velocity gradient, and represents the molecular aligning effect of the flow: the linearity in D is a simplifying assumption. Here L L L(Q) depends (not necessarily linearly) on Q, and L L L(Q)D is frame-indifferent with respect to simultaneous coordinate choice for the flow and the molecular orientation. The coefficients ω and β are constant scalars that depend on the physical characteristics of the liquid crystal molecule as well as the flow. In this study we take ω as fixed, and regard β as a variable parameter.
In the Olmsted-Goldbart model [61] used in [12], [75] the term L L L(Q)D is simply a constant scalar multiple of D. A more detailed model for L L L(Q)D is the basis of a series of studies by the second author and co-workers [26]- [30], [48] as well as by many other authors [8], [35], [56], [63]. We draw attention also to the earlier theoretical work [52,53] assuming a general form for L L L(Q)D and where similar methods to ours are used to study equilibrium states (uniaxial or biaxial), although the question of periodic orbits in general and kayaking orbits in particular is hardly addressed, the existence of the latter having yet to be discovered.
We remark that although in this paper our underlying assumption is of spatial homogeneity there have been studies of nematic liquid crystals dynamics in a nonhomogeneous environment: see among others [13] for analytical results and [77] for numerical simulations.
A particular model of the form (1.2) that 'combines analytic tractability with physical relevance' [60] is the Beris-Edwards model [4], a basis for some more recent investigations [18], [20], [60], [76]  which helps to keep track of the analysis, and also enables the results to apply to simpler models for which one or more of the m i may be zero. For the Beris-Edwards model (1.3) the ratios are (m c : m l : m q ) = (2/3 : 1 : −2), while for the Olmsted-Goldbart model [61] the ratios are (1 : 0 : 0) and for the model in [54] they are ( 3/10 : 3/7 : 0). Moreover, in Appendix B we pursue the analysis for general L L L(Q)D, using the 7-term expression assumed for example in [52,53], and show that with the exception of one term the results are the same as those for (1.5) albeit with different interpretation of the coefficients m c , m l , m q . The exceptional term (being the symmetric traceless form of Q 2 D) also fits into our overall framework as shown in the expressions (B.18) and (B.19) with (B.2).
When β = 0 the equation (1.2) represents the co-rotational case or long time regime, as discussed in [60]. If Q * ∈ V satisfies G(Q * ) = 0 then frame-indifference of G, interpreted as equivariance (covariance) of G under the action of the rotation group SO(3) on V , implies that every element Q of the SO(3) group orbit O of Q * also satisfies G(Q) = 0. If moreover [W, Q * ] = 0 then F(Q * , 0) = 0 and so Q * is an equilibrium for (1.2): the rotational component of the shear flow leaves Q * fixed. This implies that Q * has two equal eigenvalues, and if these are less than the third (principal) eigenvalue then Q * represents a logrolling regime. Moreover, [W, Q ] is tangent to O for every Q = Q * ∈ O and so O (which is topologically a copy of RP 2 ) is an invariant manifold for the flow on V generated by (1.2) when β = 0. The dynamical orbit of every such Q ∈ O is periodic, as it coincides with the group orbit of rotations about the axis orthogonal to the shear plane: in the language of equivariant dynamics [15], [25], [43] it is a relative equilibrium. All of these periodic orbits represent kayaking regimes, except for a unique orbit representing tumbling, and they are neutrally stable with respect to the dynamics on O, as also is the logrolling equilibrium Q * . We discuss this geometry of the SO(3)-action on V in more detail below; it plays a central role in what follows, as it must do in any global study of the system (1.2), an observation of course recognised by other authors [26], [52,53].
There are a few rigorous mathematical proofs of the existence of tumbling limit cycle orbits with limiting assumptions. By positing 2D rods, both with a tensor model [48] and with the stochastic ODE [40], proofs follow from the Poincaré-Bendixson theorem; for 3D rods with a tensor model the proof in [12] uses geometric arguments on in-plane tensors. Until now, there has been no proof of existence of (stable) kayaking orbits, and the purpose of this paper is to provide a proof for second-moment tensor models (1.2), (1.5) at low rates of molecular interaction (although not necessarily low shear shear rates). We thus consider a dynamical regime different from those considered by other authors in numerical simulations such as [66], [27]. A regime analogous to ours in considered in the theoretical work [52,53] using very similar methods, but in that case the molecules are assumed biaxial and it is equilibria rather than periodic orbits that are sought.
The approach we take is to regard β as a small parameter and view (1.2) as a perturbation of the co-rotational case. This enables us to use tools from equivariant bifurcation theory [15], [33], [43], [69,70] and in particular Lyapunov-Schmidt reduction over the group orbit O to obtain criteria for the persistence or otherwise of the periodic orbits of the corotational case after perturbation, and to determine the stability or otherwise of the resulting logrolling, tumbling and kayaking dynamics. Our general results are independent of the choice of the interaction field G, given that it is frame-indifferent and the logrolling state is an equilibrium: G(Q * ) = 0 (Assumptions 1 and 2 in Section 2) and also that the eigenvalues λ , µ of the linearisation of G at Q * normal to O are real and nonzero (Assumption 3 in Section 3). In addition we require a natural condition of frame-indifference for the perturbing field L L L(Q)D (Assumption 4 in Section 3). Finally, the stability results require λ , µ < 0 (Assumption 5 in Section 7). However, our methods do not allow us to make deductions when β is large compared with the rotational coefficient ω. Other limit cycles are possible, and indeed are routinely observed numerically.
Our main result is Theorem 7.7 with Remark 7.9, showing that the existence of a limit cycle kayaking orbit after perturbation depends on the ratio λ /µ as well as the size of the product λ µ relative to the rotation coefficient ω. We show also in Corollary 7.8 that for the Beris-Edwards and Olmsted-Goldbart models the kayaking orbit is linearly stable without further assumption. This paper is organised as follows. In Section 2 we discuss symmetries of the model and key features of the action of SO(3) on V that it inherits from the usual action on R 3 . Of particular importance are the tangent and normal subspaces to the group orbit O. Section 3 gives initial results showing the persistence of log-rolling and tumbling regimes after perturbation, and introduces the rotating coordinate system convenient for further analysis. In Section 4 a natural Poincaré section for the (dynamical) flow near O is described and relevant first-order derivatives of the associated Poincaré map are calculated and shown to vanish. Lyapunov-Schmidt reduction is applied in Section 5 to obtain a real-valued bifurcation function defined on a meridian of O. This function happens to vanish to first order in β and so we are obliged to pursue the β -expansion to second order. In Section 6 we choose L L L(Q)D explicitly as (1.5) and evaluate these second order terms. Finally, in Section 7 the zeros of the bifurcation function are found and the conditions for existence and stability of kayaking motion are determined. For the specific cases of the Beris-Edwards and Olmsted-Goldbart models with Landau-de Gennes free energy the criteria for existence and stability of kayaking orbits are stated explicitly. Following a brief concluding section there are Appendices giving some technical results arising from symmetries that simplify the main calculations, as well as a discussion of how a fully general form of the molecular alignment term L L L(Q)D fits into the framework of our analysis.

Geometry and symmetries of the system
The molecular interaction field G is independent of the coordinate frame and therefore equivariant (covariant) with respect to the action of the rotation group SO(3) on V by conjugation induced from the natural action on R 3 . Therefore our first working assumption in this paper is the following. Further discussion of equivariant maps, in particular relating to the action of SO(3) on V that we shall use extensively in this paper, is given in Appendix A.
Choosing coordinates (x, y, z) ∈ R 3 so that the shear flow velocity field has the form k(y, 0, 0) for constant k = 0 the velocity gradient tensor is k   0 1 0 0 0 0 0 0 0   with symmetric and anti-symmetric parts kD/2 and −kW /2 respectively, where Without loss of generality we take k = 2 since the coefficients ω and β in (1.2) are at present arbitrary. The rotational component W corresponds to infinitesimal rotation about the z-axis.
A nonzero matrix Q ∈ V is called uniaxial if it has two equal eigenvalues less than the third, in which case it is invariant under rotations about the axis determined by the third eigenvalue. Matrices with three distinct eigenvalues are biaxial. In this paper an important role is played by the uniaxial matrix where 0 < a < 1/3 for which the principal axis (largest eigenvalue) is the z-axis and about which Q * is rotationally invariant. We take a > 0 to ensure that Q * is uniaxial, and the upper bound on a is imposed for physical reasons since the second moment of the probability distribution defining the Q-tensor has eigenvalues in the interval [0, 1] and so those of Q are no greater than 2/3: see [3] for example. We exclude a = 1/3 as we shall need to work in a neighbourhood of Q * .
Our second underlying assumption is that this phase is an equilibrium for the system (1.2) in the absence of flow, that is when ω = β = 0. In other words With this assumption, the equivariance property of G implies that G vanishes on the entire SO(3)-orbit O of Q * in V , and O is an invariant manifold for the flow on V generated by (1.2) with β = 0. The dynamical orbits on O coincide with the group orbits of rotation about the z-axis under which Q * remains fixed, this being the only fixed point on O since if Q ∈ O and [W, Q ] = 0 then Q is a scalar multiple of and hence equal to Q * .

Rotation coordinates: the Veronese map
For calculation purposes it is natural and convenient to take coordinates in V geometrically adapted to O. We do this in a standard way by representing the orbit O of Q * as the image of the unit sphere S 2 ⊂ R 3 under the map where again t denotes matrix (or vector) transpose. Here V is the projection to V of the case n = 3 of the more general Veronese map construction R n → R m with m = n 2 and it represents O as a Veronese surface in R 5 : see for example [34] or [39]. It is straightforward to check that V is equivariant with respect to the actions of SO(3) on R 3 and V , that is if where {e 1 , e 2 , e 3 } is the standard basis in R 3 , and that V (e 1 ) and V (e 2 ) are obtained from Q * by permutation of the diagonal terms.
On V we have a standard inner product given by H, K = tr(H t K) = tr(HK). However, the Veronese map is quadratic and does not preserve inner products. Nevertheless, up to a constant factor, its derivative does preserve inner products on tangent vectors to S 2 .
with the dot denoting usual inner product in R 3 , from which it follows that for z ∈ S 2 and u, v ∈ R 3 orthogonal to z Observe that the restriction of V to S 2 is a double cover for all z ∈ R 3 . Through V the familiar latitude and longitude coordinates on S 2 go over to a corresponding coordinate system on O. Any z = e 3 ∈ S 2 can be written using spherical coordinates as for unique θ mod π and φ mod 2π, where R j (ψ) denotes rotation by angle ψ around the jth axis in R 3 , j = 1, 2, 3, so that in particular Hence by (2.6) and equivariance (2.3) any Z ∈ O can be written (not uniquely) as for some z ∈ S 2 , as the counterpart of (2.6) using rotations R on V in place of R on R 3 . We shall make frequent use of this notation throughout the paper.
By analogy with S 2 we call each closed curve θ = const = 0 mod π on O a latitude curve and each curve φ = const on O a meridian. It follows from (2.5) that all latitude curves are orthogonal to all meridians. The case θ = 0 mod π corresponds to Q * , and so we think of Q * as the north pole of O ∼ = RP 2 .

Isotypic decomposition
The rotation symmetry of O about the north pole Q * plays a fundamental role in our analysis of (1.2) for sufficiently small nonzero β , and enables us to choose coordinates in V that are strongly adapted to the inherent geometry of the problem. More generally, for any z ∈ S 2 let denote the isotropy subgroup of z (namely the group of rotations about the z-axis) under the natural action of SO(3) on R 3 . Equivariance of V implies that Σ z also fixes Z = V (z) in O under the conjugacy action, and moreover Z is an isolated fixed point of Σ z on O since z is an isolated fixed point of Σ z on S 2 .
At this point it is convenient to develop some further machinery from the theory of linear group actions to describe key features of the geometry highly relevant to our analysis. Introductions to the theory of group actions and orbit structures can be found for example in [1], [14], [58]. We shall make much use of the further fact that corresponding to the action of Σ z on V there is an isotypic decomposition of V (for theoretical background to this notion see for example [15], [25], [33]) into the direct sum of three Σ z -invariant subspaces on each of which Σ z acts differently: the element R z (ψ) ∈ Σ z denoting rotation about the z-direction through angle ψ acts on V Z k by rotation through kψ for k = 0, 1, 2. In particular, with z = e 3 and Z = Q * writing where the mutually orthogonal matrices E 0 , E 1 (α), E 2 (α) are given by and we set Here R 3 (φ ) acts on V * 1 and V * 2 by where we keep in mind that E 2 (α) is defined in terms of 2α. For Z = Z(θ , φ ) as in (2.7) we use the notation and A consequence of SO(3)-equivariance is that for Z ∈ O the derivative DG(Z) : V → V respects the decomposition (2.8) and commutes with the Σ z -rotations on each component. A further important consequence that simplifies several later calculations is the following.
The only linear map V → R invariant under all rotations of V Z 1 and of V Z 2 must be zero on those components. ⊓ ⊔

Alignment relative to the flow
Since the element R 3 (π) ∈ SO(3) acts on V * k by a rotation through kπ it follows that V * 0 ⊕V * 2 is precisely the fixed-point space for the action of R 3 (π) on V . Thus Q = (q i j ) ∈ V is fixed by R 3 (π) if and only if q 13 = q 23 = 0, in which case q 33 is an eigenvalue with eigenspace the z-axis and the other eigenspaces lie in (or coincide with) the x, y-plane. It is immediate to check that if Q = pE 0 + qE 2 (α) then the eigenvalues of Q are 2p/ √ 6 and (−p ± √ 3q)/ √ 6 and so Q has two equal eigenvalues precisely when In the first case Q = pE 0 , while in the second case the eigenvalues are 2p/ √ 6 (repeated) and −4p/ √ 6 so that if p < 0 then Q is uniaxial with principal axis lying in the x, y-plane. From the point of view of the liquid crystal orientation relative to the shear flow such matrices Q are called in-plane; nonzero matrices which are not in-plane are called outof-plane. This agrees with standard terminology where tumbling and wagging dynamical regimes are described as in-plane (see [21], [64] for example), while logrolling and kayaking are out-of-plane.
Let C denote the equator {θ = π/2} of S 2 , and let C = V (C) ⊂ O which we also call the equator of O. It is straightforward to check that With z = (cos φ sinθ , sin φ sin θ , cos θ ) t in usual spherical coordinates we find z · E 1 (α)z = (1/ √ 2) sin2θ cos(φ − α) which vanishes for all α just when sin 2θ = 0, that is θ = 0 or θ = π/2 corresponding to Z = Q * or Z ∈ C respectively. ⊓ ⊔ since G(Q) = 0 for Q ∈ O, giving solution curves t → R 3 (ωt)Q each of which has least period 2π/ω apart from the equilibrium Q * and the equator C : this has least period π/ω, the equator C of S 1 being a double cover of C via the Veronese map. A matrix Q ∈ C is inplane and its dynamical orbit corresponds to steady rotation of period π/ω about the origin in the shear plane, and so C represents a tumbling orbit. All latitude curves of O other than the equator C represent kayaking orbits of period T 0 = 2π/ω and of neutral stability on O and so most of them are unlikely to persist for β = 0. The geometry can be visualised as follows: removing the poles at z = ±e 3 from S 2 leaves an (open) annulus foliated by circles of latitude, so that removing Q * from O leaves a Möbius strip foliated by closed latitude curves each of which traverses the strip twice since Z(π/2 + θ , φ ) = Z(π/2 − θ , φ + π), except for the 'central curve' C given by θ = 0 which traverses it only once.

Tangent and normal vectors to the group orbit O
However, for Z = Q * the tangent space T Z is also spanned by the tangents at Z to the meridian and latitude curve of O through Z. Proof If φ = 0 the vectors R 2 (θ ) e 1 and e 2 = R 2 (θ ) e 2 are respectively tangent to the meridian and latitude of S 2 through z ∈ S 2 , and so applying R 3 (φ ) gives that the vectors R z e 1 and R z e 2 are respectively tangent to the meridian and latitude through z in the general case. Therefore the corresponding tangent spaces at and so applying R z gives the result. ⊓ ⊔ which we shall make use of below.
From Corollary 2.5 it follows that the normal space N Z to O at Z (the orthogonal complement in V to the tangent space T Z ) is given by (2.23)

The dynamical system after perturbation
Since G is SO(3)-equivariant and so in particular is equivariant with respect to the action of the isotropy subgroup Σ z on V , the fact that Σ z fixes Z means that the derivative DG(Z) : V → V respects the decomposition (2.8). Moreover, Assumption 2 and equivariance imply that G vanishes on the entire orbit O and so DG(Z) vanishes on T Z = V Z 1 . Let λ denote the eigenvalue of DG(Z) on V Z 0 = span{Z}, which by equivariance is independent of Z ∈ O. Since DG(Z) commutes with the rotation action of Σ z on V Z 2 its two eigenvalues on V Z 2 are complex conjugates and again independent of Z; we assume them to be real (as they will be in the gradient case, of most interest to us) and denote them by µ (repeated).
Even without the assumption µ ∈ R but with λ and ℜ(µ) both nonzero the manifold O is normally hyperbolic and therefore it persists as a unique nearby smooth flow-invariant manifold O(β ) for (1.2) for sufficiently small |β | > 0: see [24], [42] for the general theory invoked here. Our interest is to discover which periodic orbits on O persist as periodic orbits after such a perturbation.

Remark 3.1
The same approach is used in [52,53] to detect steady states (equilibria) bifurcating from more general group orbits. The geometry of the tangent and normal spaces to all orbits of SO(3) in V is exploited there in a significant way, although using constructions slightly different from ours.
We now make explicit the assumption of linearity and frame-indifference of the contribution to (1.2) from the non-rotational component of the shear flow. The frame-indifference is natural for a physical model, while the linearity is generally assumed for simplicity: see for example [51] and compare equation (4) in [36]. It is immediate to check that Assumption 4 holds for (1.5). As a consequence we have the following elementary result. Proof Using (2.23) we see from Section 2.3 that N * is the fixed-point subspace for the From the symmetry and Corollary 3.3 we have two immediate results: the north pole Q * equilibrium (logrolling) and the equator C periodic orbit (tumbling) persist after perturbation. (ii) a smooth family of periodic orbits C (β ) in N * with C (0) = C and with period tending to π/ω as β → 0.
(ii) The equator C lies in O ∩ N * and is an isolated periodic orbit in N * with characteristic multipliers there e πλ /ω and e π µ/ω (repeated). We seek a fixed point for the first-return map on a local Poincaré section. Since the multipliers differ from 1, the Implicit Function Theorem applied on N * gives the result. ⊓ ⊔

Rotated coordinates
The effect of the perturbation β L L L(Q)D on the system (1.2) when β = 0 is most usefully understood in terms of a co-moving coordinate frame that rotates with the unperturbed system (β = 0), since in these coordinates the rotation term [W, Q ] vanishes (cf. [60, Section 2]). Explicitly, with W = W 3 and the substitution using Assumption 4 on frame-indifference of L L L(Q). Thus in (3.1) and with Q R again written as Q the rotation term [W, Q ] has been removed from (1.2) at a cost of replacing D by the time-dependent term R 3 (−ωt)D.
For given β we denote the flow of (1.2) by ϕ t (·, β ) : V → V , and denote the time evolution map of the nonautonomous system (3.1) by To simplify notation in what follows we choose t 0 = 0 and write for Q ∈ Ṽ Observe in particular that for T 0 = 2π/ω

Local linearisation: the fundamental matrix
An important role will be played by the linear transformation (fundamental matrix) that satisfies the local linearisation of (1.2) (also called the variational equation [49,Ch.VIII], [37, p.23]) along theφ-orbit of Q when β = 0 , namelẏ with respect to the same decomposition (2.8). In particular we have the following key fact.
In what follows we shall make much use of this result, which states that the tangent space T Z = V Z 1 to O at Z consists of equilibria of the variational equation at Z.

The Poincaré map
All points Z ∈ O satisfy ϕ T 0 (Z, 0) = Z for T 0 = 2π/ω. Our aim is to discover which of these periodic orbits persist for sufficiently small |β | > 0, and to discern their stability. Systems of the form (3.1) (not necessarily with symmetry) have a long pedigree in the differential equations literature; in our application the symmetry plays a crucial role. The method we use is to apply Lyapunov-Schmidt reduction to a Poincaré map to obtain a 1-dimensional bifurcation function, and to look for its simple zeros when β = 0: by standard arguments as in [5], [7], [10], [16], [32] for example, these correspond to persistent periodic orbits. The existence of zeros Z for small |β | is established by taking a series expansion of the bifurcation function in terms of β with coefficients functions of Z. Expressions for these coefficients in a general setting are given in [7], and in principle we could simply set out to evaluate these expressions in our case. However, in so doing we could lose sight of important geometric features of V that are fundamental to the shear flow problem, and therefore instead we re-derive the relevant terms explicitly in our symmetric setting. (2.7) with Z = Q * . Let B * denote the orthonormal basis for V :

Poincaré section
where E 0 and E i j for i, j ∈ {1, 2} are defined in (2.12) and (2.13). Let B Z denote the rotated basis (also orthonormal) with notation as in (2.16).
so that V = T Z ⊕ N Z by Corollary 2.5, and so for sufficiently small ε 0 > 0 the union To construct a Poincaré section for the flow of (1.2) we restrict Z to lie on a chosen meridian is a smooth 4-manifold that intersects O transversely along M . Moreover, F(Z, 0) is nonzero and orthogonal to U M is a global (along M ) Poincaré section for all the (periodic) orbits through M \Q * generated by the unperturbed vector field F(·, 0). The least period for Z ∈ M \Q * is T 0 = 2π/ω, with the exception that if Z lies on the equator C then the least period is T 0 /2 = π/ω. We next show that there exists 0 < ε ≤ ε 0 such that the corresponding U ε M is in an appropriate sense a Poincaré section for all orbits close to O generated by the perturbed vector field F(·, β ) including those lying in Q * + N * ε .
M a tubular neighbourhood of O restricted to M constructed using the normal bundle as in (4.4). Then there exists β 0 > 0 and 0 < ε ≤ ε 0 and a smooth function A key part of Proposition 4.1 is the smoothness of T on all of its domain including Q * + N * ε × (−β 0 , β 0 ), since there F(·, β ) lies in N * by Corollary 3.3 and so T is not strictly a 'time of second return'.
for j = 2, 3. We show that there is a positive smooth function on a neighbourhood of M that coincides withφ away from Q * + N * ε , so that time t can in effect be replaced there by angle φ .
Inspecting the terms on the right hand side of (4.5) we find from (4.6) and using (4.7) Next, we take the inner product of (4.5) with orthogonal to the right hand side of (4.8) and toU ∈ V * 2 : since R z preserves inner products this annihilates theθ andU terms in (4.5) and leaves The next step is to replaceQ in (4.10) by the right hand side F(Q, β ) of (1.2). Since G respects the isotypic decomposition (2.8) we have by equivariance G(Z +U Z ) ∈ N Z and so as in (4.9), (4.10). Finally, (4.13), and so we focus on D 0 T . We have by elementary matrix evaluation Therefore from (4.13) and (4.16) where Substituting (4.11), (4.12) and (4.17) into (4.10 Taking ε > 0 small enough so that b(a, u) > 0 and dividing (4.19) through by b(a, u) sin θ , for β sufficiently small we haveφ > ω/2 for θ = 0 and we observe thatφ extends smoothly to θ = 0, corresponding to Q ∈ N * .
Consequently in (θ , φ ,U) coordinates for sufficiently small |β | and |u| the flow has positive component in the φ direction. Since O (given by u = 0) in invariant under the flow of (1.2) when β = 0 and is given by φ -rotation only, it follows that for ε and |β | sufficiently small and Q = Z + U Z ∈ U ε we can define T (Q, β ) to be the time-lapse from φ = φ 0 to φ = φ 0 + 2π if Z = Q * and to be T 0 = 2π/ω when Z = Q * . ⊓ ⊔ Now we are able to define a Poincaré map close to O and for sufficiently small |β |.
where T (Q, β ) is as defined in Proposition 4.1.
The bifurcation analysis that follows proceeds by expanding P(Q, β ) in terms of the perturbation parameter β .

First order β -derivatives
Differentiating (4.20) with respect to β gives where here and throughout we use ′ to denote differentiation with respect to the second component β . At Recall from Lemma 2.4 that the tangent space to M at Z is spanned by E Z 11 (0) while the tangent space to the latitude curve through Z is spanned by E Z 12 (π/2). It follows that the derivative Dm(Z) : V → R annihilates E Z 11 (0) and is an isomorphism from span{E Z 12 (π/2)} to R, so that in particular Dm(Z)F(Z, 0) = 0 (4.25) since from (2.22) we see We next introduce a variable that plays a central role in subsequent calculations.

Lyapunov-Schmidt reduction
Our aim in this section is to seek solutions Q = Q(β ) ∈ U ε M for sufficiently small |β | > 0 to the equation M is as in (4.20), and to determine the stability of the T (Q(β ), β )periodic orbit of (1.2) that each of these represents. Of particular interest are out-of-plane solutions, corresponding to kayaking orbits. We apply Lyapunov-Schmidt reduction to (5.1) along M exploiting the SO(3)-invariant tangent and normal structure to O.
Lyapunov-Schmidt reduction is a fundamental tool in bifurcation theory, and amounts to a simple application of the Implicit Function Theorem. Accounts can be found in many texts such as [ [11], [38], [57]. Although the method is local in origin, it can be applied globally on a manifold on which a given vector field vanishes, or on which given mapping is the identity, by piecing together local constructions and invoking the uniqueness clause of the Implicit Function Theorem. This is the version we use here, which fits into the general framework of [5,7,50] and has significant overlap with the geometric methods of [52,53].
Let Q ∈ N Z . Then Q N = Q and Q T = 0 where the suffices N, T will denote projections to N Z , T Z respectively. Hence (5.1) is equivalent to the pair of equations When β = 0 the equation (5.2) is satisfied by Q = Z, and by (3.9) the Q-derivative is an isomorphism. It follows by the Implicit Function Theorem and the (smooth) local triviality of the normal bundle, as well as the compactness of M , that for all sufficiently small |β | there exists a smooth section of the normal bundle of O restricted to M such that for sufficiently small |β | the map has the property that It therefore remains to solve the equation (5.3) along M given (5.4), that is to solve the reduced equation or bifurcation equation for (Z, β ) ∈ M × R and for |β | sufficiently small. Since T Z = V Z 1 = span{E Z 11 , E Z 12 } and by construction the Poincaré map P has no component in the direction of the vector field E 12 , the bifurcation equation (5.5) can by Lemma 2.4 be written more specifically as with P 11 = p Z 11 P where p Z 11 denotes projection to span{E Z 11 }. We thus seek the zeros of the bifurcation function F (·, β ) : M → R where for sufficiently small |β | > 0. We shall find these by taking a perturbation expansion of F (Z, β ) in terms of β .

Perturbation expansion of the bifurcation function
First, we need a β -expansion of the Poincaré map P which we write as We also make use of the 'approximate' Poincaré map for i = 0, 1, and P 2 (Z) −P 2 (Z) ∈ span{F(Z, 0)}.

First order term of the bifurcation function
Here we denote as in (5.7), and likewise for the second derivatives. . Thus both terms on the right hand side of (5.16) vanish.

⊓ ⊔
A geometric interpretation of Proposition 5.3 is that to first order in β the SO(3)-orbit O, on which every dynamical orbit (other than the fixed point Q * ) is 2π/ω periodic, perturbs to an invariant manifold with the same dynamical property, so that neutral stability of all periodic orbits is preserved.

Second order term of the bifurcation function
Given that the first order term in the β -expansion of F (Z, β ) vanishes by Proposition 5.3 we turn to the second order term. Differentiating P 11 (Z, β ) twice with respect to β at β = 0 we obtain from the left hand side of (5.6) where we write P i 11 for p Z 11 P i , i = 0, 1, 2. To evaluate (5.18) a significant simplification can be made.
Proposition 5.5 P may be replaced byP in all terms on the right hand side of (5.18).
Proof By Proposition 5.1 and Proposition 5.2 each term differs from its counterpart withP by a scalar multiple of F(Z, 0), which is annihilated by p Z 11 . ⊓ ⊔ We next investigate in turn each of the terms of (5.18) withP in place of P.

Second Q-derivative ofP
so that in particularφ t 0 (Q) =φ t (Q, 0), we see from (3.1) with β = 0 that D 2φ t 0 (Q) satisfies the equation and so we obtain from the variation of constants formula and (3.4) Since σ ′ (Z, 0) ∈ N Z and so by (3.9) also M(s, using Corollary 3.5 and the bilinear property of D 2 G(Z) given in Proposition A.7.1.

First Q-derivative ofP 1
By definition of the solution to (3.1) through Q at t = 0 we havė Differentiating with respect to β at β = 0 we obtaiṅ To evaluate the term involving DP 1 in (5.18) we must next substitute H = σ ′ (Z, 0) ∈ N Z into (5.27). We write to emphasise the tangent and normal character of these projections.
where y T (t, Z) := p Z T y(t, Z) and B Z 11 := p Z 11 B Z . Finally, to complete the evaluation of (5.18) we make explicit the term involvingP 2 (Z) in that equation.

Explicit calculation of the bifurcation function
For explicit calculation of the second order term F ′′ (Z, 0) we now take Z = Z(θ , φ ) and express the bifurcation function (5.7) in terms of θ and φ . The choice of φ is arbitrary so we expect the existence and stability results for periodic orbits to be independent of φ , but nevertheless we retain φ at this stage as a check on the calculations.
Up to this point our analysis has assumed little more than the SO(3)-equivariance (that is, frame-indifference) of the vector field G and the perturbation term L L L(Q)D in the system (1.2) and the fact that Q * is an equilibrium for the unperturbed (β = 0) system. To proceed further and evaluate F 2 (Z) we now need to make an explicit choice for the form of L L L(Q)D.

Choices for the perturbing field L L L(Q)D
We consider in turn the three terms comprising the field L L L(Q)D in (1.5), that is where D as in (2.1) represents the symmetric part of the shear velocity gradient and we recall the notation (1.4). From (2.14) we obtain: Lemma 6.1 In the co-moving coordinate frame as in Section 3.1 the perturbation terms become respectively H. ⊓ ⊔ We next need expressions for the components of E 2 (π/4 − ωt) in the basis B Z as in (4.2). These could formally be found using 5 × 5 Wigner rotation matrices describing the action of SO(3) on V as in physics texts such as [68], but in our case it will be simpler to calculate directly.

Calculation of y(t, Z)
Armed with these coefficients we are now in a position to calculate y(t, Z) and subsequently χ(t, Z), needed in order to evaluate (5.36). We consider in turn the three cases (i),(ii) and (iii) of Section 6.1, denoting the corresponding y by y c , y l , y q respectively.
From (4.30) and using (3.9) we have For convenience we now introduce the polar coordinate notation (ν, 2ω) = r ν (cos 2γ ν , sin2γ ν ) (6.12) for ν = λ , µ, as well as the abbreviations r ν e νt cos(2φ + 2γ ν ) − cos(2ωt + 2φ + 2γ ν ) (6.14) with the limiting cases The cases when t = T 0 = 2π/ω will also be important: Using these we obtain from Corollary 6.4 the following expression for y c (t, Z) in terms of the basis B Z . where y c 01 = √ 2c 01 (θ )S(t, φ , λ ) Proof By Corollary 6.4 and (6.11) and The calculations for y c 12 , y c 21 and y c 22 are very similar. ⊓ ⊔ since Proposition A.2 shows that L L L l (t, Z)D differs from L L L c (t, Z)D only in that the coefficients multiplied by 2a, a, −2a respectively. Hence in this case the result is the following.
Since χ is linear in y (see (6.28)) we immediately deduce from Proposition 6.6 the relations Again since χ is linear in y it follows from (6.26) that Z), 0, 0). (6.34)

The bifurcation function
We are now ready to calculate the terms appearing in the expression (5.36) that determine the bifurcation function. With y T = y 1 and χ N = χ 0 + χ 2 the first term is using the bilinearity of B Z . Substituting χ c 0 (t, Z) from (6.29) and using Proposition 6.5 and Corollary A.11 we find (6.38) since by (6.15) T 0 0 sin(4ωt + 4φ + 2γ λ ) + sin2γ λ dt and sin2γ λ = 2ω/r λ as in (6.12). Here we have introduced the notation Likewise from (6.31) with Proposition 6.5 and Corollary A.11 we have Observe that as anticipated the expressions (6.38) and (6.39) do not depend on the meridianal angle φ .
We now turn to the second term appearing in the expression (5.36) for the bifurcation function, namely where χ 01 = χ 01 (t, Z) etc. denote the coefficients of χ(t, Z) in the basis B Z . We now take χ = χ c and evaluate (6.40) by integrating (6.41) from t = 0 to t = T 0 . Straightforward trigonometrical integrals using (6.29) and Corollary 6.4 give since λ = r λ cos 2γ λ . Moreover, by (6.32) and Corollary 6.4 and also while from (6.31) Similarly we find Thus (6.41) and (6.42)-(6.47) give Using the above calculations, we can now evaluate the bifurcation function for the term L L L(Q)D as a linear combination (1.5) of Cases (i),(ii),(iii). The corresponding y term has the form y(t, Q) = m c y c (t, Q) + m l y l (t, Q) + m q y q (t, Q) (6.49) with y c ,y l and y q as in Proposition 6.5 with (6.21) and (6.23),(6.25) respectively, while with the relevant components given by (6.29), (6.31), (6.32) for χ c , by (6.33) for χ l and by (6.34) for χ q .
First take the restricted case m q = 0. Here we find from (5.36) since D L L L c (t, Z) = 0. Writing the bifurcation function F (Z, β ) in coordinates as here dropping the redundant variable φ , we observe from (6.38) and (6.39) (for B Z ) and (6.48) (for D L L L) that each term in f 2 (θ ) is a linear combination of the two terms: the notation reflecting the fact that it is only the components χ 0 and χ 2 of χ that play any role here.
The coefficients of these arising from the various terms that appear in (6.51) are respectively as follows: and so collecting up terms in (6.51) gives Now consider terms involving m q , not yet included. Since y q 1 = 0 from (6.26) and χ q N = χ q 0 from (6.34) the only terms that arise from B Z are from (6.38), and likewise from (6.21) Regarding terms arising from D L L L, we have from (6.41) and (6.34) (6.60) and so using (6.42). Observe also that from Proposition 6.2 (iii) and (6.5), (6.32) as in (6.43), and so likewise Therefore the contribution of the m q -term in L L L(Q)D is merely to replaceΛ i in (6.56) by and Λ 2 =Λ 2 , the m l m q terms from (6.59) and (6.61) fortuitously cancelling.

Corollary 7.3 For the Beris-Edwards model and the Olmsted-Goldbart model the second order term f 2 (θ ) of the bifurcation function has no zeros
These models both have m c = 0. If m c = 0 with m l = 0 (so L L L(Q)D has linear but no constant term) then Λ 0 = 0 while Λ 2 = 0 and we see from (7.7) that f 2 (θ ) does not vanish for any θ = 0, π/2 mod π.

Periodic orbits
Since as in (6.52), the Implicit Function Theorem implies that if θ = θ 0 is a simple zero of f 2 then for sufficiently small |β | > 0 there exists a unique θ β close to θ 0 such that the right hand side of (7.9) vanishes at θ = θ β and θ β → θ 0 as β → 0. Thus θ β corresponds to a solution In fact we know by Proposition 3.4 that the solutions θ = 0, π/2 corresponding to the north pole Q * and equator C do persist for sufficiently small |β | > 0, and we verify that and so the north pole solution is always a simple solution, while the equator solution is a simple solution provided Λ 0 τ λ = Λ 2 τ µ . In general, if θ = θ 0 := π/2 ±Θ is another zero of which is nonzero since Λ 0 τ λ ,Λ 2 τ µ have the same sign by Corollary 7.2, and so θ 0 is also a simple solution.
When Λ 0 = Λ 2 as in the Beris-Edwards or Olmsted-Goldbart models we thus have the following result on periodic orbits after perturbation.  for fixed λ , µ with τ λ /τ µ < 1 the equator C is the unique periodic orbit on O (other than the equilibrium Q * ) that persists for sufficiently small |β | > 0; its period is close to π/ω. For τ λ /τ µ > 1 there is in addition β 0 > 0 and a smooth path {Q(β ) : |β | < β 0 } in V with Q(0) = Z(θ , φ ) ∈ M φ where θ = π/2 ±Θ as in Corollary 7.3 such that there is a periodic orbit of (1.2) through Q(β ) with period T (Q(β ), β ) → T 0 = 2π/ω as β → 0. ⊓ ⊔ The perturbed equator represents a periodic orbit close to tumbling, possibly with a small kayaking and/or biaxial component. The periodic orbit through Q(β ) represents a kayaking orbit that (for fixed λ , µ) arises from a particular kayaking orbit on O persisting after perturbation. The two values θ = π/2±Θ correspond to the two intersections of the same periodic orbit with the Poincaré section: see the geometric description at the end of Section 2.3. Thus if (sufficiently small) β = 0 is fixed and τ λ /τ µ increases through 1, the equator tumbling orbit generates a kayaking orbit through a period-doubling bifurcation.

Stability
So far the discussion has rested on Assumption 3 ensuring the normal hyperbolicity of the SO(3)-orbit O under the dynamics of the system (1.2) when β = 0. In this section we investigate dynamical stability of the periodic orbits on O that persist close to O for sufficiently small |β | > 0. A necessary condition for stability is that O itself be an attracting set, and so we make now the following further assumption: Assumption 5: The eigenvalues λ , µ of DG(Q * ) are negative.
Consequently the perturbed flow-invariant manifold O(β ) is normally hyperbolic and attracting for sufficiently small |β | > 0, therefore the stability of any equilibrium or periodic orbit lying on O(β ) is determined by its stability or otherwise relative to the system (1.2) restricted to O(β ). The manifold O(β ) can be seen as the image of a section of the normal bundle of O, its intersection with U ε M being the image M (β ) of a sectionσ (·, β ) of this normal bundle restricted to M = M φ . The 1-manifold M (β ) is invariant under the Poincaré map P(·, β ), the restriction of P(·, β ) to M (β ) determining a 1-dimensional discrete dynamical system on M (β ) whose fixed points correspond to periodic orbits (or fixed points) of F(·, β ) on O(β ).
In our analysis, rather than useσ (·, β ) which is harder to compute, we have used σ (·, β ) and the method of Lyapunov-Schmidt to construct a vector field on M whose zeros correspond to the periodic orbits (or fixed points) of F(·, β ) on O(β ). It follows from the general Principle of Reduced Stability [45], [74] that stability of periodic orbits on O(β ) corresponds to stability of the corresponding zeros of the vector field F (Z, β )E Z 11 on M in the present context where dim(M ) = 1. However, we now show this directly, using a simple geometric argument taken from [16,Section 9.4]. Recall that in terms of the θ -coordinate for Z on M we have F (Z, β ) = f (θ , β ). Proof Suppose this fails for a given fixed value of β , so that (without loss of generality) Q(β ) is attracting on M (β ) while θ 0 is repelling on M . In particular this means that there is an interval (θ − , θ 0 ) such that all corresponding points on M (β ) are moved to the right (greater θ -value) by the Poincaré map P(·, β ), and there is also an interval (θ 0 , θ + ) on which f (θ , β ) > 0. Now consider a perturbation of the system (1.2) which adds a vector field of the form Q → ζ (Q)E Z 11 where ζ : V → R is a smooth non-negative bump function with ζ (Q(β )) > 0 and vanishing outside a sufficiently small neighbourhood U of Q(β ) in V . Note that such a perturbation will be far from SO (3)

Stable kayaking orbits
We are now able to describe the global dynamics close to O for the Beris-Edwards model, under the standing Assumptions 1-5. From Corollary 7.4 and (7.10),(7.11) we deduce the following stability result.
Theorem 7.7 For the Beris-Edwards model first suppose Λ 2 > 0. Then for τ λ /τ µ < 1 and sufficiently small |β | > 0 the perturbed equator C (β ) is an attracting limit cycle (close to tumbling) on the invariant manifold O(β ) that is the perturbed SO(3)-orbit O, its basin of attraction on O(β ) being the whole of O(β ) apart from the perturbed equilibrium Q * (β ) (log-rolling). For τ λ /τ µ > 1 the perturbed equator C (β ) is a repelling limit cycle, and there is precisely one other limit cycle on O(β ): this limit cycle (kayaking) is attracting, and has period approximately twice that of C (β ). If Λ 2 < 0 the attraction/repulsion is reversed. ⊓ ⊔ For the simpler Olmsted-Goldbart model we have Λ 2 = m 2 c > 0 and so stability of the kayaking orbit (when it exists) automatically holds. In general we have and so the stability condition Λ 2 > 0 holds precisely when w < −4a or w > 2a where w := m c /m l , supposing m l = 0. If m l = 0, m c = 0 the kayaking orbit is automatically stable if it exists, while if m l = 0, m c = 0 there is no kayaking orbit. For the Beris-Edwards model we have w = 2/3, and in this case stability depends on the coefficient a and holds automatically given that a < 1/3. Thus, to summarise: Corollary 7.8 For the Beris-Edwards and Olmsted-Goldbart models, if the SO(3)-orbit of the logrolling equilibrium Q * is normally hyperbolic and attracting (so that Q * is a stable equilibrium state in the absence of the shear flow, up to rigid rotations) then the kayaking orbit, when it exists, is an asymptotically stable limit cycle. Remark 7.9 Given Assumption 5 the condition τ λ /τ µ > 1 is the same as τ λ < τ µ , that is where f : R m → R is a smooth function and {X 1 , X 2 , . . ., X m } are a basis for the ring of SO(3)-invariant polynomials on V . It is well known in the liquid crystal literature (see for example [52, eq.(4.9)]) that such a basis is given by {X,Y } where a proof being given in [33, Ch.XV, §6] via reduction to the group of symmetries of an equilateral triangle. Note that for Q ∈ V the Cayley-Hamilton Theorem shows immediately that tr Q 3 = 3 det Q.
With f X , f Y denoting the partial derivatives of f we find that the functions g,ḡ of (A.9) are then given by and so also for their derivatives 14) The equilibrium condition (A.10) is respectively. The eigenvalues of DG(Q * ) are λ , µ and 0 where by (A.14) and (A.15), with ∆ f * X := D f X (Q * )Q * and likewise ∆ f * Y . For the particular and important case of the Landau -de Gennes potential and The equilibrium condition (7.15) is thus that the coefficient a > 0 should satisfy τ + 6a 2 c − ab = 0 (7.23) and the eigenvalues λ , µ are given by Here µ is automatically negative, and it is straightforward to check that (7.23) has two real solutions 0 < a 1 < a 2 provided 0 < τ < b 2 /(24c). Then a 2 > 1 2 (a 1 + a 2 ) = b/(12c) and so λ < 0 for a = a 2 and we choose a = a 2 in the definition of Q * .

Corollary 7.11
In this setting the result of Theorem 7.10 giving the condition for the existence of kayaking orbits becomes with stability for a = a 2 < 1/3. ⊓ ⊔ It is natural to ask for what range of values of b, c, τ and ω these conditions can simultaneously hold.
Proposition 7.12 A necessary condition for the existence of stable kayaking orbits is b < 4c.

Conclusion
The geometry of uniaxial and biaxial nematic liquid crystal phases is most naturally expressed in terms of the action of the rotation group SO(3) on the 5-dimensional space V of (symmetric, traceless) Q-tensors. In this paper we have used techniques from bifurcation theory related to symmetry, applied to a rather general class of ODEs on V widely used to model a homogeneous nematic liquid crystal in a simple shear flow, in order to prove the existence under certain conditions of an asymptotically stable limit cycle representing a 'kayaking' orbit, where the principal axis of molecular orientation of the ensemble of rigid rods lies out of the shear plane and rotates periodically about the vorticity axis. Our key assumption, however, is that the dynamical effect of the symmetric part of the flow-gradient tensor should be small compared to that of the anti-symmetric (rotational) part, so that the system we study is viewed as a perturbation of the co-rotational case which involves only the (frame-indifferent) molecular interaction field in addition to the rotation of the fluid. The results require expansion to second order in the perturbation parameter, as a consequence of the assumed linearity of the molecular aligning effect of the flow in terms of its velocity gradient. In cases where the molecular interaction field is the negative gradient of a free energy function, such as the Landau-de Gennes fourth order potential, we give explicit criteria on the coefficients to ensure the existence of the stable kayaking orbit for sufficiently small contribution from the symmetric part of the flow gradient. The admissible size of this contribution is not estimated, so that care must be taken in interpreting experimental or numerical verification.

A.1 Bilinear maps
From (2.14) and equivariance it follows that the element R z (π) ∈ Σ z acts on each isotypic component V Z i by for v i ∈ V Z i , i = 0,1,2, and so from (A.5) we see that any Σ z -equivariant bilinear map B : Thus R z (π) fixes B(v i ,v j ) when i + j is even and multiplies it by −1 when i + j is odd. As a consequence we have the following result, extremely useful for simplifying calculations. With G(Q * ) = 0 this gives λ = 2aḡ * + ∆ g * + 2a∆ḡ * = 2aḡ * + ∆ĝ * (A.14) using Proposition A.2 and (A.10), where g * denotes g(Q * ) and ∆ g * := Dg(Q * )Q * etc.. In the main text we need to evaluate the component of this expression tangent to the SO(3)-orbit O of the uniaxial matrix Q * at points Z ∈ O. Here we calculate this for Z = Q * making significant use of Proposition A.3 and Corollary A.5, and will be able to transfer the result to a general Q = Z ∈ O by applying the SO(3) action.