Anisotropic mesoscale turbulence and pattern formation in microswimmer suspensions induced by orienting external fields

This paper studies the influence of orienting external fields on pattern formation, particularly mesoscale turbulence, in microswimmer suspensions. To this end, we apply a hydrodynamic theory that can be derived from a microscopic microswimmer model (Reinken et al 2018 Phys. Rev. E 97, 022613). The theory combines a dynamic equation for the polar order parameter with a modified Stokes equation for the solvent flow. Here, we extend the model by including an external field that exerts an aligning torque on the swimmers (mimicking the situation in chemo-, photo-, magneto- or gravitaxis). Compared to the field-free case, the external field breaks the rotational symmetry of the vortex dynamics and leads instead to strongly asymmetric, traveling stripe patterns, as demonstrated by numerical solution and linear stability analysis. We further analyze the emerging structures using a reduced model which involves only an (effective) microswimmer velocity field. This model is significantly easier to handle analytically, but still preserves the main features of the anisotropic pattern formation. We observe an underlying transition between a square vortex lattice and a traveling stripe pattern. These structures can be well described in the framework of weakly nonlinear analysis, provided the strength of nonlinear advection is sufficiently weak.


Introduction
Active matter exhibits a variety of large-scale self-organized structures that arise due to the interactions between the moving constituents. As shown in experiments, these structures are ranging from dynamical clustering [1,2] and giant number fluctuations [3] to vortices and swirling [4][5][6]. To describe and understand the fascinating collective behavior from a theoretical point of view, models on different levels of detail have been applied, from studies simulating large numbers of individual particles [7,8] to phenomenological approaches [9][10][11][12][13][14] standing on the opposite side of the spectrum. To bridge the gap, many efforts have been made to derive coarse-grained hydrodynamic theories from microscopic models [15][16][17][18][19]. For a selection of recent reviews on the full scope of active matter see [20][21][22][23][24][25][26][27][28].
While many of the self-organized structures in systems of active constituents are well understood, the impact of external fields on the spatiotemporal pattern formation is mostly unexplored. For example, the motion of biological microswimmers such as bacteria or algae cells, is strongly determined by the response to external stimuli. These stimuli are of various origins: for example, swimmers react to concentration gradients (chemotaxis) [29,30], a light source (phototaxis) [31,32] or magnetic and gravitational fields (magnetotaxis [33][34][35][36] and gravitaxis [37,38]). The interplay between externally applied fields and internally generated motion is not only interesting from a fundamental perspective, but also essential for a variety of potential applications. For example, external fields offer the possibility to control the suspension in order to exploit the coherent motion of the swimmers. This includes tasks like cargo-delivery [39][40][41] (e.g. drug transport for medical purposes [42,43]), powering microfluidic devices [44,45] or swimmer-induced mixing on the small scales, where high Reynolds numbers are not accessible [46].
In the present paper, we explore the impact of an external field on a prominent example of active pattern formation, labeled and known as mesoscale turbulence [47]. This state can be observed in a variety of systems, including bacterial suspensions [48][49][50][51] and Janus particles [52]. Mesoscale turbulence is characterized by chaotic vortex structures similar to inertial turbulence occurring in passive fluids at high Reynolds numbers [53], but it has two very unique features: first, it arises in the low Reynolds number regime of bacterial swimming (Stokes flow). Second, in contrast to inertial turbulence, it does not exhibit a spectrum of length scales but displays one characteristic vortex size controlled by the microswimmer details [47,50,54].
From the theoretical side, the main features of mesoscale turbulence have been successfully reproduced by a phenomenological continuum theory for the effective microswimmer velocity [47,[55][56][57][58][59][60][61]. More recently, mesoscale turbulence and vortex lattices were also reported in a model of self-propelled particles with competing alignment interactions [62,63]. Going beyond the purely phenomenological approach, we have recently presented a derivation of a hydrodynamic theory starting from Langevin equations for a generic microswimmer model [18,19]. This theory consists of a dynamic equation for the polar order parameter field coupled to the solvent flow, the latter determined by a modified Stokes equation. In the limit of weak coupling between orientational order and solvent flow, our hydrodynamic theory reduces to the phenomenological model. Here, we extend the theory [19] towards the effect of an orienting external field, the aim being to describe situations where swimmers are subject to an externally applied torque (stemming, e.g. from a magnetic or gravitational field). To this end, we incorporate the field in the Langevin model and derive the additional terms in the dynamic equation for the polar order parameter. This is outlined in section 2 (and appendix A) of the paper.
The remainder of the paper separates into two main parts. In the first part, we investigate the impact of an orienting external field in the full model consisting of the polar order parameter dynamics and an explicit equation for the solvent flow (i.e. a modified Stokes equation). Corresponding numerical results are presented in section 3. We show that, at intermediate field strengths, the emerging patterns become strongly anisotropic as a consequence of the externally broken symmetry. Performing a linear stability analysis (section 4), we find that this is due to the suppression of modes that are perpendicular to the field. For even higher external fields, the mesoscale-turbulent state is completely suppressed and we observe a homogeneous stationary polar state. The analytical results are then used to construct a state diagram of the system. Previously, there has been related work on pattern formation in anisotropic systems with broken rotational symmetry, such as convection in liquid crystals [64] and chemical waves in catalytic surface reactions [65]. Experimental and theoretical studies of such anisotropic systems showed novel phenomena such as ordered arrays of topological defects [66,67], anisotropic phase turbulence [68] and stratified spatiotemporal chaos exhibiting very anisotropic correlations [69,70].
In the second part of the paper (section 5), we switch to a reduced model equivalent to the phenomenological theory, where the only dynamical variable is the effective microswimmer velocity field. The reduced model is significantly easier to handle analytically, but still exhibits the emergence of anisotropic patterns. Using the framework of weakly nonlinear analysis, we investigate the underlying transition from a square vortex lattice to a traveling stripe pattern that occurs upon an increase of the external field. Finally, we present conclusions and an outlook in section 6. The paper is supplemented by seven appendices providing technical details of the calculation.

Hydrodynamic theory
Recently, a phenomenological model [47,[55][56][57][58][59][60][61] has been proposed that reproduces the main features of mesoscale turbulence including the emergence of a characteristic length scale, i.e. vortex size [49]. It is a fourthorder field theory for the divergence-free collective microswimmer velocity and combines the Toner-Tu equation [71] with Swift-Hohenberg-like pattern formation. The main feature of the Toner-Tu equation is the transition from a disordered to a polar state corresponding to collective movement of the swimmers. Higher order gradient terms, as introduced in the Swift-Hohenberg equation, lead to a finite-wavelength instability that is responsible for the collective length scale.
In earlier publications [18,19] we have shown that an equivalent hydrodynamic theory for the microswimmer velocity can be derived via a Fokker-Planck equation approach starting from a generic Langevin model (similar to [15,16]) for a system of microswimmers of length ℓ, diameter d and self-swimming speed v 0 with constant density ρ (see appendix A for a summary of the derivation). The Langevin model includes two types of interactions: first, there are short-range contributions stemming from an activity-driven polar interaction characterized by strength γ 0 and range r [72]. Second, far-field hydrodynamic effects are included via a coupling to the solvent flow field u. The latter is determined via an averaged Stokes equation supplemented by an appropriate ansatz for the active stress tensor containing gradient terms of up to fourth order. The resulting coarse-grained dynamics is then given by an equation for the polar order parameter field P (characterizing the swimmer's mean local orientation) coupled to the solvent flow field u. The effective velocity of the microswimmers is calculated as the sum v 0 P+u. In the limit of weak coupling between the solvent flow and the polar order, the dependence on u can be neglected and the dynamics is adequately described by one field [19]. In contrast to the phenomenological approach, the coefficients of the field equation are directly linked to the parameters of the microscopic Langevin model [19] (see also appendix A).
For the first part of this article, however, we will consider the full model consisting of both the dynamics of the polar order parameter P and the solvent flow u. As shown in [19], the dynamical equation for P(x, t) can be conveniently written in potential form The relaxation term is given as functional derivative of q P P P P P 1 2 Note that, for reasons of brevity, equations (1)-(3) are already written in nondimensionalized form (see below for the definitions of the coefficients and appendix A and [19] for further details). We also note that the potential given in equation (2) should not be regarded as an ansatz for a free energy. Rather, all terms can be derived from the microscopic model as discussed in [19]. The derived equations (1) and (2) exhibit the same features as the phenomenological model: the isotropic-polar transition occurs when α changes sign from positive to negative. The coefficient β of the cubic term determines the saturated value of the polar solution a b -. For sufficiently strong activity Γ 2 becomes negative, which leads to a finite-wavelength instability of the homogeneous state, yielding a typical length scale of 2 2 4 2 p L = -G G . Turbulent dynamics is introduced to the model via the nonlinear advection term on the left-hand side of equation (1), with λ 0 giving the strength of the advection term. Finally, q(x) is a local Lagrange multiplier enforcing the incompressibility condition The assumption of incompressibility is well suited for dense suspensions where density fluctuations become small [47].
In contrast to the phenomenological model, the solvent flow field u enters explicitly through the first term in equation (1). Its coupling to the polar order parameter field P is contained in the generalized derivative where the vorticity tensor and the deformation rate are given by u u , respectively. The modified Stokes equation (see [19]) that determines the solvent flow field u reads where p is an effective pressure and pass s and act s refer to the two parts of the stress tensor, σ, specifically, pass s is an anisotropic passive contribution due to the elongated shape of the particles, and act s is the active stress. The anisotropic passive contribution is treated in detail in the theory of liquid crystals [73]. For the purpose of this study, however, we assume that the passive contribution pass s is not significant compared to the active stress, see [19] for a more detailed discussion. Performing an expansion of the active stress up to fifth order in gradients of the orientation [19] Equations (1)-(6) are rescaled using the microswimmer length ℓ as length scale, the self-swimming speed v 0 as characteristic velocity and ℓ/v 0 as time scale [19].
where P r , c I , r/ℓ, c F and a 0 are dimensionless parameters. The activity is quantified by the persistence number P v r 0 t = ℓ, giving the swimming speed compared to the reorientation time τ. The strength and range of the polar interaction are given by c I and r/ℓ, respectively. For c 1 I < the system favors a disordered, isotropic state, while for c 1 I > it favors an ordered, polar state. The coupling to the solvent flow is characterized by the coefficient c F . For the dependence of c I and c F on microscopic parameters of the microswimmer model, see appendix A. Finally, flow aligning effects are characterized by the shape parameter a 0 which depends on the swimmer aspect ratio (see equation (A.3) in appendix A).
In the present work, we extend equation (1) by terms incorporating an external field that affects the swimmer's orientations. On the microscopic level, we assume that the external field generates a potential for every swimmer which depends on the angle between the field's direction h and swimmer orientation n according to a standard ferromagnetic coupling, i.e. n h ext F µ -· . This term also occurs in some passive liquids in an orienting field (e.g. ferrofluids [74,75], which are suspensions of ferromagnetic colloids in a passive solvent) and was recently considered in the context of active fluids. Indeed, the same ansatz has been used to describe chemotactic [30] and magnetotactic [76] bacteria. It is in principle applicable to any microswimmer suspension subjected to an aligning torque exerted by an external field, as it occurs, e.g. in magnetotaxis [33][34][35][36], phototaxis [31,32] or gravitaxis [37,38]. Introducing the external potential in the Langevin equations and performing the same steps as in the field-free case, one arrives at the order parameter equation (for details see appendix A) The additional term on the right-hand side of the evolution equation is given as where B 0 denotes the dimensionless external field strength. The first term on the right-hand side of equation (9) increases the polar order in the direction of the external field, similar to what happens in passive fluids. The second term arises due to the conservation of the unit vector n, i.e. the microswimmer's orientation. It leads to a saturation of P with increasing B 0 , as we will later see (compare figure 2). Finally, the third term in equation (9) arises as a consequence of the closure scheme applied for the nematic order parameter tensor Q (see appendix A). Clearly, this term incorporates a coupling to the solvent flow field. Physically, it adds a reduction of polar order due to the flows generated by the swimmers. In principle, changes on the flow field due to the external field will couple back to the stress tensor. However, this is a higher order (feedback) effect, which we will neglect for the purpose of this study. Also note that the term describing the impact of the external field g cannot be straightforwardly integrated into F due to the closure relations used in the derivation, see appendix A and [19]. This indicates again that F cannot be simply interpreted as a free energy.

Numerical observations
The aim of the present study is to explore the impact of a spatially homogeneous stationary field on the mesoscale-turbulent state observed in the absence of a field. As a starting point, we will discuss the dynamical behavior that can be observed based on numerical solution of equation (8) for the polar order parameter field P, coupled to the Stokes equation (6) for the flow field u. The solution is performed in two-dimensional space (for the numerical methods, see appendix B), where the two dimensions are defined by the directions parallel and perpendicular to the field, i.e. x  and x^. The most instructive quantity to visualize the dynamics is the vorticity field, which is defined as a three-dimensional axial vector. In the following, we will refer to the z-component of the vorticity P z ( ) (where z is both perpendicular to x  and x^) as scalar vorticity or just vorticity. We start with the field-free case, B 0 = 0. Here, we focus on parameters where the system is in a mesoscale turbulent state. In figure 1(a) a snapshot of the (scalar) vorticity of the polar order parameter field is shown. For a more descriptive visualization of the dynamics in the absence of a field see the supplementary movie 1 available online at stacks.iop.org/NJP/21/013037/mmedia. For better visibility, a section of the field is presented in a larger version in figure 1(d). Here, we also visualize the (vectorial) order parameter field by arrows, with the arrow length indicating the magnitude of the polar order. The observed dynamical state is characterized by the formation, motion and decay of clockwise (blue) and counter-clockwise (red) rotating vortices. The emerging patterns are dominated by one characteristic length or vortex size depending on the microscopic parameters of the swimmers [18,19,47,50,54], hence the name mesoscale turbulence. This stands in contrast to inertial turbulence observed in the Navier-Stokes equation where one finds a broad spectrum of vortex sizes [53]. For a more detailed discussion on the turbulent state without the influence of external fields, see [47,58,59,61] for the phenomenological model and [18,19] for the present model. Note that, compared to most of the listed publications, the coefficient λ 0 is rather large and, therefore, the shape of the vortices is highly irregular.
Turning on the external field, B 0 >0, introduces several new features. First, the external field induces a net polar order in the system, that is propulsion speed leads to a local transport in the direction of the polar order given by v P 0 . Thus, as a consequence of the generated net polar order, we observe a net transport in the direction of the field. Second, the overall rotational symmetry is broken, which leads to the emergence of asymmetric patterns. In figures 1(b) and (c) this is illustrated by snapshots of the vorticity at finite (nonzero) field strengths. In contrast to the case B 0 =0, we here observe the formation of elongated structures in the vorticity field, or, more precisely, a highly irregular stripe pattern with numerous defects. In figure 1(c), where the field strength is larger compared to (b), the number of defects is smaller and the stripe pattern more regular. Interestingly, the defects are elongated in the direction parallel to the field. The enlarged sections in figures 1(e) and (f) additionally show the polar order in the field's direction. For a visualization of the transport of the patterns, see the supplementary movies 2 and 3. They show the dynamics for the same values of the field strength as represented in figures 1(b) and (c), i.e. B 0 =0.6 and B 0 =0.8, respectively. In the remainder of this work we will elucidate the observed dynamical features using analytical methods, particularly linear stability analysis and weakly nonlinear analysis.

Homogeneous stationary solution
If the activity is sufficiently weak or the external field sufficiently strong, we numerically observe a homogeneous stationary state. This state is the starting point of our linear stability analysis.
Clearly, the external field breaks the rotational symmetry of the system. Thus, it is useful to distinguish between components of the polarization parallel and perpendicular to the field, i.e. P P P , = ( ). To calculate the homogeneous stationary solution, we assume a quiescent state, where the solvent velocity field vanishes, i.e. u=0. Further, the pressure p=p 0 and Lagrange-multiplier q=q 0 are set constant in space and time. Then, equation (8) reduces to The real positive solution of equation (

Linear stability analysis
In order to determine the linear stability of the homogeneous stationary solution P P , 0 0 = ( ), q=q 0 , u 0, 0 = ( ), p=p 0 , we consider small perturbations, i.e. For all perturbations, we make the ansatz P P q u u p P P q u u p , , , ). We now insert equations (12) and (13) into equations (2), (4), (6), (8) and (9) and linearize with respect to the perturbations. As shown in appendix C, the perturbations of the velocity field δu, the pressure δp and the Lagrange multiplier δq can be related to the perturbations of the polar order parameter δP. As a result (see appendix C) one obtains a linearized system involving only δP, that is The imaginary part k Im s ( ), which determines the linear traveling speed c 0  , is given by It is seen that k Re s ( ) depends not only on the magnitude k k = | | of the wavevector but also on its direction. For B 0 =0 and c F =0, one reproduces the growth rate obtained in [55]. Note that, when the system transitions into a polar state in the absence of the external field, i.e. c F >1 (α<0), the direction of P is not uniquely defined (i.e. there is continuous degeneracy), but instead chosen spontaneously by the system itself. In this case, we have to set x  parallel to the direction of spontaneous order and x^perpendicular to that direction.

State diagram
Due to the explicit dependence of the growth rate k Re s ( ) on the wavevector's direction, it makes sense to distinguish between two limiting cases illustrated in figure 3: perturbations with a wavevector that is purely parallel to the external field, i.e. k k =  , k 0 = , or purely perpendicular to the field, i.e. k 0 = The incompressibility condition of the order parameter field, i.e. k P , then dictates the form of the respective perturbations. For a parallel wavevector, the component P d  must vanish. Therefore, the linearly unstable pattern is a perturbation in the perpendicular direction with periodicity in the parallel direction (see figure 3(a)). The corresponding wavelength is given by . This type of pattern manifests itself as stripes in the vorticity of the field. Analogously, for a perpendicular wavevector, the component P d^must vanish and the pattern is a perturbation along the field with periodicity in the perpendicular direction ( ). The occurrence of both, modes in the parallel and in the perpendicular direction, results in a square vortex lattice.
The results obtained from linear stability analysis and the distinction between the two limiting cases for the perturbation enable the construction of a state diagram in the plane spanned by P r and B 0 , see figure 4(a). The diagram is supplemented by plots of the growth rate as function of the wavevector, see figures 4(b)-(e). Regions where the homogeneous stationary state described in section 4.1 is stable are indicated by a white background color. Considering first the case B 0 =0, we observe an instability at a critical value (P 5 r » ) where mesoscale turbulence sets in. This is reflected by a band of unstable modes (see figure 4(c)), one of which grows the fastest and yields a typical vortex size 2 2 4 2 p L = -G G . Note that, without an external field, the full rotational symmetry is still intact and the two curves for parallel and perpendicular wavevector perturbations coincide. For a (numerical) snapshot of the order parameter field P and its vorticity P ´in the mesoscale turbulent state at Keeping the persistence number at a constant value within the zero-field turbulent state, e.g. P 8 r = , and increasing the external field from zero results in a shift of the growth rate curve σ R depending on the wavevector's direction (see figure 4(d)). It is seen that perturbations with perpendicular wavevector become suppressed above a field strength of B P Moreover, we observe a very intriguing feature in the case of a perturbation with k h  . Here, the growth rate curve is first shifted upwards for intermediate field strengths before being shifted down for stronger fields. This is due to the explicit coupling of the solvent flow to the generated polar order in the system. As discussed in a variety of publications [9,21,23], active suspensions with uniform orientational order such as active nematics [77] are intrinsically unstable. The interplay between actively generated flow and the resulting flow alignment of the elongated particles destabilizes the uniaxial order. In active nematics, this results in a bend instability for particles that generate an extensile active stress and a splay instability for particles that generate a contractile active stress [21]. In our model of polar microswimmers, a similar mechanism leads to the upwards shift of the growth rate curve of parallel modes for intermediate field strengths. The external field generates orientational order, which in turn induces a solvent flow that destabilizes the homogeneous state due to the flow alignment of . Note that the diagram is independent of the magnitude of the nonlinear advection term 0 l . the swimmers. For higher field strengths, however, this interplay between polar order and solvent flow is outweighed by the other terms in the growth rate (equation (15)) that suppress the instability. This behavior is reflected in the state diagram in figure 4(b) by the pronounced 'bulge' for persistence numbers below P 5 r » . Finally, note that the growth rate Re s (see equation (15)) does not depend on the nonlinear advection term in equation (1). The imaginary part Im s , however, is strongly dependent on to the parameter λ 0 characterizing the magnitude of the advection term (see equation (16)). Thus, in the framework of linear stability analysis, the influence of the nonlinear advection term is restricted to transporting the emerging patterns with a traveling speed of c 0  , given by equation (16), in the direction of the external field. Comparing the state diagram in figure 4 to the full numerical solution of equations (8) and (6), we find that the analytical calculations are consistent with the numerical results. They accurately predict the onset of the mesoscale-turbulent state, the complete suppression of the instability for high field strength and the region of anisotropic patterns for intermediate field strength (compare figure 1).

Reduced model
In order to characterize the emergence of patterns in more detail, it seems appropriate to perform a weakly nonlinear analysis. However, it turns out that this is an extremely difficult task when starting from the model equations discussed so far. For the remainder of this work we thus consider a reduced model, which still gives the essential physics, i.e. the transition from rotationally symmetric to asymmetric patterns and the complete suppression of the turbulent state for increasing field strength.
As we discussed in our previous publication [19], the response of the flow field to the forces exerted by every swimmer on the surrounding fluid scales with the coefficient c F . The latter is proportional to the ratio of active to viscous forces. In the limit c 1 F  , the collective dynamics of the suspension depends solely on the dynamics of the polar order parameter. Introducing an effective microswimmer velocity proportional to P then reproduces the phenomenological model briefly introduced in section 2, with the notable difference that we can calculate all the coefficients as functions of parameters of the microswimmer model. This phenomenological model, consisting of the dynamics of only one vector field, is indeed frequently used to describe mesoscale turbulence [47,[55][56][57][58][59]61]. Note that the coefficient Γ 2 has to be negative in order to obtain a mesoscale-turbulent state. Thus, additionally to the limit c 1 F  , the product P c r F still has to be sufficiently large compared to c I (see equation (7)).
The advantage of this scaling is that the number of independent coefficients is now reduced to four: First, the coefficient λ determines the strength of the nonlinear advection term. Second, the coefficient a characterizes the distance to the onset of the instability at B 0 0 = . For a<0, the system favors the isotropic homogeneous state, while for a>0 the system develops mesoscale turbulence. Note that for a>1, the homogeneous stationary solution becomes polar. In this work, we will focus on the case a<1. Third, the coefficient b of the cubic term in v is responsible for the saturation of the emerging patterns. Fourth, the external field strength is given by B 0 . The scaling of the Lagrange multiplier q enforcing the incompressibility v 0  = · is irrelevant, thus, we keep the same notation as before. For the explicit dependencies of the four remaining coefficients λ, a, b and B 0 on the coefficients of the full model for the dynamics of P given in equations (8) and (9), see appendix D.
The reduced and rescaled model is significantly easier to handle than the full version but still exhibits the emergence of anisotropic patterns due to the external field. Note that, for the sake of brevity, we use the same notation for the growth rate and related concepts as for the full model. The state diagram obtained from the stability analysis of equation (17) is shown in figure 5. As in the full model we find a region where perturbations with a parallel wavevector grow but perpendicular wavevector perturbations decay. However, compared to figure 4, the bulge on the left-hand side is absent. This is because the reduced model is solely given by equation (17) and, thus, the coupling of the externally generated polar order to the Stokes equation (6) is neglected. As discussed in section 4.3, this coupling is essential for the mechanism producing the bulge.

Weakly nonlinear analysis
To analyze in more detail the emerging patterns in the reduced model, we now perform a weakly nonlinear analysis [78,79]  and a>0), where the system's rotational symmetry is broken and modes with a wavevector parallel to the field dominate the dynamical behavior. Here, we observe the emergence of stripes stacked in the field's direction. The effective microswimmer velocity field for these two regular patterns is visualized in figures 6(a) and (b), respectively. Various technical details related to the weakly nonlinear analysis are provided in appendix E. Although equation (17) is similar to the Swift-Hohenberg equation, the dominant pattern in the field-free case is a square vortex lattice and not (as in the Swift-Hohenberg equation) hexagons or stripes. This is due to the fact that the order parameter considered in this work is a divergence-free vector field and not a scalar. For example, neighboring vortices always favor opposing directions of rotation, i.e. clockwise and counter-clockwise, and this antiferromagnetic configuration can never be realized in a hexagonal pattern. In principle, a potential for the vorticity can be constructed that favors a certain direction of rotation, which then allows for a hexagonal pattern consisting only of clockwise or counter-clockwise rotating vortices. However, for the parameter regime considered in this work, our theory does not break this symmetry and therefore we only exhibit a square lattice (see [61] for further discussion).  figure 4) the bulge for P 5 r < which corresponds to a<0 is absent. The coefficient of the cubic term is set to b=0.1. Note that the diagram is independent of the magnitude of the nonlinear advection term λ.

Case of zero external field
In the absence of an external field, the onset of the instability occurs at a=0 (see figure 5). Due to the rotational symmetry, the weakly nonlinear analysis for this case is essentially standard (see [78,79]). We introduce a small parameter ε characterizing the distance to the bifurcation via the growth rate of critical modes B k a 0, 1 2 Re 0 (18) and (19)). This allows the definition of a slow time scale T t 2 e = and corresponding spatial variable X x e = (motivated by the fact that the growth rate scales in lowest order quadratic in k). These long time and space variables characterize the scales on which the amplitude of the emerging patterns evolves. The next step is to expand the effective microswimmer velocity field v in orders of ε. In the present case, unstable modes have no preferred direction and the emerging pattern is a square vortex lattice (see figure 6(a)) which is formed by critical wavenumber modes perpendicular to each other. Although their directions are arbitrary at B 0 =0, for convenience, we choose one mode parallel and one mode perpendicular to the direction of the external field which will be set to finite values in the next section. We thus also consider the (scalar) components v  and v^of the effective velocity field. Taking into account that the incompressibility condition, i.e. v 0  = · , has to be satisfied, the expansion reduces in lowest order to where A  and A^are the amplitudes of the two modes and c.c. denotes the complex conjugate. In the exponential functions, we have used k c =1. Now, we insert the expansion equations (21) into (17) and use the definitions for the slow time scale T and corresponding spatial variable X. Matching terms with the same order of ε yields the amplitude equations for the parallel and perpendicular mode in where we already scaled back to the fast time and length scales, t and x, respectively. It is seen that the two amplitude equations (22) correspond to two coupled Ginzburg-Landau equations [80]. This result was already obtained in [61], where the pattern formation in the reduced model equation (   where we have used equation (19). The slow time scale is again defined by T t 2 e = . The scaling of the corresponding spatial variable is given by t X x v g e = -( ), where we introduced the group velocity v g for the following reason: in contrast to the case B 0 0 = , the system at B 0 0 > displays a net polarization, and the emerging stripe pattern (see figure 6(b)) is traveling in the field's direction with the traveling speed c  . We also expect modulations of the pattern to travel through the system. The corresponding velocity of the modulations is denoted as group velocity v g and is, at this point, still to be determined. The next step is to expand the effective microswimmer velocity field v with respect to the distance to the bifurcation ε. Here, we use an ansatz which incorporates only parallel modes characterized by multiples of the critical wavenumber (k c =1). In contrast to equation (21), perpendicular modes are not considered, since they decay in the considered region in the state diagram. The full expansion is given as where we already scaled back to the fast scales. We find that the coefficient g>0 (see equation Thus, the amplitude of the emerging patterns modifies the speed with which they travel through the system. A second unique feature of the symmetry-broken case is the anisotropic diffusion of pattern modulations, as reflected by the difference of the diffusions constants D  and D^, respectively, At the parameters considered, the coefficient for parallel transport is approximately one order of magnitude higher than for perpendicular transport. This leads to an interesting additional effect, which we discuss in detail in section 6.

Transition between square vortex lattice and stripe pattern
Having obtained the amplitude equations (22) and (25) for the square vortex lattice and the traveling stripe pattern, respectively, there remains the question about their predictive power. Indeed, one would expect that only situations with weak nonlinearities, i.e. situations near the onset of the respective instabilities, are adequately described in the framework of weakly nonlinear analysis [50]. These weak nonlinearities include the quadratic and the cubic term in equation (17), which are responsible for the saturation of the amplitudes, as well as the linear part of the self-advection term that is responsible for the transport of the patterns in the direction of the field. This linear part is proportional to the net velocity v á ñ and can be written as v v l á ñ  · (compare equation (17)). However, the weakly nonlinear analysis does not capture the full nonlinearity of the advection term, as it is hard to handle analytically. From classical turbulence theory it is known that this term transfers energy that is inserted into the system between different length scales [53]. This contradicts the basis of weakly nonlinear analysis, i.e. the assumption that the dominant mode is given by the critical wavenumber. For example, for two-dimensional flows, the energy cascade decreases the dominant wavenumber in the system, (see also figure F1 in appendix F and [61] for a more detailed discussion).
However, if the strength of the advection term determined by the parameter λ in equation (17) is small, the full nonlinear nature of this term is expected to be less relevant. In this case, we expect the formation of a regular square vortex lattice for B 0 0 = and traveling stripes for B B B  (17) numerically for small values of λ, we observe a regular square vortex lattice also in the presence of a finite external orienting field, i.e. B 0 0 > , provided that the field's magnitude is sufficiently small. Note that the formation of an asymmetric lattice, that is a lattice where the two directions exhibit different characteristic lengths, is not observed in simulations. This is due to the fact that the external field does not change the critical mode but only shifts the entire growth rate curve (see also section 4.3). Starting from the vortex lattice and increasing the field strength, a transition to regular stripes traveling in the field's direction occurs at a critical strength B 0 * . This transition can be observed in figure 7 where the maximum vorticity (in the rescaled model simply given by the sum of amplitudes of parallel and perpendicular modes, ) is plotted over the external field strength. Numerical results, obtained by solving the full dynamical equation for v (equation (17) Interestingly, the transition field strength B 0 * is not given by B 0^, that is, the field strength where, according to linear stability analysis of the homogeneous stationary solution, perpendicular modes start to grow (see section 4.2). The reason for this discrepancy is that perpendicular modes become suppressed by parallel modes already at field strengths smaller than B 0^. The transition field strength B 0 * is obtained by taking the coupling between the modes into account and performing a linear stability analysis of the stripe pattern with respect to perpendicular modes (see appendix G). In order to check for dependencies on the initial values in the numerically obtained data points in figure 7, we performed the numerical solution twice, starting from a regular vortex lattice and a stripe pattern. We found no impact of the chosen initial values indicating the absence of hysteretic behavior. As visible in figure 7, once the vortex lattice is fully developed, the saturated value is given by equations (22), which means the external field has no influence on the amplitude. Instead, preliminary numerical results show that the external field only distorts the lattice. This intriguing effect will be discussed elsewhere.

Conclusions
This article studies the impact of a homogeneous, stationary external field on pattern formation in suspensions of microswimmers that exhibit mesoscale turbulence in the field-free case. Based on a numerical solution of the full model presented in section 2 and a linear stability analysis, we observe that, in the presence of an orienting field, the patterns become anisotropic. From a mathematical point of view, the growth rate of modes becomes dependent on the wavevector's direction due to the broken rotational symmetry. In particular, modes perpendicular to the applied external field are suppressed for intermediate field strengths leading to the dominance of modes parallel to the field. The resulting structures can be described as stripes that travel along the field's direction. For even higher field strengths, the instability is completely suppressed and we observe a homogeneous stationary polar state.
In the second part of the paper, we have presented a weakly nonlinear analysis of a reduced model for the effective microswimmer velocity. The model is significantly easier to handle analytically, yet still exhibits anisotropic pattern formation, which we have analyzed by deriving the amplitude equations (22) and (25). Upon an increase of the external field, there is a transition form a square vortex lattice to a traveling stripe pattern. The regularity of these patterns is strongly dependent on the coefficient λ, that determines the strength of the nonlinear advection term. When λ is increased, defects are generated and the patterns become less regular. In this case, the amplitude equations, (22) and (25), lose their validity. This boils down to the problem already discussed in section 5.2: the nonlinear energy transfer between scales leads to a shift of the dominating mode contradicting the ansatz leading to the amplitude equations (see also appendix F).
However, we can still observe one unique feature of the amplitude equation (25) for the stripe pattern, even for large values of λ: the equation is strongly anisotropic as is reflected by the difference of the diffusion coefficients, D  and D^. Indeed, the transport of modulations of the pattern in the direction parallel to the external field is approximately one order of magnitude faster than perpendicular to the field (see also equation (27)). As a consequence, defects in the patterns, generated by the nonlinear advection term, will be elongated by a factor of D D in the field's direction. This is consistent with numerical observations shown in figure 1(c) for the full model.
As explained above, the amplitude equation (25) does not incorporate the full nonlinear advection term. This is apparent from the fact that it corresponds to the real-valued Ginzburg-Landau equation, which does not generate defects. The complex two-dimensional Ginzburg-Landau equation, however, exhibits a variety of chaotic states including one denoted as defect turbulence [80]. A preliminary comparison of this state with the dynamics of the amplitude modulations of the stripe patterns observed in the present system shows quite similar spatial structures. The derivation of an amplitude equation valid even for higher values of λ presents an interesting future challenge.
As discussed, the external field induces polar order in the system, which leads to net transport in the field's direction. This net transport can be modified by emerging patterns, as was demonstrated in a recent model for magnetic microswimmers, where band-like structures reflected in the inhomogeneous density of swimmers decrease the net polarization [76]. We also find an influence of emerging patterns on transport properties in our model that assumes a constant density of microswimmers on the coarse-grained level: the emerging stripe pattern in the polar order parameter field influences the traveling speed (see equation (26)). Moreover, preliminary numerical results show that the mean transport in the system is impacted in a quite complex manner. Further exploring this feedback will be especially relevant for applications, where transport properties are essential.
Anisotropic patterns emerging in microswimmer suspensions subjected to external fields have indeed already been observed in magnetotactic bacteria [81,82]. However, in these experiments, band-like structures in the swimmer density were found, whereas in our case, we are dealing with a purely orientational effect. The stripes, visible in the vorticity, correspond to wave-like patterns in the polar order parameter field.
The overdamped motion of a swimmer σ in an ensemble of σ=1, ..., S identical swimmers is given by the Langevin equations for the position X s and the orientation N s , respectively, where the functions x s and h s denote Gaussian white noise generating diffusion in the translational and rotational motion. Note, that the overdamped Langevin equations are already divided by the translational and rotational friction coefficients, respectively. Self-propulsion is introduced via the first term on the right-hand side of equation (A.1), involving the self-swimming speed v 0 . Additionally, the swimmer is transported by advection with the averaged surrounding flow field t u x, ( ). The orientational motion given by equation (A.2) is determined, first, by rotation and alignment with the flow field. This is reflected by the terms involving the vorticity u u ( )] and deformation rate u u , respectively. In analogy to Jeffrey's theory for oscillatory tumbling motion of elongated particles in flow [83][84][85], the shape parameter a 0 is given as a function of the aspect ratio, defined by length ℓand diameter d of the swimmers, Further, the projector I N N P = -s s (where I is the unit matrix) is introduced to conserve the length of the unit vector N s . The second deterministic contribution to the orientational motion is the conservative potential Φ involving all swimmers. In the present study we consider both, pair interactions and an external contribution, that is, where 0 g is the strength of the interaction, and x 1 Q = ( ) for x 0  and zero otherwise. Note that this is a simplified ansatz that describes the polar alignment of neighboring swimmers due to near-field hydrodynamics interactions [72], an effect that has been observed experimentally [3]. Finally, the external potential for swimmer μ is given by Coming back to equations (A.1) and (A.2), far-field hydrodynamic interactions are taken into account by the coupling terms involving the average surrounding flow field u. At low Reynolds numbers, this field is determined by the Stokes equation, augmented by an active contribution to the stress tensor, which, in turn depends on the order parameter. This coupling eventually leads to equation (6). Due to the lengthy derivation of the active stress we refer to our previous publication [19] for details. The next step is to obtain the Fokker-Planck equation for the one-particle probability density function t x n , , ( ), which gives the probability to find a swimmer at position x with orientation n at time t. To this end, we assume a constant swimmer density t x, r r = ( ) and employ a mean-field approximation to treat the twoparticle correlations stemming from conservative interactions. We then project on to orientational moments n, nn, ... of t x n , , ( ), which are directly connected to the polar order parameter P and the nematic order parameter Q via P n Q nn I , 3 .
Following our previous studies [18,19] we apply a closure relation for the nematic order parameter, where the coefficients q and λ K can be calculated analytically [19]. This is an extension of the Doi closure for passive particles, which incorporates the fact that active particles always generate flow gradients affecting the local nematic order. For moments higher than the second we apply the so-called Hand-closure [86]. The details of the procedure are discussed in [19]. Finally, we arrive at a closed dynamic equation for the polar order parameter P, given by equation (8)  where f 0 denotes the strength of the force dipole exerted by every swimmer on the surrounding fluid, and μ is the effective viscosity of the suspension [19].

Appendix B. Numerical methods
Our numerical results are obtained by solving the dynamical equations (8) or (17) using the Runge-Kutta-Fehlberg method (RKF45) [87]. We use a finite-difference discretization of the spatial derivatives on a periodic grid consisting of 256×256 points. For the purpose of visualization, we have increased the resolution in figure 1 by cubic interpolation. The incompressibility condition is enforced by a pressure-correction method [88]. In the case of the full model, the Stokes equation (6) is solved applying a stream-function approach [88]. If not otherwise stated, we apply the homogeneous stationary solution introduced in section 4.1 as initial conditions and add small random variations.

Appendix C. Details of the linear stability analysis
In the following, we present the calculation of the complex growth rate σ (see equations (15) and (16)) determining the stability of the homogeneous stationary solution. To this end, we will first eliminate the dependence of the linearized system on the perturbations u dˆ, p dˆand q dˆ, thus reducing the set of variables to just the perturbation of the polar order parameter, P dˆ.
As a first step we insert the perturbed solution (see equation (12) and (13) in the main text) into the Stokes equation (6) and linearize, yielding We multiply equation (C.1) by the wavevector k and employ the incompressibility conditions which, for the perturbed system, yield k u 0 d = ·ˆand k P 0 d = ·ˆ. From this, we find that the pressure perturbation amplitude must satisfy p 0 d = . Equation (C.1) then yields the perturbation u dˆas a function of the perturbation of the polar order parameter, Equation (C.4) still involves the perturbation q dˆof the Lagrange multiplier. Utilizing again the incompressibility condition for the polar order parameter field P, yielding k P 0 d = ·ˆ, and equation (C.4), we find that the perturbation q dˆmust satisfy Inserting equation (C.9) into (C.4) we can identify the projector k I kk k 2 P = -( ) | | . With this, we obtain equation (14) in the main text, giving a linearized system involving only perturbations of P. The complex growth rate can now be readily calculated as solution of the eigenvalue problem (equation (14)) via ), the derivatives in the dynamic equation (17)  . This is due to the rotational symmetry and, thus, the absence of net transport in the system. In the present case, however, the rotational symmetry is broken and, near the onset of the instability at B 0  , only perturbations with a parallel wavevector grow. Therefore, we expand the effective velocity field v v v , = ( )in orders of ε incorporating only parallel modes characterized by multiples of the critical wavenumber. The full expansion is given in equation (24) in the main text. Similarly, we expand the Lagrange multiplier q that enforces the incompressibility constraint using only parallel modes, , e c.c.  (20)). Now, we insert all expansions (equations (24), (E.2) and (E.3)) into equation (17) and replace all derivatives via equation (E.1). Sorting the resulting terms according to powers of ε (n) and modes (m) we obtain solvability conditions. In the following, we successively go through the resulting equations: • In zeroth order, 0  e ( ), we recover equation (18) which determines the homogeneous stationary solution V 0 .
• In first order, where we used the definition of ε according to equation (23). Thus, we find that v 1,1 actually contributes in 3  e ( )to the amplitude equation obtained below.
• In third order, Appendix F. Shift of the dominating mode As discussed in sections 5.2, the nonlinear advection term leads to energy transfer between different scales once its strength, λ, is large enough. As a consequence, the dominating mode in the system shifts. This shift is towards larger scales, i.e. smaller wavenumbers, for a two-dimensional system [53,61]. To illustrate this point, we plot the dominating mode as function of the strength of the advection term λ in figure F1. The data points are obtained numerically on the basis of the reduced and rescaled model equation (17).

