Asymmetric pattern formation in microswimmer suspensions induced by orienting 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 [Phys. Rev. E 97, 022613 (2018)]. 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 asymmetric 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 spatio-temporal 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 asymmetric 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.
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 asymmetric 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 fourth-order field theory for the divergence-free collective microswimmer velocity and combines the Toner-Tu equation [64] 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 [65]. 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. The dynamical equation for P(x, t) can be conveniently written in potential form The relaxation term is given as functional derivative of 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 −α/β. 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 . 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 ∇ · P = 0. 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 Ω = 1 2 [(∇u) T − (∇u)] and Σ = 1 2 [(∇u) T + (∇u)], respectively. The modified Stokes equation (see [19]) that determines the solvent flow field u reads where p is an effective pressure, and the gradient terms of P stem from an expansion of the active stress. Equations (1) to (4) 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]. The coefficients can then be written as where P r , c I , r/ℓ, c F and a 0 are dimensionless parameters. The activity is quantified by the persistence number P r = v 0 τ /ℓ, 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 I < 1 the system favors a disordered, isotropic state, while for c I > 1 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 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., Φ ext ∝ −n · h. This term also occurs in some passive liquids in an orienting field (e.g., ferrofluids [66,67], 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 [68] 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 righthand side of equation (7) 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 (7) 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.

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 (6) for the polar order parameter field P, coupled to the Stokes equation (4) for the flow field u. The solution is performed in two-dimensional space (for the numerical methods, see Appendix B).
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. 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 P = A −1 dxP(x, t) = 0 (where A is the area). The swimmer's self-propulsion speed leads to a local transport in the direction of the polar order given by v 0 P. 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 figure 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 figure 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 figure 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 (6) reduces to The real positive solution of equation (8) defines the homogeneous stationary state P = (P 0 , 0). This solution P 0 is plotted versus the external field strength B 0 in figure 2.
As expected, the polar order grows upon increasing the field strength. Dividing by B 0 and evaluating the limit B 0 → ∞, equation (8) simplifies to The solution of equation (9) defines the saturation value P sat

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, δû , δû ⊥ , δp) e σt+ik·x , with wavevector k = (k , k ⊥ ). We now insert equations (10) and (11) into equations (2), (3), (4), (6) and (7) 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 where the projector Π(k) = I − kk/|k| 2 arises as a consequence of the incompressibility of the order parameter field, and the 2×2-matrix M(k) is the Jacobian. The components of M(k) are given in equation (C.5) in Appendix C. We obtain the complex growth rate σ = σ Re + iσ Im as a function of the wavevector k by calculating the eigenvalues of the matrix Π(k) · M(k). The real part σ Re (k), which determines the actual growth of a mode k, is given by The imaginary part σ Im (k), which determines the linear traveling speed c 0 , is given by It is seen that σ Re (k) 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].

