Three-Dimensional Spontaneous Flow Transition in a Homeotropic Active Nematic

We study the three-dimensional spontaneous flow transition of an active nematic in an infinite slab geometry using a combination of numerics and analytics. We show that it is determined by the interplay of two eigenmodes -- called S- and D-mode -- that are unstable at the same activity threshold and spontaneously breaks both rotational symmetry and chiral symmetry. The onset of the unstable modes is described by a non-Hermitian integro-differential operator, which we determine their exponential growth rates from using perturbation theory. The S-mode is the fastest growing. After it reaches a finite amplitude, the growth of the D-mode is anisotropic, being promoted perpendicular to the S-mode and suppressed parallel to it, forming a steady state with a full three-dimensional director field and a well-defined chirality. Lastly, we derive a model of the leading-order time evolution of the system close to the activity threshold.


I. INTRODUCTION
Active matter is a class of materials that lie outside of thermodynamic equilibrium due to the conversion of energy consumed by the constituent particles to mechanical work [1][2][3].Active matter can be considered an active nematic whenever when the constituent particles display orientational order akin to a nematic liquid crystal [4].Such systems can be natural, such as cell colonies [5][6][7], epithelial tissues [8,9], bacterial suspensions [10][11][12], and microtubule and motor protein mixtures [13], or artificial, such as vibrated granular rods [14,15].A key property of active matter is the emergence of spontaneous, collective motion on scales much larger that that of the individual constituents.This has important real-world implications.In biology, for example, collective motion plays a role during organ formation and development [16] and wound healing [17].There is also potential to harness the self-generated flows of active nematic materials to create self-operating microfluidic devices that do not rely on external forcing, or to incorporate other aspects of passive liquid crystals, such as utilising colloidal inclusions [18][19][20].
Active nematic systems may be modelled by adapting the well-established dynamical equations of passive nematic liquid crystals [21,22] to include active terms [4,23].One key triumph of the theory of active nematics is the prediction that such systems will spontaneously transition to a flowing state on their own accord due to their fundamental hydrodynamic instability [23].In unbounded systems, this instability sets in at arbitrarily long perturbation wavelengths and the system eventually transitions to a chaotic state known as active turbulence [4,10,12,13,24].Confinement of active nematic systems can suppress the onset of active turbulence and instead the hydrodynamic instability acts to produce non-chaotic flows, first predicted theoretically by Voituriez et al. [25] and later confirmed in simulations performed by Marenduzzo et al. [26].The spontaneous flow transition has been observed in experiments on spindle-shaped cells in confined strips [27], demonstrating potential relevance to cell transport in development or cancer.
Of particular interest to us are the spontaneous flow transitions within rectangular channels.Different flow states can be found, depending the boundary conditions and parameters.Flows can be roughly separated into two categories: streaming flow states and swirling flow states [28].These two categories can be further subdivided.For example, the streaming flow category can be subdivided into Poiseuille-like flows [33,43], shear-like flows [27], oscillatory flows [33,43], grinder train flows and double helix-like flows [48].The latter two are only seen in three dimensions and possess non-zero helicity and are yet to be seen experimentally.
Here, we study the spontaneous flow transition for an active nematic in a three-dimensional cell with normal anchoring boundary conditions.This geometry is analogous to the Frederiks transition in a homeotropic cell.We find that the transition leads to a twisted director field and a spontaneous flow that has both Poiseuille-like and shear-like components.The twist is right-handed or lefthanded with equal probability and represents a spontaneous chiral symmetry breaking, in addition to the spontaneous rotational symmetry breaking of the direction of the Poiseuille-like flow.We identify the reason for this as the degeneracy of two eigenmodes of the linear stability operator for the system at the threshold of instability.We believe that this is an accidental degeneracy, rather than arising due to some underlying symmetry.We label these modes the S-mode and the D-mode.The degeneracy is accidental, rather than arising from an underlying symmetry, and clarifies some aspects of the existing literature for planar anchoring.We develop a hierarchical perturbative analysis of the growth of both modes above threshold that reproduces all aspects of the instability in excellent agreement with full numerical simulations.

