Pattern selection and the route to turbulence in incompressible polar active fluids

Active fluids, such as suspensions of microswimmers, are well known to self-organize into complex spatio-temporal flow patterns. An intriguing example is mesoscale turbulence, a state of dynamic vortex structures exhibiting a characteristic length scale. Here, we employ a minimal model for the effective microswimmer velocity field to explore how the turbulent state develops from regular, stationary vortex patterns when the strength of activity resp. related parameters such as nonlinear advection or polar alignment strength—is increased. First, we demonstrate analytically that the system, without any spatial constraints, develops a stationary square vortex lattice in the absence of nonlinear advection. Subsequently, we perform an extended stability analysis of this nonuniform ‘ground state’ and uncover a linear instability, which follows from the mutual excitement and simultaneous growth of multiple perturbative modes. This extended analysis is based on linearization around an approximation of the analytical vortex lattice solution and allows us to calculate a critical advection or alignment strength, above which the square vortex lattice becomes unstable. Above these critical values, the vortex lattice develops into mesoscale turbulence in numerical simulations. Utilizing the numerical approach, we uncover an extended region of hysteresis where both patterns are possible depending on the initial condition. Here, we find that turbulence persists below the instability of the vortex lattice. We further determine the stability of square vortex patterns as a function of their wavenumber and represent the results analogous to the well-known Busse balloons known from classical pattern-forming systems such as Rayleigh–Bénard convection experiments and corresponding models such as the Swift–Hohenberg equation. Here, the region of stable periodic patterns shrinks and eventually disappears with increasing activity parameters. Our results show that the strength of activity plays a similar role for active turbulence as the Reynolds number does in driven flow exhibiting inertial turbulence.