State diagram
Due to the explicit dependence of the growth rate σ Re (k) 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, k ⊥ = k. The incompressibility condition of the order parameter field, i.e., k δP + k ⊥ δP ⊥ = 0, then dictates the form of the respective perturbations. For a parallel wavevector, the component δP 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 Λ = 2π/k. This type of pattern manifests itself as stripes in the vorticity of the field. Analogously, for a perpendicular wavevector, the component δP ⊥ must vanish and the pattern is a perturbation along the field with periodicity in the perpendicular direction (Λ ⊥ = 2π/k) [see figure 3(b)]. 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) to (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 r ≈ 5) 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 . 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 B 0 = 0, see figure 1(a). Keeping the persistence number at a constant value within the zero-field turbulent state, e.g., P r = 8, 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 ⊥ 0 (P r ), see dashed line in figure 4(a). In contrast, perturbations with parallel wavevector still grow until a larger field strength of B 0 (P r ) is reached, see solid line in figure 4(a). Between the two "critical" field strengths B 0 and B ⊥ 0 , the growth of parallel modes leads to asymmetric patterns which are elongated in the perpendicular direction. This is consistent with our numerical observations shown in figure 1(b) and (c). Interestingly, the characteristic length given by the critical mode k c , i.e., the maximum of σ Re (k), is not influenced by the external field [compare figure (4) 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 [69] 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 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 13] 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 r ≈ 5.
Finally, note that the growth rate σ Re [see equation (13)] does not depend on the nonlinear advection term in equation (1). The imaginary part σ Im , however, is strongly dependent on to the parameter λ 0 characterizing the magnitude of the advection term [see equation (14)]. 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 (14), in the direction of the external field. Comparing the state diagram in figure 4 to the full numerical solution of equations (6) and (4), 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 asymmetric 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.
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 F ≪ 1, 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 [47,[55][56][57][58][59]61] briefly introduced in section 2, with the notable difference that we can calculate all the coefficients as functions of parameters of the microswimmer model. Note that the coefficient Γ 2 has to be negative in order to obtain a mesoscale-turbulent state. Thus, additionally to the limit c F ≪ 1, the product P r c F still has to be sufficiently large compared to c I [see equation 5].
In this work, we rescale space and time by the length and time scale of the pattern formation, i.e., the inverse of the critical mode k c = −Γ 2 /(2Γ 4 ), and the inverse of the corresponding maximum of the growth rate σ Re (α = 0, B 0 = 0) = Γ 2 2 /(4Γ 4 ). Further, we rescale by the saturation value for the polar order parameter at B 0 → ∞ and obtain the scaled collective microswimmer velocity via v = P/P sat 0 , where P sat 0 = 5/(3c I ). The rescaled equation for v is then given by 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 atB 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 byB 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 andB 0 on the coefficients of the full model for the dynamics of P given in equations (6) and (7), see Appendix D.
The reduced and rescaled model is significantly easier to handle than the full version but still exhibits the emergence of asymmetric patterns due to the external field. The homogeneous stationary solution Analogous to the full model, performing a linear stability analysis yields the complex growth rate with real part Similar to the corresponding function σ k of the full model [compare equation (13)], equation (17) yields a finite-wavelength instability. Due to the present scaling, the critical mode is now given by k c = 1. Further, we obtain the linear traveling speed as 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 (15) 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 (15) and, thus, the coupling of the externally generated polar order to the Stokes equation (4) 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 [70,71] for two different values of B 0 : In section 5.1.1 we start withB 0 = 0 (and a > 0), where the unstable modes have no preferred direction. The emerging pattern in this case is a square vortex lattice. In section 5.1.2, we then consider the caseB ⊥ 0 <B 0 <B 0 (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 figure 6(a) and (b), respectively. Various technical details related to the weakly nonlinear analysis are provided in Appendix E.

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 [70,71]). We introduce a small parameter ε characterizing the distance to the bifurcation via the growth rate of critical modes ε 2 = σ Re (B 0 = 0, k = 1) = a [compare equation (16) and (17)]. This allows the definition of a slow time scale T = ε 2 t and corresponding spatial variable X = εx (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 v (x, t, X, T ) = εA ⊥ (X, T )e ix ⊥ + c.c. , v ⊥ (x, t, X, T ) = εA (X, T )e ix + c.c. , 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 equation (19) into equation (15) 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 O(ε 3 ), where we already scaled back to the fast time and length scales, t and x, respectively. It is seen that the two amplitude equations (20) correspond to two coupled Ginzburg-Landau equations [72]. This result was already obtained in [61], where the pattern formation in the reduced model equation (15) in the absence of an external field is investigated. Physically, equations (20) describe a relaxation towards a uniform, stationary state characterized by the homogeneous stationary solution A sat = A sat ⊥ = a/(5b), corresponding to a regular square vortex lattice.

Finite external field
Being interested in the symmetry-broken state, we now move on to perform a weakly nonlinear analysis in the region of the state diagram (figure 5) where perturbations with a parallel wavevector grow but perpendicular wavevector perturbations decay. Starting from the homogeneous stationary state at large fieldsB 0 >B 0 , the onset of the instability (at a > 0) occurs atB 0 (a), i.e., the solid line in figure 5. In analogy to the analysis forB 0 = 0, we introduce a small parameter ε quantifying the distance to the bifurcation. In the present case, we define ε via where we have used equation (17). The slow time scale is again defined by T = ε 2 t. The scaling of the corresponding spatial variable is given by X = ε(x − v g t), where we introduced the group velocity v g for the following reason: In contrast to the caseB 0 = 0, the system atB 0 > 0 displays a net polarization, and the emerging stripe pattern [see figure 6(a)] 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 (19), 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 (E.12)], thus, the bifurcation atB 0 =B 0 is supercritical. As for the vortex lattice atB 0 = 0, the obtained amplitude equation (23) corresponds to a real-valued Ginzburg-Landau equation [72] and describes a relaxation process to a uniform amplitude, given by the homogeneous stationary solution A sat = σ Re /g. However, there are several features not present in equations (20): First, the stripe pattern and modulations of it are transported along the direction of the external field. This transport is characterized by the speed c and the group velocity v g = λV 0 , which is equal to the linear traveling speed c 0 . In addition to the result from linear stability analysis, the traveling speed c contains contributions of higher order in ε [see equation (E.9) in Appendix E]. Up to second order, we have 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 (20) and (23) 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 (15), 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 [compare equation (15)]. 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 (15) 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 forB ⊥ 0 <B 0 <B 0 [see figure 6(a) and (b), respectively]. Interestingly, when solving the dynamical equation (15) 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 strengthB * 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, (∇ × v) max = 2A + 2A ⊥ ] is plotted over the external field strength. Numerical results, obtained by solving equation (15), are denoted by red dots. The green lines are given by the stationary solutions of the amplitude equations for the square vortex lattice and the traveling stripes, equation (20) and (23), respectively. Interestingly, the transition field strengthB * 0 is not given byB ⊥ 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 thanB ⊥ 0 . The transition field strengthB * 0 is obtained  (20) and (23), respectively. The red dots are obtained by numerically solving the full dynamical equation (15). Decreasing the external field strength from the onset of the instability atB 0 we observe a stripe pattern that becomes unstable with respect to the formation of a square vortex lattice atB * 0 . Due to the suppression of perpendicular modes by parallel modes, the field strengthB ⊥ 0 obtained by linear stability analysis of the homogeneous stationary solution is not relevant.
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 (20), 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 and asymmetric. 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 asymmetric pattern formation, which we have analyzed by deriving the amplitude equations (20) and (23). 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, (20) and (23), 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 (23) 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 (25)]. 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 (23) 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 [72]. 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 [68]. 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 24]. 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 [73,74]. 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.