II. SPONTANEOUS FLOW TRANSITION
We consider an extensile, uniaxial active nematic confined between two infinite, parallel plates with a fixed cell gap, d.We assume no slip boundary conditions and strong homeotropic anchoring on both plates.For this anchoring condition and with the normal of the plates being e z , the ground state (i.e. the state that the system is in below threshold) director field is n = e z which possesses evident rotational symmetry around the z axis.The setup is shown in Fig. 1.
We establish the basic character of the active instability and spontaneous flow transition by performing numerical simulations with random initial perturbations to the ground state.We find that there is a threshold in activity, below which the system remains in the ground state and above which the system spontaneously starts flowing.The flow field consists of a Poiseuille-like component and a shear-like component perpendicular to it.The Poiseuille-like flow component results in a net flux within the system, the direction of which is random and represents spontaneous rotational symmetry breaking.The director field is twisted with either a right or left handedness, occurring with equal probability.Hence, the system also undergoes spontaneous chiral symmetry breaking.We note that the shear-like component of the flow is reversed between the two possible twist configurations.The director and flow fields are shown in Fig. 1.
This twisted flow state arises from the coupled evolution of two degenerate eigenmodes that both become unstable at the activity threshold.We believe that this is an accidental degeneracy, rather than arising due to some underlying symmetry.The degeneracy may be lifted by applying a generic perturbation, such as the application of an electric field.We label these modes the S-mode and the D-mode.The different chiralities emerge from the fact that the D-mode can evolve in one of two possible directions perpendicular to the S-mode, with each direction being equally probable.The flow components associated with the S-mode and D-mode are the Poiseuille-like and shear-like flows respectively.

A. Active Nematic Hydrodynamics
Active nematic systems can be modelled by the active Beris-Edwards equations [4,21] which describe the coupled evolution of the fluid density, ρ, velocity, v, and the nematic order paramerter, Q.We solve the Beris-Edwards equations numerically using a hybrid lattice Boltzmann algorithm [49], with full details given in §VII.The activity is incorporated into (2) in the usual way by adding an additional contribution to the stress, Π a = −ζ LB Q, modelling a force dipole at the microscopic level with a strength given by the phenomenological activity parameter, ζ LB .Extensile activity corresponds to ζ LB > 0 and contractile activity to ζ LB < 0.
In the analytical analysis, we work in terms of the director field, reducing the Beris-Edwards nematodynamic equations to the Ericksen-Leslie form [26]. Assuming low Reynolds number, constant density, and a uniaxial form for the nematic order parameter, 2 n i n j −δ ij /3 , with constant S, one writes: In ( 5), σ el denotes the elastic stresses coming from the nematic director and the active stress is σ a = −ζnn, where ζ = 3S 2 ζ LB .We consider only the flow aligning regime, where the flow aligning parameter ν < −1 [22].Further relevant aspects of the correspondence between the Beris-Edwards and Ericksen-Leslie equations are given in §VII.