Appendix G. Stability of stripe pattern
The emerging pattern near the bifurcation at B 0  manifests itself as stripes in the vorticity. With decreasing external field strength, the stripe pattern becomes unstable against the formation of a square vortex lattice. The corresponding critical field strength can be obtained by performing a suitable linear stability analysis. In contrast to section 4.2, we here add a perturbation to the stripe pattern characterized by the stationary solution

A g
Re s =  of the amplitude equation (25) (and not to the homogeneous stationary solution V 0 ). The full stripe pattern solution on the level of the collective microswimmer velocity v is given as Figure F1. Dominating mode k dom in the reduced system (equation (17)) as function of the strength λ of the advection term for the field-free case, B 0 0 = . The remaining parameters are a=0.8 and b=0.1. The data points are obtained by numerically solving (equation (17)) and determining the first maximum max L of the spatial correlation function v x v x x á + D ñ ( ) · ( ) (where ... á ñdenotes the spatial and temporal average). The dominant mode then follows as k 2 dom max p = L . To reduce the number of transient defects remaining from the initial conditions, we start from a regular vortex lattice with small random variations. For small λ, the characteristic length scale of the patterns is given by the critical mode k 1 c = . For larger λ, turbulent motion sets in and energy is transfered to larger scales, i.e. smaller wavenumbers. , we find v 0 d = . We insert equation (G.2) into the dynamic equation for v (equation (17)) and obtain, after linearization, the growth rate of perpendicular modes with critical wavenumber k c =1, The critical field strength B 0 * defining the transition from the stripe pattern to the vortex lattice is then obtained by setting 0 s = . Note that s^depends on the amplitude of the parallel modes, A. This coupling is the reason why the field strength at the transition, B 0 * , is smaller than B 0^.