I. INTRODUCTION
One of the recurring questions in fluid dynamics research is how turbulence arises from a laminar base flow upon increase of the Reynolds number [1,2].In driven fluid systems described by the Navier-Stokes equation, the geometry of the problem and the specific setup have a crucial influence on the relevant mechanisms.For example, in Taylor-Couette flow dominated by inner-cylinder rotation, the development of turbulence follows from a series of subsequent linear instabilities [3][4][5].This is in contrast to dominant outer-cylinder rotation, where linear instabilities are absent and, instead, finite perturbations are necessary to trigger turbulent flow [5].The situation is similar to channel [6], pipe [2,4], and plane Couette flow [7][8][9] where mechanisms beyond linear stability are at work [10].Recent largescale experimental and numerical studies of such systems [6,8,[11][12][13] have revealed that the mechanism is better understood in terms of a "statistical phase transition".For pipe flow, the dynamics of turbulent puffs, e.g., appears analogous with the much-studied nonequilibrium phase transitions in the directed percolation universality class [14].
Turbulent states are also a central topic in active matter, a research field that focuses on living or artifical systems consisting of moving constituents [15][16][17][18][19][20].Thus, active matter systems are intrinsically out of equilibrium.As a consequence, they display diverse pattern formation, including turbulent flows commonly labelled as "active turbulence" [21] that were found experimentally in suspensions of microswimmers [22,23] and in active nematics [24].Contrary to the above-mentioned turbulence in "passive" fluids, which arises at high Reynolds numbers, active turbulence occurs in the overdamped regime.Still, similarly to some of the passive systems, a recent study on the transition to active turbulence in a model for active nematic liquid crystals (such as microtubule protein mixtures) has revealed analogies to directed percolation as well [25].Recent studies of detailed "microscopic" models for the collective dynamics of interacting self-propelled microswimmers [26,27] as well as a qualitative model with competing alignment interac-tions [28,29] have also shown turbulent states with qualitatively similar energy spectra and correlation functions as continuum models for active polar fluids [23,30,31].
Besides active turbulence, stable vortex lattices with different symmetries have been reported in earlier experiments with swimming sperm cells [32] and an active microtubule-myosin motor mixture [33] and attributed to chiral motion of the active agents.Subsequently, many simulation studies of continuum models for active nematics [34][35][36] as well as of agent-based simulations [28,29] and continuum models for active polar fluids [37] exhibited vortex patterns of different symmetries.While most of continuum model studies assume constant density, recent work also has shown that vortex lattices precede turbulent clustering also in a model for compressible active fluids that allows for density fluctuations [38].A stable vortex lattice could also be realized experimentally by introducing a periodic array of small obstacles into a bacterial suspensions that exhibits active turbulence [39].This phenomenon could be reproduced with simulations in a model of an active polar fluid [40,41].More recently, several numerical studies have suggested control schemes to tame and control turbulence in active nematics [42,43].
Most experimental studies of transitions between regular patterns and turbulence in active systems have been carried out exploiting confinement effects, i.e., gradually changing dimensions and aspect ratio of the experimental system.In two-dimensional systems, a single vortex was shown to be stable in experiments of small systems [44] often connected with edge currents [45] and transitions to more complex multiple-vortex states when increasing the system size [46][47][48][49].In three-dimensional channels several experimental and theoretical studies with active nematics have shown that widening the smallest extension of the channel can lead to a transition between coherent patterns and apparently turbulent flows [50][51][52].For the polar active fluid, simulations of the transition between turbulent states and regular vortex patterns, which are externally stabilized by small obstacles in the flow, have been shown to exhibit features of a non-equilibrium phase transition in the Ising universality class [41].The latter study has utilized a model for suspensions of microswimmers (such as Bacillus Subtilis) that supports regular structures even in the absence of external stabilization, as well as turbulent states [40].To summarize, while some knowledge has been gained on the transition to turbulence in active fluids, we are far from a general understanding, including effects of the type of order parameter and the role of confinement and geometrical constraints.In particular, the situation for (polar) active fluids in the absence of spatial constraints is essentially unclear.This problem is addressed in the present paper.In particular, we tackle the question of how the turbulent state arises from an underlying regular pattern.
The utilized model is formulated in terms of an effective microswimmer velocity field v(x, t) and was first introduced on phenomenological grounds in Ref. [23].It captures the main features of mesoscale turbulence, a state of highly dynamic vortex patterns that has been observed in suspensions of bacteria [23,53] and artifical swimmers [54].The emergent flow patterns are characterized by the presence of a characteristic length scale much larger than the size of single swimmers but smaller than the system size [23,53].Recently, the model was substantiated by a derivation from a generic microscopic theory [55,56].Although there are a large number of follow-up studies using the model [31,37,40,[57][58][59][60][61][62], it is still unclear how mesoscale turbulence emerges from the underlying regular structures that the model supports.To answer this question, we here focus on the stability of possible structures, which include stripe patterns, as well as square, triangular and hexagonal lattices of vortices.
As in the Navier-Stokes equation, the emergence of turbulent states is facilitated by nonlinear advection.In analogy to the Reynolds number in classical turbulence, we identify the strength of self-advection (advection of the velocity field via the field itself) as the central parameter that determines the stability of regular flow patterns and, thus, the emergence of the active turbulent state.Analytical investigations of this connection, however, are severely complicated by the fact that the regular base state is already nonuniform.By means of an extended stability analysis, we are able to show that the mesoscale-turbulent state arises as a consequence of a linear instability of the underlying base pattern.The key finding here is that the route to turbulence includes multiple modes that only grow simultaneously via mutual excitement.Analogously to a critical Reynolds number, we then determine the critical advection strength, which signifies the onset of the instability.Our results are in excellent agreement with numerical solution of the full nonlinear equation (see Fig. 1(b) for examples).In parallel, we unravel the role of active polar alignment.As it turns out, the alignment strength can serve as an alternative control parameter to drive the onset of turbulence.To give a first impression, we present in Fig. 1(a) a schematic overview of the behavior at fixed advection strength.Upon increasing the alignment strength from small values, there is first a transition between the disordered (isotropic) state and the regular vortex lattice.The latter becomes unstable at a critical alignment strength, beyond which a turbulent state emerges.Going back, we find a hysteresis loop that, as our results show, represents a characteristic feature related to the onset and disappearance of mesoscale turbulence.Combining the results from the extendend linear stability analysis with a numerical investigation of the hysteretic behavior, we construct a state diagram of the system.We also discuss the stability of vortex lattices with arbitrary wavenumbers on the basis of stability balloons.We here employ a as the control parameter and distinguish between different spatiotemporal states using the root-mean-square velocity observed in the system.Solid lines represent stable branches, dashed lines indicate unstable branches.Increasing a, the isotropic state becomes unstable against the emergence of a regular square vortex lattice at a = 0.This regular pattern, in turn becomes unstable at ac and the system develops turbulence.Following the turbulence branch while decreasing a, we observe a hysteresis loop as indicated by the arrows.(b)/(c) Representative snapshots of the vorticity field ω(x, t) in the square vortex lattice and in the mesoscale-turbulent state at a = 0.5 and λ = 3.0 obtained via numerical solution of Eq. ( 1).The superimposed arrows show the velocity field v(x, t) and the size of the snapshots is 12π × 6π.

II. MODEL
We utilize an established model for dense suspensions of microswimmers [23,30,31,[55][56][57][58][59][60] where the dynamics is described via the effective microswimmer velocity field v(x, t) [56].The model can be derived from microscopic Langevin dynamics including coupling to the solvent flow [56] and has been shown to capture experiments on bacteria in unconfined bulk flow [23] as well as in the presence of geometric confinement, e.g., lattices of obstacles [40].The dynamics of v(x, t) can be conveniently written as where the functional F is given by Additionally, we assume that v(x, t) is divergence-free, ∇•v = 0. Thus, the pressure-like quantity q acts as a Lagrange multiplier ensuring the incompressibility.Note that F plays a role similar to a free energy density, however, as the activity of the swimmers makes this a fundamentally non-equilibrium system, it is in fact not a free energy, but rather a more general non-equilibrium potential.This model is sometimes referred to as the Toner-Tu-Swift-Hohenberg model [21], because characteristic terms from both the Toner-Tu and the Swift-Hohenberg equations appear in the above equations.As Eq. ( 1) conveniently shows, the dynamics of v(x, t) is governed by a competition between gradient dynamics determined by the functional F and nonlinear advection λv • ∇v.We denote the coefficient λ as advection strength, which can be related to mesoscopic parameters such as the self-propulsion speed [41,56].The coefficient a can be positive or negative and encodes the strength of polar alignment of the microswimmers.Generally, it increases with the activity (self-propulsion speed and active stresses) [41,56], which leads to pattern formation for a > 0, as we will see in the next section.The quartic term in the functional F, see Eq. ( 2), yields a cubic term in the evolution equation (1).This term has a negative sign and, thus, leads to the saturation of the emerging patterns.Note that in most studies [23,31,37,40,[57][58][59][60][61][62], there is an additional positive coefficient in front of the cubic term.However, this coefficient can be scaled out by an appropriate rescaling of v and λ.
Before starting with an in-depth analysis of the model, let us briefly discuss the mesoscale-turbulent state obtained via numerical solution of Eq. (1), see Appendix A for information on the numerical methods.Figure 1(c) shows a typical snapshot of the vorticity field ω = (∇×v) z in the turbulent state at a = 0.5 and λ = 3.16.The disordered vortex patterns display a characteristic length scale or vortex size, as observed in experiments on bacterial suspensions [23,53].

III. STATIONARY PATTERN SELECTION
Before discussing the onset of the turbulent state from an analytical point of view, we start by establishing the kinds of stationary states that can be expected in our model.First, let us note that Eq. ( 1) always has an isotropic, stationary homogeneous solution, v = 0. Depending on the value of a, one further finds a polar stationary homogeneous solution, i.e., |v| = √ a − 1.The transition to the polar state at a = 1 is analogous to the transition to a state of collective motion in the Toner-Tu equations [63,64].The emergence of such a polar state breaks the rotational symmetry of the system.In this work, we restrict ourselves to the case a < 1, hence we can assume rotational invariance.
Performing a linear stability analysis of the isotropic solution of Eq. ( 1) is straightforward, see Refs.[30,60] for more details.We start by adding a small perturbation to the isotropic solution, i.e, where δv and δ q are the perturbation amplitudes, k denotes the wavevector and σ is the growth rate to be determined.We here take into account a perturbation of q because the incompressibility condition is part of the dynamical system.After linearization around the isotropic state, we obtain the growth rate as a function of the wavevector, i.e., As a first observation, note that the growth rate is real-valued.This indicates that the instability is not oscillatory and instead stationary patterns start to grow.Due to the chosen scaling of Eq. ( 1), the maximum of σ(k) is located at a critical wavenumber k c = 1.The sign of the coefficient a determines whether the isotropic state is linearly stable.For a < 0, i.e., low activity, the growth rate remains negative for all k.Increasing a above the critical value a = 0, the growth rate becomes positive at k c , signifying the onset of a finite-wavelength instability.Subsequently, for high activity, i.e., a > 0, a band of unstable modes develops, limited by the wavenumbers Note that, for the here considered case a < 1, the growth rate only depends on |k| and not on the direction of the wavevector due to rotational invariance.Building on these results, we continue the analysis by determining what kind of patterns might actually emerge.As we are dealing with a stationary instability, we will restrict ourselves to stationary patterns.Further, we will set λ = 0 for now, thus eliminating the nonlinear advection term.As a result, the instabilities that lead to the turbulent state are disregarded and the emerging states are completely determined as minima of F, see Eq. (1).In the following, we will investigate different patterns with the help of derived amplitude equations to determine if they are indeed minima.
We start by constructing possible emergent states in the system on the basis of the result from the linear stability analysis.As we have seen, for a > 0 the isotropic solution v = 0 is unstable to perturbations characterized by a band of wavenumbers, with k c denoting the critical as well as the fastest-growing mode.Therefore, close to the instability, it is reasonable to assume that the emerging patterns will be dominated by a characteristic length scale given by the wavelength Λ c = 2π/k c .As the system is translationally invariant, we can arbitrarilly set the phase and, thus, a growing mode can be represented in two dimensions as where A x and A y are the amplitudes in x-and y-direction, c.c. denotes the complex conjugate and k • k = k 2 c .The incompressibility condition, ∇ • v = 0, imposes a constraint on the relation of the amplitudes A x and A y .To meet the condition we set Thus, the amplitude can be expressed by a single quantity A. In the following we will analyze the stability of patterns that can be characterized by the superposition of N modes of the form given in Eqs. ( 6) and (7).Each mode i = 1, . . .N is characterized by its amplitude A i and a polar angle φ i relative to the x-axis.Thus, the wavevector corresponding to that mode can be written as where k c = 1.The patterns are then represented via In the following, we will restrict the considerations to three modes, i.e., N = 3, representing different directions of k.
As a consequence, we only search for the minimum of the functional F in this predefined space.However, as we will see, three modes are sufficient to represent a wide variety of patterns.Due to rotational invariance, we can set one of the angles to an arbitrary value without loss of generality.Here, we use φ 1 = 0.The remaining five parameters A 1 , A 2 , A 3 , φ 2 and φ 3 now completely describe the pattern under investigation.For example, we can depict a stripe-like pattern via This is shown in Fig. 2(a), where we plot the vorticity field ω = (∇ × v) z as color map and superimpose the velocity field (v x , v y ) using arrows.Fig. 2(b) visualizes a "square lattice", obtained by setting φ 2 = π/2 and A 3 = 0, i.e., The so-described square pattern consists of clockwise-and counterclockwise-rotating vortices that are arranged in an antiferromagnetic manner.Thinking of cogwheels, such an antiferromagnetic arrangement makes intuitive sense due to the mutual friction of rotating vortices.Further, a "triangular lattice" is obtained, e.g., by setting φ 2 = 2π/3 and φ 2 = 4π/3.Again, as is visualized in Fig. 2(c), neighboring vortices always display opposing senses of rotation.Note that there are different possibilities to represent such square and triangular patterns.The common property of these representations is that the angle between the two modes in the square pattern is a multiple of π/2, whereas the angles between the modes in the triangular pattern is a multiple of π/3.Further, note that by phase-shifting one of the modes in the triangular lattice by π/2, we obtain a pattern that exhibits the symmetries of a Kagome lattice, a structure that has been explored in the context of spatially constrained microswimmer suspensions in an earlier work [40].
We will now exploit the fact that the dynamics without nonlinear advection (λ = 0) is completely determined by the relaxation to the minimum of the functional F, compare Eq. ( 1).This feature can be used to derive a system of evolution equations for the amplitudes A i and angles φ i .To this end, we insert the representation of three superimposed modes given by Eq. ( 9) into Eq.( 2).We will refrain from writing down the full form resulting for f as it is quite involved.However, let us state a few points to motivate the final result.First, we already ensured incompressibility, i.e., ∇ • v = 0, see Eq. (7).Thus, the first term in f [given by Eq. ( 2)] vanishes.Second, the term |(1 + ∇ 2 )v| 2 vanishes for modes characterized by k c , as is also demonstrated by the fact that σ(a = 0, |k| = k c ) = 0 [see Eq. ( 4)].This also holds for a linear superposition of multiple modes, provided all have wavenumber k c .The remaining terms appearing in f , i.e., −a|v| 2 /2 + |v| 4 /4, are of even order.Thus, inserting Eq. ( 9) into Eq.( 2) yields only even modes of the form e 2ilk•x , where l is an integer number, i.e., In the following, we neglect higher order modes with |l| ≥ 1.This simplification is justified by the wavenumber dependence of the growth rate [see Eq. ( 4)], which exhibits strong damping for larger wavenumbers due to the quartic term, i.e., −|k| 4 .As a result of this simplification, f loses its dependence on the spatial variable x.Instead of the functional F[v, ∇ 2 v], we now have in lowest order the function f 0 (A 1 , A 2 , A 3 , φ 2 , φ 3 ) determining the dynamics, see Appendix B for the explicit form.This approach has thus reduced the phase space to five dimensions.Recall that the five variables A 1 , A 2 , A 3 , φ 2 and φ 2 were introduced in the ansatz for the emerging patterns as a superposition of modes characterized by k c [see Eq. ( 9)].The precise form of f 0 thus depends on the type of patterns considered.
Building on this reduction of phase space, we are now able to derive evolution equations for the amplitudes and angles.To this end, we first express the dynamic equation for v(x, t) [Eq.( 1)] in terms of f 0 instead of F. In this context, we again set λ = 0.As the explicit function this equation can be utilized to obtain evolution equations for A i and φ i , see Appendix B. The final equations are: Possible stationary patterns emerging in our original system described by the evolution equation for the effective microswimmer velocity field are now given as stationary solutions of Eqs.(13), provided λ = 0. Recall that Eq. ( 1) is governed by gradient dynamics, i.e., a relaxation to a minimum of the potential F, which, in the framework of the amplitude equations, corresponds to a linearly stable solution.Thus, the next step is to investigate the stability of the stationary solutions of Eqs.(13), which represent the isotropic state and the three possible patterns already introduced: The stripe pattern consisting of one mode, the square vortex lattice consisting of two and the triangular lattice consisting of three modes.
To this end, we determine the linear stability of each of these patterns with respect to small perturbations of the amplitudes and angles.As Eqs. ( 13) constitute a set of ordinary differential equations, we restrict ourselves to purely time-dependant perturbations that can be expressed in the reduced five-dimensional phase space consisting of A 1 , A 2 , A 3 , φ 1 and φ 2 .(More complex inhomogeneous perturbations will be later considered in Section IV.)The detailed calculations are shown in Appendix B. First, we find that the stripe pattern is unstable with respect to the growth of an additional mode perpendicular to the first.This already hints at the next result, which is that the square lattice in fact constitutes a stable solution of Eq. ( 13).In contrast, the hexagonal pattern is unstable against perturbations that change both amplitudes and angles.
Thus, in our predefined set of states represented by the superposition of modes characterized by the critical wavenumber k c , only the square lattice pattern is actually linearly stable for a > 0. Therefore, it represents the only minimum of F. This result of course relies on the validity of the representation of the velocity field v(x, t) in terms of up to three modes according to Eq. ( 9).It is not guaranteed that a higher number of modes does not also yield a linearly stable state and thus represents a minimum of F as well.However, the fact that already the three-mode, triangular state is linearly unstable indicates otherwise.A final answer would require an analysis going beyond the scope of this work.Meanwhile, numerically solving the full dynamic equation (1), we not only observe a square vortex lattice for λ = 0 [62], but also for small values of λ below the threshold to turbulence.In fact, the numerical calculations show that λ has no impact on the observed patterns (including the lattice constant) in this stationary regime whatsoever, see Appendix C.
To analytically explore the relevance of the square vortex lattice as given in Eq. ( 11) for the full nonlinear dynamics governed by Eq. ( 1), we perform a weakly nonlinear analysis, the details of which are given in Appendix C. To this end, we employ the standard procedure utilizing multiple scales [65,66] and define the small expansion parameter ε characterizing the distance to the linear instability of the uniform state using the growth rate of perturbations, i.e, ε 2 = σ(|k| = k c ) = a, see Eq. ( 4).The basis of the analysis is an expansion of the velocity field with respect to ε, i.e., close to the instability.The full expansion is given as where c.c. denotes the complex conjugate.Here, the indices m and n are used to introduce higher harmonics that may become increasingly relevant when moving further from the instability.For the sake of a simplified notation, we use A 1 = A y,1,1,0 and A 2 = A x,1,0,1 for the dominant modes of the square vortex lattice in the following.This mirrors the notation used above and in Appendix B. We have further introduced the slow time, T , and long spatial scale, X, of amplitude modulations.These are given as T = ε 2 t and X = εx (motivated by the fact that the growth rate scales in lowest order quadratic in k).Inserting the expansion ( 14) into Eq.( 1) and matching terms with the same order of ε (see Appendix C for details), we obtain evolution equations for the two dominating amplitudes A 1 and A 2 in O(ε 3 ), where we already scaled back to the fast scales.Compared to the amplitude equations that we derived above to investigate the stationary patterns selected by the functional F, Eqs. ( 15) additionally contain gradient terms.This is because we here consider also modulations of these amplitudes in space.Eqs. ( 15) are real-valued Ginzburg-Landau equations [67] and thus describe simple relaxation towards the solution given by |A 1 | = |A 2 | = a/5.Remarkably, λ has no impact on the amplitude equations.This result is in line with the full numerical solution of Eq. ( 1), where we also observe no impact of λ below the threshold to turbulence.Representing the square vortex lattice solely by two perpendicular modes [see Eq. (11)] is of course an approximation.In this context, the weakly nonlinear analysis enables us to calculate corrections to this approximation in the form of higher harmonics, which result from the nonlinearities in Eq. (1).However, as we discuss in detail in Appendix C, the amplitudes of these higher harmonics are orders of magnitude smaller than the amplitudes of the dominating modes.We conclude that we can safely continue our analysis using the two-mode approximation of the square vortex lattice.In the following, we will denote this lowest-order square pattern as v s (x), i.e., v s (x) = a 5 e iy e ix + c.c. , and employ it as an approximative solution of Eq. ( 1)

IV. EXTENDED STABILITY ANALYSIS
Having established the presence of the square vortex lattice for small values of λ, we now turn to the stability of this pattern with respect to an increase in λ.As the discussion above has shown, amplitude equations derived via standard procedures are only of limited use in this context.Instead, we will again employ a linear stability analysis of the full nonlinear equation (1).Taking into account the effect of the advection term is, however, not straightforward.This can be shown by naively considering space-and time-dependant perturbations denoted as δv and δq given by a single mode [analogously to Eq. ( 3)], i.e., where δv and δ q are the perturbation amplitudes, σ is the growth rate, k the wavevector of the perturbative mode and v s is given by Eq. ( 16).
Linearizing the evolution equation for v(x, t) [Eq.( 1)] around v s (x) yields where the Jacobian J excludes the term in Eq. ( 1) involving q.As we have specified the spatial and temporal functional form of the perturbation [see Eq. ( 17)], we are able to evaluate the spatial and temporal derivatives.Thus, after canceling out the factor e ik•x , Eq. ( 18) can be written in terms of the perturbation amplitudes δv and δ q as We introduce the matrix B(x) as where B 0 contains all terms independent of the spatial variable x and h.o.t.denotes higher-order terms, which still involve x, see Appendix D for the explicit components of B(x).To continue, we neglect these higher-order terms and obtain an equation independent of x, As a final step, we deal with the perturbation of the quantity q.To this end, we utilize the scalar product of k with Eq. ( 21).Using the incompressibility condition, in Fourier space given as k • δv = 0, we obtain Inserting Eq. ( 22) into Eq.( 21), we finally arrive at the familiar eigenvalue problem where the matrix M is calculated as Determining the eigenvalues of M, we obtain Compared to the growth rate of perturbations to the isotropic state [Eq.( 4)], the growth rate in Eq. ( 25) is again real-valued but contains the additional term −8A 2 .This term is a result of the nonlinear cubic term in Eq. ( 1) and describes a strong damping of perturbations dependent on the amplitude A of the square vortex lattice pattern that is already present.Importantly, Eq. ( 25) does not contain any contribution from the nonlinear advection term, i.e., is independent of λ.This is because all contributions to B that stem from nonlinear advection represent higher-order terms that are neglected in the analysis above, see Eqs. (D1) to (D4) in Appendix D. These terms result from the coupling between the perturbative mode and the modes constituting the vortex lattice.In order to consistently take these contributions into account we have to extend the linear stability analysis.In the following, we will outline our approach to achieve this.At this point, it is important to understand how modes in the system are coupled as a result of the mathematical form of the dynamic equation and the already present vortex lattice state.As this state is not uniform, a perturbative mode can excite or dampen other modes via the nonlinear terms even in the linearized case.This is clearly demonstrated by the components of the matrix B, which is introduced above and explicitely given in Appendix D. As was previously pointed out [68], a single-wavevector perturbation is thus not sufficient to capture this effect.We here use an ansatz that can not be represented by a single wavevector but rather consists of a more complex subset of Fourier space.To find an appropriate subset, it makes sense to visualize the wavevectors corresponding to the square vortex pattern in Fourier space, see Fig. 3(a).The primary modes that constitute this lattice, i.e., 2 cos(x) = e ix + e −ix and 2 cos(y) = e iy + e −iy are shown as blue squares in the figure.By examining the combination of the nonuniform base pattern and the kinds of nonlinear terms that are present, we determine the subset of wavevectors that captures the coupling between modes relevant for the linear stability analysis.In our case, due to the form of the square vortex lattice and the presence of the cubic and the advection term, we have to take into account wavevectors that are arranged on a grid in Fourier space with a spacing of k c = 1, oriented the same way as the vortex lattice modes.To illustrate the suitability of this grid of modes, let us consider one of the emerging terms as an example.After inserting the sum of the square vortex lattice and the perturbation given via Eq.( 17), we linearize the nonlinear advection term v • ∇v, which yields, among other contributions, a term ∝ v sx ∂ x δv x .Here, v sx and δv x are the x-components of the square vortex lattice and the perturbation, respectively.Evaluating this term gives a contribution ∝ ik x e ikxx+i(ky+1) , i.e., a mode that is shifted by ∆k = (0, k c ) = (0, 1) compared to the original mode of the perturbation k = (k x , k y ).When combining the non-uniform base state with a perturbative mode, all nonlinear terms lead to similar contributions characterized by such a shift.As a result, modes that are additionally excited or damped lie on a grid in Fourier space.Taking a closer look at the terms in the matrix B in Appendix D further demonstrates the relevance of this particular subset.
The ansatz for the perturbation is visualized in Fig. 3(a) via the yellow dots and can be written as The wavevectors k mn consituting the grid are expressed relative to the wavevector of the center mode k 0 = (k 0x , k 0y ) at m = 0 and n = 0 via Note that, for k 0 = 0, the ansatz (26) just describes vectors of the reciprocal lattice of the original square lattice in real space.In order to describe a perturbation from this pattern, k 0 must be non-zero.The total number of perturbative modes described in Eq. ( 26) is given by N p = (2M + 1) 2 .For the remainder of this paper, we choose M = 1, thus considering a perturbation that can be represented by 9 wavevectors, see Fig. 3(a).Now, inserting the perturbation [Eq.( 26)] into the dynamic equation ( 1), the nonlinear coupling between modes can be taken into account even after linearization.To continue with the stability analysis, we write down a system of 2N p linearized equations by collecting terms of O(e ikmn•x ) from the equations for v x and v y , respectively.Note that, again, the incompressibility condition, ∇ • v = 0, has to be satisfied by employing a projection operator as in Eq. (24).Finally, we obtain the eigenvalue problem where the vector δ V now contains the perturbation amplitudes of all modes in both x-and y-direction, i.e., By calculating the eigenvalues of the matrix M * , we now determine the linear stability of the square vortex lattice with respect to the simultaneous growth of multiple modes.The calculation outlined above is quite involved as the matrix M * has dimensions 2N p × 2N p .This necessitates the use of symbolic computation software such as the Python library Sympy [69].Writing down a closed form for the eigenvalues is also not feasible due to the high dimensionality.Rather, we evaluate the matrix M * for a specific set of parameters a, λ and k 0 and calculate the eigenvalues numerically.
The stability of the square lattice base pattern is then determined by the largest real part of any of these eigenvalues, which we denote as σ max in the following.Fig. 3

(b) to (d) shows the largest real part σ max
Re as a function of the grid offset k 0 at a = 0.5 and different values of λ.Note that the four-fold rotational symmetry is expected due to the symmetry of the square lattice pattern.Taking a closer look at Fig. 3(b), we observe that σ max Re stays negative for all k 0 at smaller advection strength λ.However, increasing the advection strength, σ max Re becomes positive for certain k 0 , see Fig. 3(c) and (d).In particular, the perturbations that start to grow first are characterized by k 0 located on the axes k 0 = (k 0x , 0) or k 0 = (0, k 0y ).Picking one of these axes, we plot σ max Re as a function of k 0x for k 0y = 0 in Fig. 3(e) to be able to visualize the dependency on λ in more detail.Here, we see that the growth rate becomes positive at some critical value λ c (2.0 < λ c < 2.8) and develops a clear maximum at k 0x ≈ 0.25 for larger values of λ.This picture stays qualitatively the same even for larger values of a.We always observe the development of clear maxima on the positive and negative axes [k 0 = (k 0x , 0) and k 0 = (0, k 0y )], respectively.
At this point, let us briefly discuss how the number of perturbative modes might affect the stability analysis.A subset of Fourier space consisting of a grid of 3 × 3 modes represents the lowest-order choice that allows for an investigation of instabilities caused by the nonlinear advection term, as we show in this work.In principle, the accuracy of the matrix M * in Eq. ( 28) and is plotted here as a function of k0.When increasing λ, the square lattice becomes unstable to perturbations with offset k0 on the axes k0 = (k0x, 0) and k0 = (0, k0y), as indicated by the emergence of red and yellow colors.To illustrate, the maximum growth rate is shown for k0y = 0 in (e) at different values of λ.
of the results might be increased by allowing for a larger grid of perturbative modes.However, these additional modes are characterized by larger and larger wavevectors.As a result, these will be increasingly damped by the biharmonic term in Eq. ( 1) and thus only have a diminishing influence.Further, when aiming for a consistent stability analysis including such higher-order perturbations, the representation of the square lattice as given in Eq. ( 16) should also be amended by the inclusion of higher-order modes.Note that an infinite grid of perturbative modes implies a translation symmetry in Fourier space for wavevector shifts of ∆k 0x = 1 or ∆k 0y = 1.Thus, it is sufficient to explore the region given by −k c /2 < k 0x < k c /2 and −k c /2 < k 0y < k c /2, as we have done in Fig. 3.As a perturbation consisting of 3 × 3 modes already leads to quite accurate results concerning the stability as compared to the numerical solution of the full Eq.( 1), we restrict our study to this case.In order to corroborate our results, we have additionally performed a limited analysis using 5 × 5 perturbative modes, where we find qualitatively the same results.The quantitative differences are small, in particular close to the onset of the instability.
So far we have concentrated on the real part of the eigenvalues, as this gives access to the growth rate.Appendix E further discusses the imaginary part of the eigenvalue with largest real part σ max .We find that σ max Im vanishes for perturbations with k 0 located on the axes k 0 = (k 0x , 0) or k 0 = (0, k 0y ), but is nonzero in large regions of k 0 -space.In particular, σ max Im ̸ = 0 for most unstable perturbations close to the fastest-growing perturbation, indicating the emergence of oscillatory behavior or even more complex dynamic states.The background color is used to show numerical results: dark blue for a stationary vortex lattice and light red for mesoscale turbulence.These results are obtained by numerically solving the evolution equation for v(x, t) directly, starting from a square vortex lattice with small perturbations and observing if the lattice becomes unstable and turbulence develops.The dashed red line shows the results from the linear stability analysis, i.e., the point where the real part of any of the eigenvalues of the matrix M * in Eq. ( 28) becomes positive.The area shaded in green shows the region where both stable vortex lattice and turbulent state can be encountered depending on the initial state of the system.We thus observe hysteretic behavior.This region is determined via additional numerical simulations discussed in Section VI.

V. ONSET OF TURBULENCE
The instability of the square vortex lattice with respect to the simultaneous growth of multiple modes signifies the start of the development of mesoscale turbulence in our model.Having established the procedure of our extended stability analysis, we now turn to the determination of the critical point where perturbative modes start to grow.To this end, we calculate all eigenvalues of the matrix M * as outlined above for a set of parameters a, λ and k 0 .If all of the eigenvalues for all k 0 have a negative real part, the square vortex lattice is stable.However, if any of the eigenvalues for any k 0 has a positive real part, it is unstable and mesoscale turbulence is expected to develop via the mutual excitement of multiple modes.On this basis we can determine the critical advection strength λ c as a function of the coefficient a.
Figure 4 depicts the state diagram in a-λ-space.The dashed red line shows the result obtained via the extended linear stability analysis.The background color indicates results obtained via direct numerical solution of Eq. ( 1) for specific sets of parameters λ and a.Here, we start the numerical calculations from an already ordered state and observe if initially small perturbations start to grow.In this way, we determine whether the system develops a turbulent state (indicated by the light blue color) or the lattice structure remains stable (indicated by the dark blue color).Remarkably, the extended linear stability analysis predicts the onset of the turbulent state quite accurately.Figure 4 demonstrates that the vortex lattice is stable for much higher values of λ when a is small.This observation is explained by the difference in the stationary amplitude A s of the modes of the vortex lattice, which is given by A s = a/5, see Section III.Increasing a leads to a larger amplitude.The nonlinear nature of the advection term implies that its impact grows with increasing amplitude of the patterns, thus explaining the earlier onset of turbulence for larger a.
Up until this point, we have only looked at the eigenvalues of the matrix M * .However, the corresponding eigenvectors also contain valuable information.In particular, we can determine which of the modes that constitute the perturbation will start growing the fastest once the instability sets in.To this end, we take a closer look at the eigenvector δV c corresponding to the eigenvalue with largest real part.This eigenvector contains the growth amplitudes of all N p modes of the perturbation in terms of components δv mn,x and δv mn,y .Calculating the absolute value of the components of δV c corresponding to a particular mode, we determine the growth of that mode once λ is increased above λ c .
As an illustration, the absolute eigenvector components are plotted in Fig. 5(a) as the saturation of the red coloring at the respective positions in the grid of perturbative modes in Fourier space.The coefficients are set to a = 0.5 and λ = 3.16 and we have chosen k 0 = (0.35, 0), which corresponds to the maximum of the growth rate for this parameter set.As the red color indicates, we find that six modes grow particularly quickly once the instability sets in.This result is easily verified by comparison with Fourier spectra obtained via direct numerical solution of Eq. (1).To this end, we initialize the numerical calculations with a square lattice pattern according to Eq. ( 16) and add very small random perturbations not favoring any particular wavevector.Taking a closer look at the evolution of the Fourier spectrum during the onset of the instability, we observe which of the wavevectors start growing first.The Fourier spectrum at such an intermediate time step is shown in Fig. 5(b).There, we plot the absolute value of the vorticity transformed to Fourier space, |ω|(k) for the same values of the coefficients used in Fig. 5(a).The reciprocal lattice corresponding to the square vortex lattice is clearly visible (indicated by the square boxes).The plot also shows which other modes are starting to grow at this point.For comparison, we have overlaid the spectrum with the grid of the perturbative modes shown in Fig. 5(a) as a frame indicated by the circles.The spectrum agrees very well with the analysis based on the eigenvector δV c , corroborating the results obtained via the extended linear stability analysis.Note that the four-fold rotational symmetry in Fourier space is clearly visible in the spectrum, which was already discussed in the context of Fig. 3(b) to (d).Thus, in addition to the six modes determined in Fig. 5(a), there are 18 more modes that are equally relevant at the onset of the instability, yielding a total of 24 modes.To visualize these elaborate perturbations, we transform back to real space in Fig. 5(d).Here, we show the vorticity field of the square vortex lattice on the left-hand side.Moving to the right, we add the perturbative modes according to the eigenvector The blue line denotes the square vortex lattice branch, which becomes unstable at some critical value of λ, as discussed in Section IV.Here, the system develops turbulence and vrms drops significantly.
The red data points are determined via numerical solution of the full Eq.( 1) by starting in the turbulent regime and slowly decreasing λ.(b) Root-mean-square velocity vrms as a function of a at λ = 2.1 for different spatiotemporal states.The isotropic state is denoted by the green line, the square vortex lattice by the blue line.Increasing a leads to growing amplitude and thus also growing vrms, until the lattice becomes unstable.Again, the turbulence branch of the hysteresis loop is obtained via numerical calculations by starting in the turbulent regime and slowly decreasing a.In both plots, the shaded area denotes the region where both stable vortex lattice and turbulent state can be encountered.
shown in Fig. 5(a) (plus the grid of perturbative modes rotated by 90 • , 180 • and 270 • ).Increasing the amplitude of the perturbation, we can thus visualize how the square vortex lattice becomes distorted right after the onset of the instability leading to the turbulent state, thereby mimicking the time evolution occurring in the full system.This is further visualized in the Supplementary Movie, where we show the evolution of both Fourier spectrum and vorticity field obtained via numerical solution of Eq. ( 1) at the onset of the instability.Fig. 5(c) further includes the Fourier spectrum observed at a time step in the fully developed mesoscale-turbulent state.Here, the dominant modes form a broad ring with radius |k| ≈ 1 and, as expected, the system exhibits complete rotational symmetry.

VI. HYSTERETIC BEHAVIOR
Up until now, the discussion has focused on the stability of the square vortex lattice.The instability we uncovered indicates the critical advection strength λ c (a), above which turbulence starts to develop from an initial square vortex lattice.In particular, this gives the expected behavior when λ or a are slowly increased from initially small values.Here, when numerically solving the full nonlinear equation ( 1), we observe a discontinuous jump from the stationary lattice to dynamic turbulent behavior at λ c or a c .The discontinuous nature of the transition already hints at the coexistence of different solutions and the possibility of hysteresis, i.e., different behavior when λ or a are instead slowly decreased from initially large values.
Investigating this issue within an analytical framework is not feasible due to the complexity of the turbulent state.Thus, we continue with a numerical study.Here, we solve the full nonlinear equation (1) for a particular value of λ and wait until we obtain a statistically stationary state.This state then acts as the starting point for a subsequent calculation with a slightly smaller value of λ.We continue in this way until the flow field becomes stationary and turbulence is superseded by a stationary vortex lattice.To quantify the results, we here employ the root-mean-square velocity v rms = ⟨|v| 2 ⟩, where the mean value ⟨. . .⟩ is obtained in the statistically stationary state as both a spatial and temporal average.Note that the average velocity ⟨v⟩ vanishes for both square vortex lattice and turbulence and, thus, is not a suitable measure to distinguish between different spatiotemporal states.In the stationary vortex lattice, v rms can be calculated analytically, which yields v rms = 2A s .In the turbulent state, v rms is generally smaller because the nonlinear advection term leads to energy transfer between scales and subsequent dissipation, which comes on top of the amplitude saturation via the cubic term in Eq. ( 1) [58,62].Fig. 6(a) shows the root-mean-square velocity as a function of advection strength at a = 0.5.The root-mean square velocity of the vortex lattice is simply given by a constant value in this diagram because the stationary amplitude solely depends on a.As we determined via the extended stability analysis in Section IV, at λ = λ c the vortex lattice becomes unstable to the simultaneous growth of multiple modes.In Fig. 6(a), this is indicated by the transition to the dashed line.Here, a turbulent state emerges, which manifests itself as discontinuous drop in the value of v rms .As discussed above, the turbulent branch has to be determined numerically via solution of Eq. ( 1).In Fig. 6, the results are shown by the red data points.Slowly decreasing λ, we move along the turbulent branch from the right-hand side to the left-hand side until v rms jumps to the vortex lattice branch at λ ⋆ , at which point the system has settled into a stationary state.As Fig. 6(a) shows, the corresponding value of λ is significantly smaller than the point where the square vortex lattice becomes unstable.The region in between is characterized by the existence of both stable regular vortex lattice and turbulent state.As a result, we observe hysteretic behavior.
To further illustrate the hysteresis loop, we plot v rms for constant λ as a function of a in Fig. 6(b).This diagram also gives a complete overview of the various points of interest and motivates the representative sketch in Fig. 1.First, the isotropic solution (v = 0), which is indicated by the green line, becomes unstable at a = 0.The subsequent unstable branch is shown as the dashed green line.From this instability emerges the square vortex lattice (blue line), which becomes unstable at a c , resulting in a turbulent state, compare Fig 4 .Again, the red data points show the results obtained via direct numerical solution of Eq. ( 1), starting in the turbulent regime and slowly decreasing a. Performing additional numerical calculations for different values of a, we determine the hysteresis region in a-λ-space, which is shown as the area shaded in green in Fig. 4. We obtained similar results for the lower limit of the hysteresis from numerical calculations starting from an initially isotropic state with small random perturbations [62].Due to the random initial values, the development of a disordered, turbulent state is favored.However, when the activity parameters are small enough, the system eventually settled into a regular square vortex lattice.In this way, the hysteresis region shown in Fig. 6 is reproduced independently for different initial conditions.Further, a hysteretic transition was also found in a recent investigation of the dynamics in the Toner-Tu-Swift-Hohenberg model under circular confinement [49].In that case, the transition is observed between regularly oscillating vortices and chaotic dynamics.
Since the transition from turbulence to stationary lattice is not associated with a crossing of an eigenvalue of the spectra obtained in the respective linear stability analysis above, we hypothesize that there might be a connection to a statistical phase transition analogous to classical turbulence systems like pipe flow [13].Alternatively, a long-living spatiotemporally chaotic transient [70] is possible, but can be largely ruled out due to the fact that the hysteresis region is reproduced in simulations with random initial conditions as mentioned above.This conclusion is further corroborated by long-running numerical simulations, which are presented in Appendix A and show no sign of transient behavior.Comparing Fig. 6(a) and (b), we observe that the jump in v rms at the lower end of the hysteresis region decreases from ∆v rms ≈ 0.04 at a ⋆ = 0.5 and λ ⋆ ≈ 1.48 to ∆v rms ≈ 0.02 at a ⋆ ≈ 0.13 and λ ⋆ = 2.1.Additional simulations used to determine the hysteresis region shown in Fig. 4 confirm the observation that the jump diminishes when moving to the left of the state diagram (decreasing a and increasing λ) and might even vanish altogether.More detailed numerical studies will be necessary to give a complete characterization of the transition at the end of the hysteretic branch at low activity.

VII. STABILITY OF REGULAR VORTEX LATTICES WITH ARBITRARY WAVENUMBER
Up until this point, we only considered vortex lattices comprised of two modes with wavenumber k = k c , which is the only possibility exactly at the onset of the instability at a = 0.However, increasing a, vortex patterns consisting of modes with wavenumbers k ̸ = k c might be stable as well.The extended stability analysis proposed and discussed in Section IV may also be utilized to investigate these more general structures.Motivated by the results on stationary pattern selection in Section III, we will restrict ourselves to square vortex lattices consisting of two perpendicular modes with wavevectors |k 1 | = |k 2 | = k s .Similarly to the procedure employed in Section III, evolution equations for the amplitudes of the modes, A 1 and A 2 , can be derived.Again neglecting spatial modulations and higher harmonics yields Compared to Eqs. ( 13) and ( 15), the linear coefficient now depends on the wavenumber of the lattice and is given by α(k s ) = a − 1 + 2k 2 s − k 4 s .The stationary solution is determined via The light red area represents the band of unstable modes that might grow from the isoptropic state according to Eq. ( 5).The dashed line shows the critical wavenumber kc = 1.In compliance with the emergence of turbulence, we observe that the stability balloons rapidly shrink when increasing λ.
velocity field, the square vortex lattice in this two-mode approximation is thus given for arbitrary lattice wavenumber k s as e iksy e iksx + c.c. .
Setting k s = k c = 1 reproduces the vortex lattice investigated before, see Eq. ( 16).
To determine the stability of the vortex lattice of arbitrary wavenumber, we again add a perturbation consisting of multiple modes.As before, we take into account 3 × 3 modes that are located on a grid in Fourier space, albeit with a grid spacing of k s instead of k c [compare Eq. ( 27)].Performing the same procedure as introduced in Section IV, we then calculate the maximum growth rate now depending on not only a and λ but also k s .Observing if the growth rate becomes positive, we determine whether the square vortex lattice is unstable with respect to the simultaneous growth of multiple modes.The results can be illustrated using the concept of stability balloons [71,72].A few examples are plotted in Fig. 7, where we show the regions of stable patterns in the space spanned by k s and a for various values of λ.The light red region in these plots is limited by Eq. ( 5).Here, perturbations of the isotropic, uniform state will grow, but no regular patterns are stable.As expected, we find that the region of stable patterns is largest for vanishing advection strength, λ = 0. Increasing λ, the stability balloon shrinks, which complies with the emergence of the turbulent state.Note that vortex lattices characterized by the critical wavenumber, k s = k c , are located on the dashed line in Fig. 7. Remarkably, these lattice are not necessarily the most stable.We rather find that lattices with slightly increased k s remain stable for larger values of a.
The stability of regular structures in other systems, e.g., stripe-like structures in Rayleigh-Bénard convection [72], is limited by different types of instabilities depending on the particular boundary of the stability balloon.With that in mind, Appendix F explores the specific instabilities at the boundaries for smaller and larger square lattice constant k s at λ = 0, i.e., to the left and the right of the stability balloon shown in Fig. 7(a).Here, we indeed observe very different instabilities at the left and right boundary.However, both are dominated by the growth of modes with k ≈ k c , thus leading to structures closer to the ones that emerge when starting from an isotropic state.These results indicate that the instabilities of square lattices with arbitrary wavenumber are similarly complex and diverse as the situation in other systems, such as Rayleigh-Bénard convection [72].

VIII. CONCLUSIONS
In summary, this work investigates the onset of the turbulent state in an established minimal model for active fluids, where the evolution of the effective velocity is determined as a competition between gradient dynamics governed by a functional F and nonlinear advection.We first establish that a square lattice of vortices represents a minimum of F and thus a stable solution in the absence of nonlinear advection.In analogy to high-Reynolds-number turbulence, nonlinear advection destabilizes the stationary pattern and induces a turbulent state.Here, we show, by means of an extended stability analysis that the transition follows from a linear instability.The approach is based on linearization around an approximation of the analytical vortex lattice solution.Analogous results are found if the polar alignment strength is increased above a critical value.Since both parameters increase with activity, we can state that strong enough activity is a necessary ingredient for active turbulence in the model for polar active fluids studied here and in related experiments.Remarkably, the instability is only revealed when taking into account the mutual excitement and simultaneous growth of multiple perturbative modes.Moreover, our results show that hysteretic behavior is a characteristic feature tightly connected to the onset and disappearance of mesoscale turbulence.The phenomenology observed in the hysteretic transition region between stable vortex lattices and active turbulence is reminiscent of the picture of a chaotic saddle separating laminar and turbulent regimes [73] and of a saddle-node bifurcation of coherent structures [74] to the first occurence of turbulence from the laminar state in classical pipe flow [13,75].
Our work complements earlier studies that include stability analyses of nonuniform states, e.g., a similar approach was recently used to determine the stability of traveling states in active phase field crystals [76].Nonuniform states are also sometimes investigated in an equilibrium situation, as in microphase-separated diblock copolymer melts [77], which allows for methods not available here due to the non-equilibrium nature of active fluids.
Regarding the experimental relevance of the present study, it is important to note that square and triangular lattices, see Fig. 2(b) and (c), can be stabilized in bacterial suspensions by placing appropriate arrangements of small obstacles in the flow field [39].In these experiments [39], it was also shown that the square vortex lattice is more robust against changes of the arrangement than the triangular lattice.This observation is in line with our analytical results, i.e., that only the square vortex lattice is in fact an intrinsically stable pattern in the model.The stability analysis suggests that, in principle, square vortex lattices may be achieved in bacterial suspensions without external stabilization if sufficient control of the activity is possible (e.g., by varying the availability of oxygen [53]).Further, analytical studies of a similar model for mesoscale turbulence have shown that a shear-thickening solvent fluid may induce regular vortex lattices as well [78].Note that structures similar to the earlier predicted hexagonal vortex lattice [37,58] have recently been observed in experiments [79].
In the context of stabilized vortex lattices, the same model that we investigated in this study was employed in Refs.[40,41] to explore the dynamics of microswimmer suspensions within arrangements of small obstacles.These obstacles are able to stabilize regular vortex patterns even for values of the nonlinear advection strength λ much higher than the critical value λ c determined in the present work.Interestingly, the transition between a stabilized square vortex lattice and the disordered state facilitated by increasing λ showed features of a non-equilibrium phase transition in the universality class of the two-dimenstional Ising model [41].This striking contrast to the linear instability observed in the present work in the absence of obstacles demonstrates that the route to turbulence in active polar fluids is strongly affected by the geometry of the system and by confinement, as is the transition to inertial turbulence [1,2].
In active nematics the onset of the turbulent state was recently linked to a directed percolation transition [25], similar to the route to sustained inertial turbulence in channel, Couette and pipe flows [6,8,13].In fact, Ref. [25] also investigated active nematic turbulence in a channel geometry that emerged from increasing the channel width and hence decreased the influence of spatial confinement.This is in contrast to the bulk system without spatial constraints investigated in the present study.Clearly, further work is necessary to establish whether the different routes to turbulence in these active systems are a consequence of different types of setup or possibly due to the relevance of different order parameters in these systems.shown in Fig. 5, we use a higher spatial resolution of 400 × 400 grid points and set the system size to 80π × 80π.Further, the snapshots shown in Fig. 1 are extracted from a system with 256 × 256 grid points and a size of 32π × 32π.
In the context of the calculations to explore hysteretic behavior (see Fig. 6), we set the number of grid points to 256 × 256 and the system size to 48π × 48π.To explore the hysteresis region we start from a turbulent state at high activity and then successively decrease one of the activity parameters, λ or a.The data points in Fig. 6(a) are determined by decreasing λ by ∆λ = 0.05 every ∆t = 10000 time units and then measuring v rms at the end of the respective time periods.It is possible that the persistent turbulent states observed in the hysteretic region are longliving spatiotemporally chaotic transients.In order to investigate this possibility, we further performed simulations with very large running times where we set the time periods after every successive step of decreasing the activity to ∆t = 50000.Here, we use a = 0.5 and decrease λ in steps of ∆λ = 0.02.Fig. 8 shows the evolution of the root-mean-square velocity for the last few steps before the system develops a stationary state.Using these very large time periods ∆t and smaller successive steps ∆λ, we find that the system settles into the vortex lattice at same value λ * ≈ 1.48 as observed when using smaller ∆t and larger ∆λ, compare Fig. 6(a).Thus, there are no signs of a transient phenomenon.The data points in Fig. 6(b) are determined by decreasing a by ∆a = 0.02 every ∆t = 50000 time units.Finally, in order to determine the hysteresis region in the state diagram in Fig. 4, we use ∆t = 50000 as well as steps of ∆a = 0.02 or ∆λ = 0.05.
Concerning the extended linear stability analysis described in Section IV, we determine the matrix M * with the help of the Python library Sympy [69].Subsequently, inserting specific values for the coefficients a and λ, the wavevector of the center of the grid of perturbative modes, k 0 , and, if applicable, the wavenumber of the square lattice, k s , we evaluate the eigenvalues of M * numerically.When determining the critical advection strength, λ c , in the context of the state diagram in Fig. 4, we regard the square vortex lattice as unstable when the largest eigenvalue is larger than σ max > 10 −3 .This is done in order to avoid artifacts that may result from neglecting higher harmonics.The same approach is applied to calculate the stability balloons in Fig. 7. instead of the full functional F [see Eq. ( 2)].First, we utilize the dynamic equation governing the velocity field v(x, t) [Eq.( 1)] and write which determines the amplitude equations for amplitudes A 1 , A 2 , A 3 as well as the evolution equations for the angles φ 2 and φ 3 .
The prefactors in Eq. (B3) are quadratic in the partial derivatives of v with respect to A 1 , A 2 , A 3 , φ 2 and φ 3 .The same argument as for f applies, see Eq. ( 12), and the prefactors are sums of a homogeneous term and higher order modes e 2ilk•x with integer l, i.e., Again, we neglect higher order contributions and insert Eq. (B4) into Eq.(B3).As a final step we evaluate the derivatives of f 0 with respect to A 1 , A 2 , A 3 , φ 2 and φ 3 and arrive at a closed set of ordinary differential equations for the amplitudes and angles, see Eqs. (13).
Stationary solutions of Eqs. ( 13) represent the regular patterns introduced in the main text, i.e., the isotropic solution, stripe patterns, as well as square and hexagonal lattice-like structures.In the following, we determine the linear stability of these solutions.First, for the sake of a simplified notation, we define the vector c = (A 1 , A 2 , A 3 , φ 2 , φ 3 ), which contains the phase space variables, and write the set of Eqs.(13) as where G = (G A1 , G A2 , G A3 , G φ2 , G φ3 ) contains the right-hand sides.Stationary solutions are then denoted by c 0 .In order to determine the stability, we add a small perturbation to c 0 , where δĉ = (δ Â1 , δ Â2 , δ Â3 , δ φ2 , δ φ3 ) are the perturbation amplitudes and σ denotes the growth rate.Inserting Eq. (B6) into Eq.(B5) and linearizing yields where J| c0 denotes the Jacobian J evaluated at c 0 .If J| c0 possesses eigenvalues with positive real part, the stationary solution c 0 is linearly unstable.Furthermore, the eigenvectors corresponding to the positive eigenvalues indicate the kind of perturbations expected to grow.Isotropic state.The trivial solution of Eq. (B5), A 1 = A 2 = A 3 = 0, represents the isotropic state, i.e., v = 0 everywhere.Here, the two angles φ 2 and φ 3 are kept as free parameters for the linear stability analysis.We write down the Jacobian J| c0 for c 0 = (0, 0, 0, φ 2 , φ 3 ) and calculate the eigenvalues, which yields with corresponding eigenvectors δĉ 1 = (1, 0, 0, 0, 0) , δĉ 2 = (0, 1, 0, 0, 0) , δĉ 3 = (0, 0, 1, 0, 0) , δĉ 4 = (0, 0, 0, 1, 0) , δĉ 5 = (0, 0, 0, 0, 1) .(B9) The first three eigenvalues are negative for a < 0, which means the isotropic solution is stable and represents a minimum of F. For a > 0, however, all three eigenvalues become positive, the isotropic solution becomes linearly unstable with respect to perturbations characterized by the eigenvectors corresponding to σ 1 , σ 2 and σ 3 .These eigenvectors describe growth of one of the amplitudes of the modes A 1 , A 2 and A 3 , respectively.The remaining eigenvalues σ 4 , σ 5 vanish regardless of a.This means that the orientation of the modes is irrelevant and the amplitudes start to grow for any φ 2 and φ 3 , provided a > 0. This immediately raises the question which of the possible configurations of the modes represented via Eq.( 9) will be selected in the emergent patterns, i.e., which of the stationary solutions of Eq. (B5) is linearly stable.We will investigate this in the following, starting with a state containing a single mode, i.e., a stripe pattern.FIG. 9. Stability of the stripe pattern consisting of one mode with angle φ1 = 0 and amplitude A1 = a/(3b).Shown are the growth rates σ2 and σ4 connected to mode i = 2 as a function of the angle φ2, where σ2 is connected to a growth of the amplitude A2 and σ4 is connected to a change of the angle φ2.The stripe pattern is unstable with respect to modes that are approximately perpendicular, i.e, φ2 ≈ π/2.The growth rates connected to the third mode (σ3, σ5) behave analogously with respect to the angle φ3.
Stripe pattern.For a > 0, one of the stationary solutions of Eqs. ( 13) is characterized by two of the amplitudes being zero.The other, non-zero mode yields a stripe pattern in the vorticity field, see Fig. 2(a).Without loss of generality, we set A 2 = 0 and A 3 = 0 and find A 1 = a/3.The angles φ 2 , φ 3 are undetermined for this solution and will be considered as free parameters.Writing down the Jacobian J| c0 for c 0 = ( a/3, 0, 0, φ 2 , φ 3 ) and calculating the eigenvalues yields Depending on the angles φ 2 and φ 3 , the eigenvalues σ 2 , σ 3 , σ 4 and σ 5 can be positive or negative.Taking a closer look at the eigenvectors δĉ 2 , δĉ 3 , δĉ 4 and δĉ 5 corresponding to these eigenvalues, δĉ 2 = (0, 1, 0, 0, 0) , δĉ 3 = (0, 0, 1, 0, 0) , δĉ 4 = (0, 0, 0, 1, 0) , δĉ 5 = (0, 0, 0, 0, 1) , (B11) we find that the perturbations connected to σ 2 and σ 3 describe growing amplitudes A 2 and A 3 , respectively, whereas the perturbations connected to σ 3 and σ 4 describe changing the angles φ 2 and φ 3 .Looking at a plot of the growth rates σ 2 and σ 4 as a function of the angle φ 2 , see Fig. 9, the behavior becomes clearer.For angles close to π/2, i.e., perpendicular to the first mode, σ 2 is positive, which means that the amplitude A 2 will grow.In contrast, for angles close to 0 or π, σ 2 is negative, which means that perturbations described by modes that are not approximately perpendicular will not grow.Simultaneously, σ 4 is positive for these damped modes, which means that their orientations start to rotate.The behavior of the third mode is analogous.To sum up, the stripe pattern is unstable with respect to additional modes growing, with perpendicular modes exhibiting the highest growth rate, already hinting at the next solution we will investigate: the square lattice pattern.Square lattice pattern.The square lattice pattern is characterized by two modes with the same amplitude oriented perpendicular to each other, see Fig. 2(b).Without loss of generality, we set A 1 = A 2 , A 3 = 0, φ 2 = π/2 and keep φ 3 as a free parameter.Solving Eq. ( 13) for these assumptions yields A 1 = A 2 = a/5 for the stationary value of the amplitudes.The solution under investigation is then given by c 0 = ( a/5, a/5, 0, π/2, φ 3 ).Writing down the Jacobian J| c0 and calculating the eigenvalues we find All except one eigenvalue are negative, provided a > 0. The remaining eigenvalue σ 1 = 0 warrants a closer look.The corresponding eigenvector is δĉ 1 = (0, 0, 0, 0, 1), which describes a change of the angle φ 3 of the third mode.This mode, however, vanishes for a square lattice pattern and there is also no unstable perturbation that leads to a growth of A 3 .Thus, the eigenvalue σ 1 = 0 can be disregarded when determining the stability.We conclude that the square lattice pattern is linearly stable and therefore represents a minimum of the functional F. imaginary part of the corresponding eigenvalue.To this end, we first determine the eigenvalue σ max with the largest real part and then take a closer look at the imaginary part σ max Im .We plot the absolute value of the imaginary part in Fig. 11 as a function of k 0 at the same values of λ that are used in Fig. 3(b)/(c)/(d) where the real part is shown.The imaginary part is nonzero for large regions of k 0 -space but vanishes on the axes k 0 = (k 0x , 0) and k 0 = (0, k 0y ).Recall that the maximum growth rate is observed for perturbations with k 0 located on one of these axes, see Fig. 3. Solely looking at the perturbation associated with the largest growth rate would thus indicate a stationary instability not directly connected to oscillatory or more complex dynamic behavior.However, above the onset of the instability (λ > λ c ), there is a set of unstable perturbations close to the fastest-growing perturbation, as is clearly visible in Fig. 3(c).As Fig. 11 shows, these perturbations with k 0 not on either axis [k 0 = (k 0x , 0) or k 0 = (0, k 0y )] are actually characterized by σ max Im ̸ = 0. Thus, the unstable perturbations close to the fastest-growing perturbation are instead associated with oscillatory behavior and the emergence of dynamic states.Further research is necessary to determine whether additional instabilities are involved in the emergence of fully turbulent states.exemplary values, k s = 0.87 and k s = 1.13, to explore the instabilities determining the borders of the stability region.We further set a = 0.2 and λ = 0. Again, we consider perturbations consisting of 3 × 3 modes on a grid in Fourier space as discussed in Section IV and VII.We determine the eigenvalues of the matrix M ⋆ , which are found to be reel-valued for λ = 0.The growth rate of perturbations is thus directly given by the largest eigenvalue.Fig. 12(a) and (d) show this growth rate for square lattices with k s = 0.87 and k s = 1.13.As is clearly visible, it obtains a maximum for perturbations characterized by a grid offset of k 0 = (0.42, 0) and k 0 = (0.18, 0.18) for the two lattice constants, respectively.The eigenvectors corresponding to these maxima are visualized in Fig. 12(b) and (e), where we show the perturbative modes as a grid in Fourier space.The red color indicates the relative importance of the particular modes, determined by the absolute value of the corresponding eigenvector components.The blue squares show the reciprocal lattices corresponding to the respective square lattice structure.The dashed circle indicates |k| = k c = 1, i.e., the critical wavenumber of instabilities of the isotropic state.We find that the growing modes of the perturbation are located on the circle, thus leading to emerging structures closer to k c .Although this is the case for both lattice constants, k s = 0.87 and k s = 1.13, the instabilities are indeed very different due to their orientation with respect to the square lattice.This difference becomes much clearer when visualizing the perturbations in real space.To this end, we plot the perturbed vorticity field in Fig. 12(c) and (f) for the two lattice constants.Here, a moderate perturbation according to the eigenvectors shown in (b) and (e) is added to the vorticity field in the respective square vortex lattice.For symmetry reasons, we also include the perturbative modes rotated in 90 • steps.The occurrence of different types of secondary instabilities at the boundaries of the stability region of particular patterns is also encountered in other systems.For example, the stability of stripe-like states in Rayleigh-Bénard convection is limited by the zigzag instability for smaller and the Eckhaus instability for larger wavenumbers [72].

FIG. 1 .
FIG. 1. Schematic overview of spatiotemporal states and their occurrence in the present model [Eq.(1)].(a) We here employ a as the control parameter and distinguish between different spatiotemporal states using the root-mean-square velocity observed in the system.Solid lines represent stable branches, dashed lines indicate unstable branches.Increasing a, the isotropic state becomes unstable against the emergence of a regular square vortex lattice at a = 0.This regular pattern, in turn becomes unstable at ac and the system develops turbulence.Following the turbulence branch while decreasing a, we observe a hysteresis loop as indicated by the arrows.(b)/(c) Representative snapshots of the vorticity field ω(x, t) in the square vortex lattice and in the mesoscale-turbulent state at a = 0.5 and λ = 3.0 obtained via numerical solution of Eq. (1).The superimposed arrows show the velocity field v(x, t) and the size of the snapshots is 12π × 6π.

FIG. 2 .
FIG. 2. Possible stationary patterns.The vorticity field ω(x) is visualized as color map.The superimposed arrows show the velocity field v(x).(a) Stripe pattern composed of a single mode in x-direction.(b) Square vortex lattice according to Eq. (11), which consists of clockwise-and counterclockwise-rotating vortices in an antiferromagnetic arrangement.(c) Triangular vortex pattern consisting of three modes.(d) Hexagonal vortex lattice with Kagome-like symmetry.Both vorticity and velocity are rescaled for visualization purposes.The box size for all patterns is 8π × 8π.

FIG. 3 .
FIG. 3. (a) Visualization of a perturbation consisting of multiple modes in Fourier space.The blue squares show the base vectors of the reciprocal lattice of the square lattice pattern, whereas the yellow circles show the wavevectors of perturbative modes that are located on a grid of size 3 × 3. The center of this grid is shifted by k0 from the origin.(b)/(c)/(d) Maximum growth rate of perturbations to the square lattice pattern at a = 0.1 and (b) λ = 2, (c) λ = 4 and (d) λ = 4.8.The growth rate is given by the largest real part of the eigenvalues σ max Re

FIG. 4 .
FIG.4.State diagram in the space spanned by the coefficient a and advection strength λ.The background color is used to show numerical results: dark blue for a stationary vortex lattice and light red for mesoscale turbulence.These results are obtained by numerically solving the evolution equation for v(x, t) directly, starting from a square vortex lattice with small perturbations and observing if the lattice becomes unstable and turbulence develops.The dashed red line shows the results from the linear stability analysis, i.e., the point where the real part of any of the eigenvalues of the matrix M * in Eq. (28) becomes positive.The area shaded in green shows the region where both stable vortex lattice and turbulent state can be encountered depending on the initial state of the system.We thus observe hysteretic behavior.This region is determined via additional numerical simulations discussed in Section VI.

FIG. 5 .
FIG. 5. Comparison of the eigenvector of the fastest-growing perturbation with the Fourier spectrum observed in the numerical simulations.The relative importance of eigenvector components is shown in (a) for the eigenvector corresponding to the eigenvalue with largest real part at a = 0.5 and λ = 3.16 and k0 = (0.35, 0).The saturation of the red color is determined by the absolute value of the eigenvalue components that correspond to the growth of the modes located at the respective positions on the grid and, thus, shows which of the modes are expected to grow the fastest when the instability sets in.The blue squares show the reciprocal lattice corresponding to the square lattice pattern.For comparison, (b) and (c) show the Fourier spectrum observed when numerically solving Eq. (1) at a = 0.5 and λ = 3.16.Here, (b) is obtained at an intermediate time where the initial square lattice pattern is still present (indicated by the square boxes) but the growth of other modes is already visible in the spectrum (indicated by the circles).Comparing with (a), we find good agreement with the analysis based on the eigenvector of the fastest-growing perturbation.The Fourier spectrum of fully developed mesoscale turbulence is shown in (c).(d) "Time evolution" (snapshots) of the vorticity field starting from the square vortex lattice (left) and showing increasingly perturbed patterns to the right.The form of the perturbation is given by the eigenvector shown in (a).The size of the snapshots is 12π × 12π.

FIG. 6 .
FIG.6.Hysteretic behavior.Solid lines show stable branches, while dashed lines represent unstable branches.(a) Rootmean-square velocity vrms as a function of λ at a = 0.5.The blue line denotes the square vortex lattice branch, which becomes unstable at some critical value of λ, as discussed in Section IV.Here, the system develops turbulence and vrms drops significantly.The red data points are determined via numerical solution of the full Eq.(1) by starting in the turbulent regime and slowly decreasing λ.(b) Root-mean-square velocity vrms as a function of a at λ = 2.1 for different spatiotemporal states.The isotropic state is denoted by the green line, the square vortex lattice by the blue line.Increasing a leads to growing amplitude and thus also growing vrms, until the lattice becomes unstable.Again, the turbulence branch of the hysteresis loop is obtained via numerical calculations by starting in the turbulent regime and slowly decreasing a.In both plots, the shaded area denotes the region where both stable vortex lattice and turbulent state can be encountered.

5 .FIG. 7 .
FIG.7.Stability of regular patterns with arbitrary wavenumber ks.The plots show the regions (blue) where the square vortex lattice is stable as a function of ks and a at different values of λ.The light red area represents the band of unstable modes that might grow from the isoptropic state according to Eq. (5).The dashed line shows the critical wavenumber kc = 1.In compliance with the emergence of turbulence, we observe that the stability balloons rapidly shrink when increasing λ.

FIG. 8 .
FIG.8.Time evolution of the root-mean-square velocity vrms in the simulation investigating hysteresis.In the depicted run, we start from a turbulent state and successively decrease λ in steps of ∆λ = 0.02.After every step, we wait for ∆t = 50000 time units (dashed vertical lines) to ensure a statistically stationary state.Here, we only show the last few steps before the system develops a stationary vortex lattice.The inset shows a magnified view of two subsequent time intervals, demonstrating fluctuations around a mean value of vrms, which increases for decreasing λ.

FIG. 12 .
FIG. 12. Instabilities at the borders of the stability balloon at a = 0.2 and λ = 0. (a)/(b)/(c) show plots for ks = 0.87, i.e., close to the left border of the stability balloon shown in Fig. 7(a).(d)/(e)/(f) show plots for ks = 1.13, i.e., close to the right border.(a)/(d) Maximum growth rate of perturbations to the square lattice with ks = 0.87/ks = 1.13.(b)/(e) Representation of the eigenvector corresponding to the eigenvalue with largest real part at k0 = (0.42, 0)/k0 = (0.18, 18), i.e., for the fastest-growing perturbation.The relative importance of the perturbative modes is determined by the absolute value of the corresponding eigenvector components and indicated by the saturation of the red color.The blue squares show the reciprocal lattice corresponding to the square lattice pattern.The circle indicates |k| = kc = 1.(c)/(f) Vorticity field of the square vortex lattice with ks = 0.87/ks = 1.13 with added perturbations according to the eigenvector shown in (b)/(e).The size of the snapshots is 12π × 12π.