B. Linear Stability Analysis
We start by considering the Ericksen-Leslie formalism in quasi-one-dimensional geometry where the spatial dependence is only along the cell normal (z-direction) but the flow field, v, and active nematic director, n, can be in any 3D direction.The continuity equation then implies that v z = 0 and the Stokes equation ( 5) can be integrated directly to give Here, the notation the average of the argument over the cell gap; these terms arise from the no slip condition at the two cell boundaries, z = 0, d.For the director dynamics, we will find it convenient to write n in the form n = cos φ cos θ e z + sin θ e x + sin φ e y , parameterised by two angles θ and φ, in terms of which the director dynamics (6) becomes where we have defined the unit vectors m φ = − sin φ cos θ e z + sin θ e x + cos φ e y .
Substituting the flow solution (8) for D and Ω, (10) and (11) reduce to a pair of coupled, nonlinear integrodifferential equations for the two angles, which we give in full in Appendix A. Both equations have the same linearisation, which we write only for θ, The uniform state (θ = 0) is linearly unstable when the linear integro-differential operator, L, has a positive eigenvalue, λ.A perturbation along the associated eigenfunction then grows exponentially with rate λ until it saturates at a steady state solution of the full nonlinear equations.The eigenfunctions of L separate into two symmetry classes according to whether they are odd or even about the cell midplane.
For the odd eigenfunctions, the integral terms in ( 14) vanish and L reduces to a Schrödinger-type operator whose eigenfunctions are where A S is an amplitude.The lowest mode, n = 1, becomes unstable first, which happens at the threshold activity In the flow aligning regime (ν < −1), which we restrict our attention to, the instability is for extensile activity (ζ > 0).This unstable mode is associated with a flow that is even about the cell midplane and represents a fluid flux along a spontaneously chosen direction.We refer to this unstable mode, and the steady spontaneous flow state it evolves into, as the 'S-mode' due to the appearance of the director across the cell gap.
For the even eigenfunctions of L, the integral terms in ( 14) do not vanish and we have not found closed-form expressions for all of the eigenfunctions.However, one can verify directly that is an eigenfunction with eigenvalue zero at the threshold activity, ζ = ζ th .A D is an amplitude for the mode.The associated fluid flow is shear-like and odd about the cell midplane with no net flux.As before, the direction is chosen spontaneously.We refer to this unstable mode as the 'D-mode', again due to the appearance of the director across the cell gap.Numerically, we seeded an S-mode of the form ( 15) and a D-mode of the form (18) separately and letting them evolve into steady state for activities very close to the threshold.To do this, we reduced our simulation box to 3 bulk points in the direction perpendicular to the flow, effectively reducing the system to two dimensions so as to more easily isolate the individual eigenmodes.The results of the S-and D-modes are shown with their associated flow fields in Fig. 2, along with direct comparison to analytical predictions, where there is excellent agreement.
Overall, the spontaneous flow instability with homeotropic and no-slip boundary conditions is characterised by having two degenerate modes, one in each symmetry class, that become linearly unstable at the same threshold activity, each with a spontaneously chosen inplane direction.This degeneracy in the linear instability distinguishes the active spontaneous flow transition from the Frederiks transition in a passive system, where the instability is to the fundamental mode in the even sector at a threshold well below that of the first mode in the odd sector [22].