Appendix A. Relation to the microscopic model
A detailed derivation of the continuum equations in the absence of an external field, equations (1) -(4), is given in Ref. [19]. Extending the derivation towards an external field is quite straightforward. Here, we will give a short summary of the key points.
The overdamped motion of a swimmer σ in an ensemble of σ = 1, . . . , S identical swimmers is given by the Langevin equations for the position X σ and the orientation N σ , respectively, where the functions ξ σ and η σ 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 u(x, t). 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 Ω = 1 2 [(∇u) T − (∇u)] and deformation rate Σ = 1 2 [(∇u) T + (∇u)], respectively. In analogy to Jeffrey's theory for oscillatory tumbling motion of elongated particles in flow [75][76][77], 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 σ (where I is the unit matrix) is introduced to conserve the length of the unit vector N σ . 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, The pair potential Φ int (N µ , N ν , r µ,ν ) describes an activity-driven polar alignment of two swimmers with distance r µ,ν = |X µ − X ν | ≥ r, where γ 0 is the strength of the interaction, and Θ(x) = 1 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 [65], an effect that has been observed experimentally [3]. Finally, the external potential for swimmer µ is given by where the unit vector h denotes the field's direction. Here, we already introduced the potential in such a way that the field's magnitude 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 (4). 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 P(x, n, t), 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 ρ(x, t) = ρ and employ a mean-field approximation to treat the twoparticle correlations stemming from conservative interactions. We then project onto orientational moments n, nn, ... of P(x, n, t), which are directly connected to the polar order parameter P and the nematic order parameter Q via P = n, 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 [78]. 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 (6) that includes a coupling to the Stokes equation (4). The two coefficients c I and c F , not already specified in the main text, are given by 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 (6) or (15) using the Runge-Kutta-Fehlberg method (RKF45) [79]. 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 [80]. In the case of the full model, the Stokes equation (4) is solved applying a stream-function approach [80]. 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 (13) and (14)] determining the stability of the homogeneous stationary solution. To this end, we will first eliminate the dependence of the linearized system on the perturbations δû, δp and δq, thus reducing the set of variables to just the perturbation of the polar order parameter, δP.
As a first step we insert the perturbed solution [see equation (10) and (11) in the main text] into the Stokes equation (4) and linearize, yielding −|k| 2 δû = c F i6c I k P 0 δP − |k| 2 δP + 1 28 We multiply equation (C.1) by the wavevector k and employ the incompressibility conditions which, for the perturbed system, yield k · δû = 0 and k · δP = 0. From this, we find that the pressure perturbation amplitude must satisfy δp = 0. Equation (C.1) then yields the perturbation δû as a function of the perturbation of the polar order parameter, We now consider the equation of motion for the polar order parameter P, equation (6). Performing the linearization yields We can eliminate the dependence of equation (C.3) on the velocity perturbation δû by inserting equation (C.2). The reduced system can then be written as where the components of the Jacobian matrix M(k) are given by Equation (C.4) still involves the perturbation δq of the Lagrange multiplier. Utilizing again the incompressibility condition for the polar order parameter field P, yielding k · δP = 0, and equation (C.4), we find that the perturbation δq must satisfy (C.9) Inserting equation (C.9) into equation (C.4) we can identify the projector Π(k) = I − kk/|k| 2 . With this, we obtain equation (12) 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 (12)] via For the explicit form of σ(k), see equations (13) and (14) in the main text.