IV. GROWTH RATES ABOVE THRESHOLD
Above the threshold activity both unstable modes will grow exponentially at rates given by their respective eigenvalues of the linear stability operator L. We first determine these for the S-and D-modes separately and subsequently consider how they coevolve.The S-mode ( 15) is an eigenfunction of L for all values of the activity, with eigenvalue For the D-mode, the expression ( 18) is an exact eigenfunction only at the threshold activity, ζ = ζ th , where the eigenvalue is zero.We do not have its closed form more generally.However, for activities close to threshold we can determine the eigenvalue from perturbation theory, expanding to first order in ζ − ζ th .We provide the details in §VII and state here only the result Comparing against the eigenvalue of the S-mode gives a ratio λ S /λ , from which we see that the S-mode will grow fastest above threshold.As a result, unless it is suppressed, it will dominate the initial evolution of the system post instability.
To allow for comparison with numerical simulations, we convert (20) and ( 21) to simulation units.
To acquire the growth rates numerically, we individually seeded S-and D-modes with the forms given by ( 15) and ( 18) respectively, and then analysed their evolution at different values of ζ LB .To extract their linear growth rates, we plotted the maximum mode amplitude over time on a logarithmic scale and used the gradient to extract the mode's exponential growth rate.We plot the exponential growth rates vs activity in Fig 3, from which we can extract the growth rate coefficient and the activity threshold.
We see a very good agreement with both the growth rate magnitude and the threshold for both the S-and D-mode, with percentage differences not exceeding 5%.One notable discrepancy is the threshold prediction from the graph of the D-mode's growth rate, which is larger than that of the S-mode yet is predicted to be the same from the linear stability analysis.Around the threshold, we observe a plateau before it eventually starts to grow linearly, the fit of which gives us the threshold prediction.The simulations for a seeded D-mode were significantly more difficult and unstable compared to the S-mode, especially around threshold, making D-mode's threshold result unsurprising.This same anomalous behaviour for a seeded D-mode was also observed by Marenduzzo et al. [26].As expected, the growth rates eventually deviate from a linear relation at a large enough activity above threshold.
The difference in the numerical values of the two growth rates suggests a separation of timescales that allows us to treat the instability as effectively a two-stage process.Initially, the S-mode grows fastest and attains a finite amplitude and steady state, while the D-mode remains infinitesimal.Subsequently, the D-mode evolves on top of the established S-mode.As the S-mode spontaneously breaks rotational symmetry within the cell, the problem is no longer isotropic and we consider separately growth of the nascent D-mode parallel and perpendicular to the established S-mode.
We denote by θ * (z) the steady state solution of ( 10) with φ = 0, corresponding to a fully established pure S-mode.It is given by and reduces to the quadrature where the expression applies for 0 The amplitude of the mode is A * S = θ * (d/4), which may be obtained implicitly as a function of ζ by setting z = θ * = 0 in (25).We show this dependence in Fig. 4. The behaviour close to threshold has the square root form A * S ∼ (ζ − ζ th ) 1/2 , which may be found from an expansion of (25) (with z = θ * = 0) to linear order in A S .Explicitly, to leading order we find the amplitude is We now determine the growth rate of the D-mode, to linear order in ζ − ζ th , in the presence of a steady state S-mode.This amounts to retaining all terms up to O(A * S 2 ) from the steady-state S-mode in the linearised dynamics for the D-mode, which therefore modifies the growth rates as compared to (21).We consider separately the growth of the D-mode parallel and perpendicular to the (spontaneously chosen) direction of the established S-mode.For the perpendicular case we substitute θ = θ * (z), φ = δφ D (z, t) into (11) and linearise in δφ D .The calculation of the growth rate uses the same perturbation theory as before and is given in Appendix B 2. For the parallel case we substitute θ = θ * (z) + δθ D (z, t) into (10) and linearise in δθ D ; the analysis is again given in Appendix B 2. The two growth rates are which in simulation units read The main result is that λ ∥ < 0 and λ ⊥ > 0 so that the D-mode only remains linearly unstable along the direction perpendicular to that set by the S-mode.As a result, the director evolves into a truly three-dimensional configuration with the S-and D-modes growing along orthogonal in-plane directions.This interplay between the two modes leads to twisted director fields with where in the second form we have linearised in θ and φ and taken them to have the threshold forms ( 15) and ( 18), respectively.The twist maintains a single sign (handedness) throughout the cell, vanishing only on the two boundaries.Since the S-mode spontaneously breaks rotational symmetry in the plane, its amplitude A S is always positive.In contrast, the amplitude of the D-mode, A D , can be positive or negative (corresponding to the two directions orthogonal to the established S-mode); when it is positive the twist is right-handed and when negative it is left-handed.In a nematic material we expect both to occur with equal probability and any particular realisation represents a spontaneous chiral symmetry breaking.This general mechanism for confined active nematics may also be relevant to the emergence of twist in bulk three-dimensional systems [42,50] and possibly also to the prevalence of twist loops in the statistics of their defect loops [51,52].
The growth rate λ ⊥ (27) for the orthogonal D-mode can be verified numerically by initialising the director with an S-mode along e x and a small amplitude D-mode along e y .Tracking the exponential growth of the D-mode as a function of activity allows for a fit of the growth rate and threshold activity as before.This is shown in Fig. 5.The agreement with the theoretical prediction is again excellent.We note, particularly, that we obtain better agreement for the threshold activity ζ th than we found from simulations with only the D-mode.

V. MODE EVOLUTION AND STEADY STATE
We now summarise and describe the full evolution of the instability to the steady spontaneous flow state.This can be studied systematically in numerics by seeding a small amplitude director perturbation consisting of an S-mode along e x and a D-mode along e y and tracking their amplitudes -the maximum values of θ and φover time.This is shown in Fig. 6.The evolution can be divided into three distinct regimes: in the first (I), there is exponential growth of both modes, but with the S-mode growing significantly faster.This corresponds to the independent and isotropic mode dynamics described in Fig. 3.In the second regime (II), the S-mode amplitude attains a plateau and there is an increase in the exponential growth rate of the D-mode.The S-mode amplitude at its plateau corresponds to the value A * S described in §IV and the enhanced growth rate of the D-mode is the cross-over to the rate λ ⊥ as described in Fig. 5. Finally, in the third regime (III) the D-mode amplitude attains its steady state value and promotes a small further increase of the S-mode amplitude to its steady state.
This joint evolution can be cast as a coupled dynamical system for the amplitudes A S , A D of the S-and D-modes where the growth rate functions g S , g D have the fixed point structure shown schematically in Fig. 7.We define A S to be strictly positive, meaning that A D can take either sign.There are three important fixed points: the origin and the two points corresponding to the right-and left-handed states.Below threshold, the origin is a stable fixed point, but above it becomes unstable to all S-and D-mode perturbations.Depending on the sign of A D , perturbations around the origin will either flow to the left-handed or the right-handed stable fixed points, corresponding to the handedness of the resulting flow state's chirality.These flows are shown by the blue and red dashed lines in Fig. 7.The trajectory follows closely to the A S axis until A S is large and then rapidly moves away from the axis to the fixed points.This corresponds physically to the S-mode growing to a large amplitude before there is any significant D-mode growth.Finally, we note that in the absence of a D-mode, the S-mode grows to the semi-stable fixed point labelled as A * S , which has a slightly smaller amplitude than the left-and right-handed fixed points.This fixed point is described in 25.For activities close to threshold, the fixed points are close to the origin and we can expand the growth functions as where the non-linear terms are those allowed by symmetry.We give the calculation of the Λ coefficients in Appendix C.This system connects to our previous results and reproduces the numerically observed dynamics of Fig. 6.The amplitude A * S of the S-mode plateau in regime II is given by (λ S /Λ 1 ) 1/2 and matches the value in (26).Similarly, the enhanced growth rate of the D-mode in regime II is given by λ D + (Λ 3 /Λ 1 )λ S and matches the rate λ ⊥ in (27).At this leading order, we obtain the steady state amplitudes We note that This coincides with the numerical observation that there is an increase in the S-mode amplitude when the D-mode comes into steady state.In the numerical observations, this increase in small, implying that the D-mode coupling to the S-mode evolution is small, i.e.Λ 2 /Λ 1 ≪ 1.Indeed, in simulation units, Λ 2 /Λ 1 = 0.08.With this weak coupling, the evolution of the S-mode can be approximated by which gives We can substitute this into the leading order part of g D and solve to give the approximate the evolution of A D as where 2 F 1 is Gauss's hypergeometric function.We compare (40) and (41) with numerical data of mode amplitude evolution in Fig. 6 for ζ = 0.0675.The analytical model captures the qualitative, triphasic nature of the system well.For the S-mode, the agreement is very good up to phase III, where the second plateau is not captured due to the decoupling approximation that was made.For the D-mode, the analytical prediction of the growth rate in phase I, (23), is larger than the numerics, although this is consistent with what we have already observed when we seeded an independent D-mode.Recall that the numerical isotropic growth rate for a Dmode has a shifted threshold compared to (23) which has a significant effect on the growth rate close to threshold, thus making the numerical growth rate noticeably smaller than what is analytically predicted.Nevertheless, the analytical prediction tracks the numerical data very well thereafter, albeit translated upwards due to the first phase growth rate being too large.
As we move further and further above threshold, the analytical model becomes worse and worse.If, however, θ and φ are in odd end even symmetry classes respectively, then the functions ∂ t θ and ∂ t φ are also odd and even respectively which can be checked by inspection of each term in (10) and (11).This implies that if θ and φ start off as odd and even respectively, then the functions will remain in the same symmetry class for the entirety of their non-linear, coupled time evolution.Furthermore, we can analyse the symmetry of v and upon inspection of the formula, if θ and φ are in their aforementioned symmetry classes then v x will be an even function and v y will be and odd function.Hence, we expect that the generation of a chiral director field is a general property of the system, rather than just a feature close to threshold.