Appendix D. Coefficients in the reduced and rescaled model
Four dimensionless coefficients remain in the reduced and rescaled model equation (15). These are related to the coefficients of the full model given by equations (1), (2), (3) and (7) via Finally, the traveling speed c is also expanded in orders of ε, where the zeroth component is proportional to the homogeneous stationary solution V 0 , i.e., c 0 = λV 0 [see equation (18)]. Now, we insert all expansions [equations (22), (E. 3) and (E.2)] into equation (15) 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, O(ε 0 ), we recover equation (16) which determines the homogeneous stationary solution V 0 .
• In first order, O(ε 1 ), we find that the contributions v 1,0 , v 1,1 , v ⊥ 1,0 , q 1,0 , q 1,1 and c 1 must vanish to satisfy simultaneously the dynamical equation (15) and the incompressibility condition ∇ · v = 0. For v ⊥ 1,1 we obtain where we used the definition of ε according to equation (21). Thus, we find that v ⊥ 1,1 actually contributes in O(ε 3 ) to the amplitude equation obtained below. • In second order, O(ε 2 ), using the incompressibility constraint, we find the relation v 2,1 = i∂ X ⊥ v ⊥ 1,1 . (E.5) As a solvability condition for equation (15) we obtain for the zeroth mode (m = 0) v 2,0 = D 0 |v ⊥ 1,1 | 2 , (E. 6) where Matching all terms containing the first mode (m = 1) yields where we inserted equation (E.5). The coefficient D ⊥ is given in equation (25) in the main text. Equations (E.5) -(E.8) show that the second order amplitudes v 2,0 , v 2,1 and q 2,1 are "slaved" to the first order amplitude v ⊥ 1,1 . Further, we obtain the second order contribution to the traveling speed c 2 = λD 0 |v ⊥ 1,1 | 2 . (E.9) All other higher order modes (n > 1), i.e., v 2,2 , v ⊥ 2,0 , v ⊥ 2,1 , v ⊥ 2,2 , q 2,0 , q 2,2 , . . . either violate the incompressibility constraint or the solvability conditions obtained from equation (15), or they contribute in orders higher than O(ε 3 ).  Figure F1. Dominating mode k dom in the reduced system [equation (15)] 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 (15)] and determining the first maximum Λ max of the spatial correlation function v(x) · v(x + ∆x) (where . . . denotes the spatial and temporal average). The dominant mode then follows as k dom = 2π/Λ max . 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 c = 1. For larger λ, turbulent motion sets in and energy is transfered to larger scales, i.e., smaller wavenumbers.
• In third order, O(ε 3 ), we finally obtain a dynamic equation for the amplitude v ⊥ 1,1 by matching terms containing the first mode (m = 1), (E.10) Inserting equations (E.6) and (E.8) into equation (E.10) and replacing v ⊥ 1,1 → A yields The diffusion coefficients D and D ⊥ are given in equation (25) in the main text. Further, the coefficient of the cubic term follows as Finally, scaling back to the fast time and length scales, t and x, yields the amplitude equation (23) in the main text.
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 (15).

Appendix G. Stability of stripe pattern
The emerging pattern near the bifurcation atB 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 = σ Re /g of the amplitude equation (23) (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 where A * denotes the complex conjugate of A. The form of the perturbation we consider is a perpendicular mode with critical wavenumber k c = 1, i.e., Utilizing the incompressibility condition ∇ · v = 0, we find δv ⊥ = 0. We insert equation (G.2) into the dynamic equation for v [equation (15)] and obtain, after linearization, the growth rate of perpendicular modes with critical wavenumber k c = 1, The critical field strengthB * 0 defining the transition from the stripe pattern to the vortex lattice is then obtained by settingσ ⊥ = 0. Note thatσ ⊥ 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 thanB ⊥ 0 .