VI. DISCUSSION AND CONCLUSIONS
We have studied a spontaneous flow transition in an active nematic fluid with an infinite slab geometry and normal surface anchoring.We find the existence of two independent flow instabilities, the S-mode and the Dmode, that occur at the same threshold but have different growth rates above threshold, which we calculate using perturbation theory of a non-Hermitian integrodifferential operator.Above threshold, the S-mode with its larger growth rate grows in a random direction to steady state, breaking the initial rotational symmetry of the system.Thereafter, any perturbations within the system are subject to an anisotropic environment.In particular, D-mode perturbations that are parallel and perpendicular to the direction of anisotropy decay and grow respectively.The growth of the D-mode perpendicular to the S-mode yields a steady-state with a full, threedimensional chiral director field, with spontaneously broken chiral symmetry.We analytically describe the mode growth with a leading-order model that captures the key characteristics of the system and its evolution into the spontaneously flowing state.
The first natural extension of this work is to explore the possible inhomogeneity of the flow within the cell plane, enabling the study of the chiral flowing state's stability to the Goldstone mode coming from the breaking of rotational symmetry and the umbilic defect lines associated to this.Furthermore, in large systems, the spontaneous chiral symmetry breaking could yield domains of different chirality, leading to an effective non-conserved binary mixture.We hope the novel three-dimensional flow instability we have uncovered can provide motivation for experimental research of active nematic systems with normal anchoring.
Active nematics are fundamentally analogous to passive (i.e.non-active) nematic liquid crystals, with the orientational ordering of the anisotropic material building blocks crucially determining the material dynamics, including at the surface.Today, surface anchoring in passive nematics can be realised experimentally in different configurations, ranging from uniform planar and degenerate planar to homeotropic and even tilted and tilted degenerate.Such advanced control over surface alignment was shown together with control over confinement geometries to diverse material structures and phenomena, such as realisation of colloidal and field knots [53][54][55], self-assembly [56], hexadecapolar and 32-pole field configurations [57,58], memory [59,60], static and dynamic solitons [61][62][63], tunable positioning of topological defects [64,65].Clearly, advancing the ability to control different surface anchoring regimes in combination with confinement [66] could open diverse research directions in confined active nematics, especially at the experimental level.Even rather simple surfaces like spheres imposing homeotropic anchoring on active nematics would lead to the emergence of topologically imposed and conditioned bulk Saturn ring defect states that are commonly ob-served in passive nematics but not seen in active nematics.For example, could combining three-dimensional selfassembly of (active) nano-objects (e.g., bacteria) combined with confinement lead to different effective anchoring beyond the typical planar?[67] Another possibly emerging direction from this work is the observation of the spontanenous chiral symmetry in an active (nematic) system.At least in non-active nematic materials, chiral symmetry breaking naturally leads to the emergence of different phenomena, such as chiral domains and associated topological defects such as solitonic-like nematicons [68].In active systems, chiral activity can lead to topological edge modes [69] and odd elastic responses [70].Combining spontaneous chiral symmetry breaking with activity in more advanced geometries and setups -beyond the simple homeotropic cells studied in this work -could lead to an exciting advancement of the control and design of novel active functional matter.

A. Numerical Methods
We solve the Beris-Edwards equations ( 1)-(3) using a hybrid lattice Boltzmann algorithm.In these equations, Π is the stress tensor defined by where p is the pressure, µ is an isotropic shear viscosity, χ is the flow alignment parameter, Γ is a rotational viscosity and ζ LB is the activity parameter.
are the components of the molecular field, in which is the free energy, where The bulk part of the free energy is described via phenomenological constants for phase transition A, B, C, one constant approximation for the elastic part L, and the surface part is described via homeotropic anchoring with strength W h and orientation Q 0 ij corresponding to the preferred perpendicular director on the surface plate n = ±e z .Finally, the tensor S with components describes the coupling between the nematic order parameter and the flow field.We use D and Ω to represent the symmetric and antisymmetric parts of the velocity gradient tensor respectively.
In the lattice Boltzmann simulations, the fundamental scaling parameters are nematic correlation length Bs eq + Cs 2 eq ) ( 47) and nematic time scale All of the parameter values are expressed in units of the elastic constant L. We used the following Landau-de Gennes parameters A = −0.687L/ξ Under the aforementioned choice of parameters s eq = 0.651.This corresponds to the use of Beris-Edwards parameters χ = 1, µ = 1.38/Γ and ρ = 0.031 1 LΓ 2 .We performed simulations with a cell size of 201 × 201 × 45 bulk points confined between the two plates using homeotropic anchoring with strength W h = 2 3 L/ξ n and periodic boundary conditions on the side.We used disretisation of spatial coordinates ∆x = 1.5 ξ n and a time step ∆t = 0.025 τ n .
The nematic director, n, is obtained as the eigenvector associated with the largest eigenvalue of Q, which represents the magnitude of the order parameter, S.

Eigenvalues above Threshold
We solve for the leading order correction to the eigenvalue, λ, of the stability operator equation Lθ = λθ when ζ > ζ th .To do so, we perform a perturbation expansion, L = L 0 + L 1 + . . .(and likewise for λ and θ), in powers of ζ − ζ th , with L 0 corresponding to L when ζ = ζ th (and likewise for λ 0 and θ 0 ).In the main text, L is defined by (14).However, it is important to note that the analysis provided in this section is valid for all operators, L, that linearise to when ζ = ζ th .We proceed in the usual way of expanding Lθ = λθ and equating order by order, giving . . .L 0 is non-Hermitian over the interval z ∈ [0, d] due to the ⟨∂ 2 z • ⟩ term, meaning that we cannot apply standard methods for Hermitian operators to solve for λ 1 ; we must instead take a slightly different approach that is unique to operator (51).Firstly, we take the average of all in (53) and rearrange to get Next, we multiply (53) by θ 0 and average, giving To simplify this, we make use the result where the second equality comes from the fact that ⟨θ 1 L 0 θ 0 ⟩ = 0 and ∂ 2 z θ 0 = 0. We can combine (54), ( 55) and (56) to eliminate ∂ 2 z θ 1 and solve for λ 1 , giving We note that in the case where the mode is antisymmetric about the cell midplane, all of the integral terms in (51) vanish, making L 0 Hermitian over the interval z ∈ [0, d].
As expected, (57) reduces to λ 1 = ⟨θ 0 L 1 θ 0 ⟩ / θ 2 0 , the standard result of perturbation theory on Hermitian operators.We use (57) to calculate the isotropic and anisotropic growth rates, with calculation details given in Appendix B.

Amplitude Evolution above Threshold
We obtain the leading order time evolution equations of the S-and D-mode amplitudes close to threshold.We take θ to be an S-mode and φ to be a D-mode, and then make the usual expansion θ = θ 0 + θ 1 + . . .and φ = φ 0 + φ 1 + . . . in powers of ζ − ζ th .We substitute these expansions into (10) and ( 11) and equate order by order.The leading order balance gives us the now familiar eigenfunctions which we write as θ 0 = A S (t) sin 2πz d := A S (t)ψ S (z) and φ 0 = AD(t) The subsequent analysis is given in terms of θ, but is the same for φ.At the next order in ζ − ζ th , we obtain We can now solve for dA dt in the same way as we solved for λ 1 in the previous sub-section, giving us This is the equation used to obtain the leading order parts of ( 34) and (35) in the main text.We go through the details of the calculations in Appendix C.

ACKNOWLEDGMENTS
For the purpose of open access, the authors have applied a Creative Commons Attribution (CC-BY) licence to any Author Accepted Manuscript version arising from this submission.We acknowledge support from Nordita during the Program on Current and Future Themes in Soft and Biological Active Matter.V.J.P. acknowledges funding from the EPSRC, grant EP/T51794X/1.M.R. and E.C. acknowledge funding from the Slovenian research agency ARRS grants P1-0099, N1-0195 and J1-2462, and EU ERC AdG LOGOS.
Next, we substitute these results into equation (57) to give This is the result given in equation (21) in the main text.

Anisotropic D-mode Growth Rates
In this section, we calculate the exponential growth rates of D-mode perturbations on top of an anisotropic steady state due to an S-mode established in the spontaneously chosen direction, e x .The S-mode is denoted by θ * (z) and is the solution to equation (24).We consider the two cases: a D-mode perturbation to φ, θ = θ * (z) and φ = δφ D (z, t), and a D-mode perturbation to θ, θ = θ * (z)+δθ D (z, t).In the small angle regime close to threshold, these perturbations have leading order contributions along e y and e x respectively, thus we interpret them as perpendicular and parallel perturbations to the direction of anisotropy.

a. Perpendicular Growth Rate
We substitute θ = θ * (z) and φ = δφ D (z, t) into equation (A2) and linearise about δφ D , giving The right-hand side is the stability operator, L, of δφ D .In the limit ζ → ζ th , θ * → 0 and the right-hand side of equation (B4) reduces to the form given in equation ( 51), meaning that we can apply the perturbation theory outlined in §VII.Close to threshold, θ * = A * S sin 2πz d + higher order terms, where . Therefore, upon expanding the stability operator in θ * , we obtain the expansion L = L 0 + L 1 + . . ., with terms O(θ * ) 2 contributing to L 1 .(B5) Only the leading order term of θ * contributes to this order.Thus, to calculate the growth rate, we evaluate equation (57) with δφ D,0 ∝ 1 − cos 2πz d and θ * = A * S sin 2πz d .The result is given by equation (27) in the main text.
b. Parallel Growth Rate The calculation is the essentially the same as that given in Appendix B 2 a, in which more detail is given.We substitute θ = θ * (z) + δθ D (z, t) and φ = 0 into equation (A1) and linearise about δθ D , giving (B6) From this, we expand in θ * to obtain Following the same method as in Appendix C 1, we expand (A2) up to cubic order, giving from which we extract In this case, equation ( 59) is evaluated as

FIG. 1 .
FIG. 1. Director and flow fields of the spontaneous flow transition.(a) Left: steady-state director and velocity fields after the spontaneous flow transition with left-handed chirality.Right: decomposition into the x and y components of the left-handed steady-state director and velocity fields.(b) Left: steady-state director and velocity fields after the spontaneous flow transition with right-handed chirality.Right: decomposition into the x and y components of the right-handed steady-state director and velocity fields.

FIG. 2 .
FIG. 2. Comparison of the simulation results of the director and flow profiles (solid blue and red curves) with the corresponding analytical predictions (black dashed curves) for activity in the vicinity of the threshold activity, ζLB = 0.063.The blue curves show the director (left) and velocity (right) profiles of the S-mode, so called because of its 'S'-like appearance across the cell gap.The red curves show the director (left) and velocity (right) profiles of the D-mode, again, named after the 'D'-like appearance of its director profile.The observed amplitude of the S-mode profile is AS = 0.0043 and the velocity amplitude is vS = 1.71 • 10 −4 ΓL/ξn.Observed ratios between S-and D-mode are AD/AS = 0.95 and vD/vS = 0.25.

FIG. 3 .
FIG. 3. Exponential growth rates vs activity for an individually seeded S-mode (blue), and D-mode (red).The black dashed lines are linear fits of the numerical data close to threshold, where the growth rates exhibit a linear relation with activity.

FIG. 4 .
FIG. 4. Plots of the S-mode amplitude against ζ/ζ th .The blue line shows the quadrature solution (25) and the black, dashed line the leading order part of the expansion of A *S , given by(26).The black data points show numerical data of an individual S-mode in steady state, using the reduced simulation box of 3 bulk points perpendicular to the direction of flow.The numerical data is scaled with a threshold activity of ζ th = 0.062 to fit the analytical prediction, again in good agreement with the theoretical prediction.

FIG. 5 .
FIG. 5. Plot of the D-mode's growth rate perpendicular to an established S-mode.The black dashed line is a linear fit of the numerical data in the selected range of data close to the threshold.

FIG. 6 .
FIG. 6.(a) Comparison of the numerical evolution of the S-mode amplitude and D-mode amplitude (blue and red lines, respectively), with the analytical predictions of the leading order amplitude evolution from (40) and (41) (black, dashed lines) at ζLB = 0.0675.The plot has a logarithmic scale on the vertical axis.(b) Left: S-mode profile along the cell gap at three different times of the time evolution (I, II and III, as seen in (a)).Right: D-mode profile along the cell gap at the same three time points.We used time frames 0.3 × 10 5 τ , 1.2 × 10 5 τ and 2.7 × 10 5 τ for I, II and III respectively.

FIG. 7 .
FIG.7.Schematic of the phase portrait of the spontaneous flow transition above threshold.The vertical blue line is the S-mode trajectory where the S-mode grows independently of the D-mode, which is the approximation made in(39).The red and blue dotted lines show the trajectory of the coupled evolution of the S-and D-modes to the left-and right-handed chiral states.