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

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–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–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–19]. For a selection of recent reviews on the full scope of active matter see [20–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–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–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–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–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.

Recently, a phenomenological model [47, 55–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 [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₀ 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 γ₀ 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₀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

$$\begin{eqnarray}&&{{\mathfrak{D}}}_{t}({\bf{u}}){\bf{P}}+{\lambda }_{0}{\bf{P}}\cdot {\rm{\nabla }}{\bf{P}}=-\displaystyle \frac{\delta {\mathfrak{F}}}{\delta {\bf{P}}},\qquad {\rm{\nabla }}\cdot {\bf{P}}=0.\end{eqnarray} \tag{ 1 }$$

The relaxation term is given as functional derivative of

$$\begin{eqnarray}&&{\mathfrak{F}}=q({\rm{\nabla }}\cdot {\bf{P}})+\displaystyle \frac{1}{2}\alpha | {\bf{P}}{| }^{2}+\displaystyle \frac{1}{4}\beta | {\bf{P}}{| }^{4}+\displaystyle \frac{1}{2}{{\rm{\Gamma }}}_{2}{\left({\rm{\nabla }}{\bf{P}}\right)}^{2}+\displaystyle \frac{1}{2}{{\rm{\Gamma }}}_{4}{\left({\rm{\nabla }}{\rm{\nabla }}{\bf{P}}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

such that,

$$\begin{eqnarray}&&\displaystyle \frac{\delta {\mathfrak{F}}}{\delta {\bf{P}}}={\rm{\nabla }}q+\alpha {\bf{P}}+\beta | {\bf{P}}{| }^{2}{\bf{P}}-{{\rm{\Gamma }}}_{2}{{\rm{\nabla }}}^{2}{\bf{P}}+{{\rm{\Gamma }}}_{4}{{\rm{\nabla }}}^{4}{\bf{P}}.\end{eqnarray} \tag{ 3 }$$

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 $\sqrt{-\alpha /\beta }$. For sufficiently strong activity Γ₂ becomes negative, which leads to a finite-wavelength instability of the homogeneous state, yielding a typical length scale of ${\rm{\Lambda }}=2\pi \sqrt{-2{{\rm{\Gamma }}}_{4}/{{\rm{\Gamma }}}_{2}}$. Turbulent dynamics is introduced to the model via the nonlinear advection term on the left-hand side of equation (1), with λ₀ giving the strength of the advection term. Finally, q(x) is a local Lagrange multiplier enforcing the incompressibility condition ${\rm{\nabla }}\cdot {\bf{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

$$\begin{eqnarray}&&{{\mathfrak{D}}}_{t}({\bf{u}}){\bf{P}}={\partial }_{t}{\bf{P}}+{\bf{u}}\cdot {\rm{\nabla }}{\bf{P}}-{\boldsymbol{\Omega }}\cdot {\bf{P}}-\kappa {\boldsymbol{\Sigma }}\cdot {\bf{P}},\end{eqnarray} \tag{ 4 }$$

where the vorticity tensor and the deformation rate are given by ${\boldsymbol{\Omega }}=\tfrac{1}{2}\left[{\left({\rm{\nabla }}{\bf{u}}\right)}^{{\rm{T}}}-({\rm{\nabla }}{\bf{u}})\right]$ and ${\boldsymbol{\Sigma }}=\tfrac{1}{2}\left[{\left({\rm{\nabla }}{\bf{u}}\right)}^{{\rm{T}}}+({\rm{\nabla }}{\bf{u}})\right]$, respectively. The modified Stokes equation (see [19]) that determines the solvent flow field u reads

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\bf{u}}={\rm{\nabla }}\cdot ({{\boldsymbol{\sigma }}}_{\mathrm{pass}}+{{\boldsymbol{\sigma }}}_{\mathrm{act}})+{\rm{\nabla }}p,\end{eqnarray} \tag{ 5 }$$

where p is an effective pressure and ${{\boldsymbol{\sigma }}}_{\mathrm{pass}}$ and ${{\boldsymbol{\sigma }}}_{\mathrm{act}}$ refer to the two parts of the stress tensor, σ, specifically, ${{\boldsymbol{\sigma }}}_{\mathrm{pass}}$ is an anisotropic passive contribution due to the elongated shape of the particles, and ${{\boldsymbol{\sigma }}}_{\mathrm{act}}$ 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 ${{\boldsymbol{\sigma }}}_{\mathrm{pass}}$ 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] yields

$$\begin{eqnarray}&&{\rm{\nabla }}{{\boldsymbol{\sigma }}}_{\mathrm{act}}={c}_{{\rm{F}}}\left(6{c}_{{\rm{I}}}{\bf{P}}\cdot {\rm{\nabla }}{\bf{P}}+{{\rm{\nabla }}}^{2}{\bf{P}}+\displaystyle \frac{1}{28}{{\rm{\nabla }}}^{4}{\bf{P}}\right).\end{eqnarray} \tag{ 6 }$$

Equations (1)–(6) are rescaled using the microswimmer length ℓ as length scale, the self-swimming speed v₀ as characteristic velocity and ℓ/v₀ as time scale [19]. The coefficients can then be written as

$$\begin{eqnarray}\begin{array}{rcl}{\lambda }_{0} & = & {c}_{{\rm{I}}}(3+2{a}_{0}{P}_{{\rm{r}}}{c}_{{\rm{F}}})/5,\qquad \kappa ={a}_{0}(3-{c}_{{\rm{I}}})/5,\\ \alpha & = & \left(1-{c}_{{\rm{I}}}\right){/P}_{{\rm{r}}},\qquad \beta =\displaystyle \frac{3}{5}{c}_{{\rm{I}}}^{2}{/P}_{{\rm{r}}},\\ {{\rm{\Gamma }}}_{2} & = & \displaystyle \frac{1}{10}{\left(\displaystyle \frac{r}{{\ell }}\right)}^{2}{c}_{{\rm{I}}}{/P}_{{\rm{r}}}-\displaystyle \frac{{a}_{0}}{15}{P}_{{\rm{r}}}{c}_{{\rm{F}}},\qquad {{\rm{\Gamma }}}_{4}=\displaystyle \frac{{a}_{0}}{420}{P}_{{\rm{r}}}{c}_{{\rm{F}}},\end{array}\end{eqnarray} \tag{ 7 }$$

where P_r, c_I, r/ℓ, c_F and a₀ are dimensionless parameters. The activity is quantified by the persistence number ${P}_{{\rm{r}}}={v}_{0}\tau /{\ell }$, 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}_{{\rm{I}}}\lt 1$ the system favors a disordered, isotropic state, while for ${c}_{{\rm{I}}}\gt 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₀ 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. ${{\rm{\Phi }}}_{\mathrm{ext}}\propto -{\bf{n}}\cdot {\bf{h}}$. 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–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)

$$\begin{eqnarray}&&{{\mathfrak{D}}}_{t}({\bf{u}}){\bf{P}}+{\lambda }_{0}{\bf{P}}\cdot {\rm{\nabla }}{\bf{P}}=-\displaystyle \frac{\delta {\mathfrak{F}}}{\delta {\bf{P}}}+{\bf{g}},\qquad {\rm{\nabla }}\cdot {\bf{P}}=0.\end{eqnarray} \tag{ 8 }$$

The additional term on the right-hand side of the evolution equation is given as

$$\begin{eqnarray}&&{\bf{g}}={B}_{0}{\bf{h}}\cdot \left[\displaystyle \frac{2}{3}{\bf{I}}-\displaystyle \frac{3}{5}{c}_{{\rm{I}}}\left({\bf{PP}}-\displaystyle \frac{{\bf{I}}({\bf{P}}\cdot {\bf{P}})}{3}\right)-\displaystyle \frac{2}{15}{a}_{0}{P}_{{\rm{r}}}{\boldsymbol{\Sigma }}\right],\end{eqnarray} \tag{ 9 }$$

where B₀ 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₀, 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 ${\mathfrak{F}}$ due to the closure relations used in the derivation, see appendix A and [19]. This indicates again that ${\mathfrak{F}}$ cannot be simply interpreted as a free energy.

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}_{\parallel }$ and ${x}_{\perp }$. 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 ${\left({\rm{\nabla }}\times {\bf{P}}\right)}_{z}$ (where z is both perpendicular to ${x}_{\parallel }$ and ${x}_{\perp }$) as scalar vorticity or just vorticity.

We start with the field-free case, B₀ = 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 λ₀ is rather large and, therefore, the shape of the vortices is highly irregular.

Zoom In Zoom Out Reset image size

Turning on the external field, B₀ > 0, introduces several new features. First, the external field induces a net polar order in the system, that is $\langle {\bf{P}}\rangle ={A}^{-1}\int {\rm{d}}{\bf{xP}}({\bf{x}},t)\ne 0$ (where A is the area). The swimmers' self-propulsion speed leads to a local transport in the direction of the polar order given by ${v}_{0}{\bf{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 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, 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.6 and B₀ = 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.

4.1. 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. ${\bf{P}}=({P}_{\parallel },{P}_{\perp })$. 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₀ and Lagrange-multiplier q = q₀ are set constant in space and time. Then, equation (8) reduces to

$$\begin{eqnarray}&&0=-\alpha {P}_{0}-\beta {P}_{0}^{3}+\displaystyle \frac{2}{3}{B}_{0}-\displaystyle \frac{2}{5}{B}_{0}{c}_{{\rm{I}}}{P}_{0}^{2}.\end{eqnarray} \tag{ 10 }$$

The real positive solution of equation (10) defines the homogeneous stationary state ${\bf{P}}=({P}_{0},0)$. This solution P₀ is plotted versus the external field strength B₀ in figure 2. As expected, the polar order grows upon increasing the field strength. Dividing by B₀ and evaluating the limit ${B}_{0}\to \infty$, equation (10) simplifies to

$$\begin{eqnarray}&&0=\displaystyle \frac{2}{3}-\displaystyle \frac{2}{5}{c}_{{\rm{I}}}{P}_{0}^{2}.\end{eqnarray} \tag{ 11 }$$

The solution of equation (11) defines the saturation value ${P}_{0}^{\mathrm{sat}}=\sqrt{5/(3{c}_{{\rm{I}}})}$ that is approached for ${B}_{0}\to \infty$.

Zoom In Zoom Out Reset image size

4.2. Linear stability analysis

In order to determine the linear stability of the homogeneous stationary solution ${\bf{P}}=({P}_{0},0)$, q = q₀, ${\bf{u}}=(0,0)$, p = p₀, we consider small perturbations, i.e.

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\parallel } & = & {P}_{0}+\delta {P}_{\parallel },\qquad {P}_{\perp }=\delta {P}_{\perp },\qquad q={q}_{0}+\delta q,\\ {u}_{\parallel } & = & \delta {u}_{\parallel },\qquad {u}_{\perp }=\delta {u}_{\perp },\qquad p={p}_{0}+\delta p.\end{array}\end{eqnarray} \tag{ 12 }$$

For all perturbations, we make the ansatz

$$\begin{eqnarray}&&(\delta {P}_{\parallel },\delta {P}_{\perp },\delta q,\delta {u}_{\parallel },\delta {u}_{\perp },\delta p)=(\delta {\hat{P}}_{\parallel },\delta {\hat{P}}_{\perp },\delta \hat{q},\delta {\hat{u}}_{\parallel },\delta {\hat{u}}_{\perp },\delta \hat{p})\ {{\rm{e}}}^{\sigma t+{\rm{i}}{\bf{k}}\cdot {\bf{x}}},\end{eqnarray} \tag{ 13 }$$

with wavevector ${\bf{k}}=({k}_{\parallel },{k}_{\perp })$. 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

$$\begin{eqnarray}&&\sigma \delta \hat{{\bf{P}}}={\boldsymbol{\Pi }}({\bf{k}})\cdot {\bf{M}}({\bf{k}})\cdot \delta \hat{{\bf{P}}},\end{eqnarray} \tag{ 14 }$$

where the projector ${\boldsymbol{\Pi }}({\bf{k}})={\bf{I}}-{\bf{kk}}/| {\bf{k}}{| }^{2}$ arises as a consequence of the incompressibility of the order parameter field, and the 2 × 2-matrix ${\bf{M}}({\bf{k}})$ is the Jacobian. The components of ${\bf{M}}({\bf{k}})$ are given in equation (C.5) in appendix C. We obtain the complex growth rate $\sigma ={\sigma }_{\mathrm{Re}}+{\rm{i}}{\sigma }_{\mathrm{Im}}$ as a function of the wavevector k by calculating the eigenvalues of the matrix ${\boldsymbol{\Pi }}({\bf{k}})\cdot {\bf{M}}({\bf{k}})$. The real part ${\sigma }_{\mathrm{Re}}({\bf{k}})$, which determines the actual growth of a mode ${\bf{k}}$, is given by

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{\mathrm{Re}}({\bf{k}}) & = & -\alpha -{{\rm{\Gamma }}}_{2}| {\bf{k}}{| }^{2}-{{\rm{\Gamma }}}_{4}| {\bf{k}}{| }^{4}-\beta {P}_{0}^{2}\left({k}_{\parallel }^{2}+3{k}_{\perp }^{2}\right)/| {\bf{k}}{| }^{2}\\ & & +3{c}_{{\rm{F}}}{c}_{{\rm{I}}}{P}_{0}\left({P}_{0}(1+\kappa )-\displaystyle \frac{2}{15}{a}_{0}{P}_{{\rm{r}}}{B}_{0}\right){k}_{\parallel }^{2}/| {\bf{k}}{| }^{2}\\ & & -\displaystyle \frac{{c}_{{\rm{I}}}}{5}{B}_{0}{P}_{0}\left(3{k}_{\parallel }^{2}+4{k}_{\perp }^{2}\right)/| {\bf{k}}{| }^{2}.\end{array}\end{eqnarray} \tag{ 15 }$$

The imaginary part ${\sigma }_{\mathrm{Im}}({\bf{k}})$, which determines the linear traveling speed ${c}_{0}^{\parallel }$, is given by

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{\mathrm{Im}}({\bf{k}}) & = & \left[-{\lambda }_{0}{P}_{0}+\displaystyle \frac{{c}_{{\rm{F}}}}{2}\left({P}_{0}(1+\kappa )-\displaystyle \frac{2}{15}{B}_{0}{a}_{0}{P}_{{\rm{r}}}\right)\left(1-\displaystyle \frac{| {\bf{k}}{| }^{2}}{28}\right)\right]{k}_{\parallel }\\ & = & -{c}_{0}^{\parallel }{k}_{\parallel }.\end{array}\end{eqnarray} \tag{ 16 }$$

It is seen that ${\sigma }_{\mathrm{Re}}({\bf{k}})$ depends not only on the magnitude $| {\bf{k}}| =k$ of the wavevector but also on its direction. For B₀ = 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 }$ parallel to the direction of spontaneous order and ${x}_{\perp }$ perpendicular to that direction.

4.3. State diagram

Due to the explicit dependence of the growth rate ${\sigma }_{\mathrm{Re}}({\bf{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}_{\parallel }=k$, ${k}_{\perp }=0$, or purely perpendicular to the field, i.e. ${k}_{\parallel }=0$, ${k}_{\perp }=k$. The incompressibility condition of the order parameter field, i.e. ${k}_{\parallel }\delta {P}_{\parallel }+{k}_{\perp }\delta {P}_{\perp }=0$, then dictates the form of the respective perturbations. For a parallel wavevector, the component $\delta {P}_{\parallel }$ 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 ${{\rm{\Lambda }}}_{\parallel }=2\pi /k$. This type of pattern manifests itself as stripes in the vorticity of the field. Analogously, for a perpendicular wavevector, the component $\delta {P}_{\perp }$ must vanish and the pattern is a perturbation along the field with periodicity in the perpendicular direction (${{\rm{\Lambda }}}_{\perp }=2\pi /k$) (see figure 3(b)). The occurrence of both, modes in the parallel and in the perpendicular direction, results in a square vortex lattice.

Zoom In Zoom Out Reset image size

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₀, 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, we observe an instability at a critical value (${P}_{{\rm{r}}}\approx 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 ${\rm{\Lambda }}=2\pi \sqrt{-2{{\rm{\Gamma }}}_{4}/{{\rm{\Gamma }}}_{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 ${\bf{P}}$ and its vorticity ${\rm{\nabla }}\times {\bf{P}}$ in the mesoscale turbulent state at B₀ = 0, see figure 1(a).

Zoom In Zoom Out Reset image size

Keeping the persistence number at a constant value within the zero-field turbulent state, e.g. ${P}_{{\rm{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}^{\perp }({P}_{{\rm{r}}})$, see dashed line in figure 4(a). In contrast, perturbations with parallel wavevector still grow until a larger field strength of ${B}_{0}^{\parallel }\,({P}_{{\rm{r}}})$ is reached, see solid line in figure 4(a). Between the two 'critical' field strengths ${B}_{0}^{\parallel }$ and ${B}_{0}^{\perp }$, the growth of parallel modes leads to anisotropic patterns which are elongated in the perpendicular direction. This is consistent with our numerical observations shown in figures 1(b) and (c). Interestingly, the characteristic length given by the critical mode k_c, i.e. the maximum of ${\sigma }_{\mathrm{Re}}({\bf{k}})$, is not influenced by the external field (compare figures 4(c)–(e)).

Moreover, we observe a very intriguing feature in the case of a perturbation with ${\bf{k}}\parallel {\bf{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 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}_{{\rm{r}}}\approx 5$.

Finally, note that the growth rate ${\sigma }_{\mathrm{Re}}$ (see equation (15)) does not depend on the nonlinear advection term in equation (1). The imaginary part ${\sigma }_{\mathrm{Im}}$, however, is strongly dependent on to the parameter λ₀ 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}^{\parallel }$, 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).

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}_{{\rm{F}}}\ll 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 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–59, 61]. Note that the coefficient Γ₂ has to be negative in order to obtain a mesoscale-turbulent state. Thus, additionally to the limit ${c}_{{\rm{F}}}\ll 1$, the product ${P}_{{\rm{r}}}{c}_{{\rm{F}}}$ still has to be sufficiently large compared to c_I (see equation (7)). 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}_{{\rm{c}}}=\sqrt{-{{\rm{\Gamma }}}_{2}/(2{{\rm{\Gamma }}}_{4})}$, and the inverse of the corresponding maximum of the growth rate ${\sigma }_{\mathrm{Re}}(\alpha =0,{B}_{0}=0)={{\rm{\Gamma }}}_{2}^{2}/(4{{\rm{\Gamma }}}_{4})$. Further, we rescale by the saturation value for the polar order parameter at ${B}_{0}\to \infty$ and obtain the scaled collective microswimmer velocity via ${\bf{v}}={\bf{P}}/{P}_{0}^{\mathrm{sat}}$, where ${P}_{0}^{\mathrm{sat}}=\sqrt{5/(3{c}_{{\rm{I}}})}$. The rescaled equation for v is then given by

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}{\bf{v}}+\lambda {\bf{v}}\cdot {\rm{\nabla }}{\bf{v}} & = & -{\rm{\nabla }}q+a{\bf{v}}-b| {\bf{v}}{| }^{2}{\bf{v}}-{\left(1+{{\rm{\nabla }}}^{2}\right)}^{2}{\bf{v}}\\ & & +{\tilde{B}}_{0}{\bf{h}}\cdot \left[{\bf{I}}-\displaystyle \frac{1}{2}\left(3{\bf{vv}}-{\bf{I}}({\bf{v}}\cdot {\bf{v}})\right)\right].\end{array}\end{eqnarray} \tag{ 17 }$$

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 ${\tilde{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 ${\tilde{B}}_{0}$. The scaling of the Lagrange multiplier q enforcing the incompressibility ${\rm{\nabla }}\cdot {\bf{v}}=0$ is irrelevant, thus, we keep the same notation as before. For the explicit dependencies of the four remaining coefficients λ, a, b and ${\tilde{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. The homogeneous stationary solution ${\bf{v}}=({v}_{\parallel },{v}_{\perp })=({V}_{0},0)$ of the reduced model is calculated via

$$\begin{eqnarray}&&0=(a-1){V}_{0}-{{bV}}_{0}^{3}+{\tilde{B}}_{0}-{\tilde{B}}_{0}{V}_{0}^{2}.\end{eqnarray} \tag{ 18 }$$

Analogous to the full model, performing a linear stability analysis yields the complex growth rate with real part

$$\begin{eqnarray}&&{\sigma }_{\mathrm{Re}}({\bf{k}})=a-1+2| {\bf{k}}{| }^{2}-| {\bf{k}}{| }^{4}-{{bV}}_{0}^{2}\displaystyle \frac{{k}_{\parallel }^{2}+3{k}_{\perp }^{2}}{| {\bf{k}}{| }^{2}}-\displaystyle \frac{1}{2}{\tilde{B}}_{0}{V}_{0}\displaystyle \frac{3{k}_{\parallel }^{2}+4{k}_{\perp }^{2}}{| {\bf{k}}{| }^{2}}.\end{eqnarray} \tag{ 19 }$$

Similar to the corresponding function ${\sigma }_{{\bf{k}}}$ of the full model (compare equation (15)), equation (19) 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

$$\begin{eqnarray}&&{c}_{0}^{\parallel }=-\displaystyle \frac{{\sigma }_{\mathrm{Im}}({\bf{k}})}{{k}_{\parallel }}=\lambda {V}_{0}.\end{eqnarray} \tag{ 20 }$$

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.

Zoom In Zoom Out Reset image size

5.1. Weakly nonlinear analysis

To analyze in more detail the emerging patterns in the reduced model, we now perform a weakly nonlinear analysis [78, 79] for two different values of B₀: in section 5.1.1 we start with ${\tilde{B}}_{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 case ${\tilde{B}}_{0}^{\perp }\lt {\tilde{B}}_{0}\lt {\tilde{B}}_{0}^{\parallel }$ (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).

Zoom In Zoom Out Reset image size

5.1.1. 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 ${\varepsilon }^{2}={\sigma }_{\mathrm{Re}}({B}_{0}=0,k=1)=a$ (compare equations (18) and (19)). This allows the definition of a slow time scale $T={\varepsilon }^{2}t$ and corresponding spatial variable ${\bf{X}}=\varepsilon {\bf{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, 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}_{\parallel }$ and ${v}_{\perp }$ of the effective velocity field. Taking into account that the incompressibility condition, i.e. ${\rm{\nabla }}\cdot {\bf{v}}=0$, has to be satisfied, the expansion reduces in lowest order to

$$\begin{eqnarray}\begin{array}{rcl}{v}_{\parallel }({\bf{x}},t,{\bf{X}},T) & = & \varepsilon {A}_{\perp }({\bf{X}},T){{\rm{e}}}^{{\rm{i}}{x}_{\perp }}+{\rm{c}}.{\rm{c}}.,\\ {v}_{\perp }({\bf{x}},t,{\bf{X}},T) & = & \varepsilon {A}_{\parallel }({\bf{X}},T){{\rm{e}}}^{{\rm{i}}{x}_{\parallel }}+{\rm{c}}.{\rm{c}}.,\end{array}\end{eqnarray} \tag{ 21 }$$

where ${A}_{\parallel }$ and ${A}_{\perp }$ are the amplitudes of the two modes and ${\rm{c}}.{\rm{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 ${ \mathcal O }({\varepsilon }^{3})$,

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}{A}_{\parallel } & = & {{aA}}_{\parallel }-b\left(3| {A}_{\parallel }{| }^{2}+2| {A}_{\perp }{| }^{2}\right){A}_{\parallel }+4{\partial }_{{x}_{\parallel }}^{2}{A}_{\parallel },\\ {\partial }_{t}{A}_{\perp } & = & {{aA}}_{\perp }-b\left(3| {A}_{\perp }{| }^{2}+2| {A}_{\parallel }{| }^{2}\right){A}_{\perp }+4{\partial }_{{x}_{\perp }}^{2}{A}_{\perp },\end{array}\end{eqnarray} \tag{ 22 }$$

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 (17) in the absence of an external field is investigated. Physically, equations (22) describe a relaxation towards a uniform, stationary state characterized by the homogeneous stationary solution ${A}_{\parallel }^{\mathrm{sat}}={A}_{\perp }^{\mathrm{sat}}=\sqrt{a/(5b)}$, corresponding to a regular square vortex lattice.

5.1.2. 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 fields ${\tilde{B}}_{0}\gt {\tilde{B}}_{0}^{\parallel }$, the onset of the instability (at a > 0) occurs at ${\tilde{B}}_{0}^{\parallel }\,(a)$, i.e. the solid line in figure 5. In analogy to the analysis for ${\tilde{B}}_{0}=0$, we introduce a small parameter ε quantifying the distance to the bifurcation. In the present case, we define ε via

$$\begin{eqnarray}&&{\varepsilon }^{2}={\sigma }_{\mathrm{Re}}^{\parallel }={\sigma }_{\mathrm{Re}}({k}_{\parallel }=1,{k}_{\perp }=0)=a-{{bV}}_{0}^{2}-\displaystyle \frac{3}{2}{\tilde{B}}_{0}{V}_{0},\end{eqnarray} \tag{ 23 }$$

where we have used equation (19). The slow time scale is again defined by $T={\varepsilon }^{2}t$. The scaling of the corresponding spatial variable is given by ${\bf{X}}=\varepsilon ({\bf{x}}-{{\bf{v}}}_{{\rm{g}}}t)$, where we introduced the group velocity v_g for the following reason: in contrast to the case ${\tilde{B}}_{0}=0$, the system at ${\tilde{B}}_{0}\gt 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}_{\parallel }$. 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

$$\begin{eqnarray}\begin{array}{rcl}{v}_{\parallel }({\bf{x}},t,{\bf{X}},T) & = & {V}_{0}+\displaystyle \sum _{n=1}^{\infty }\displaystyle \sum _{m=0}^{n}{\varepsilon }^{n}{v}_{n,m}^{\parallel }({\bf{X}},T){{\rm{e}}}^{{\rm{i}}m({x}_{\parallel }-{c}_{\parallel }t)}+{\rm{c}}.{\rm{c}}.,\\ {v}_{\perp }({\bf{x}},t,{\bf{X}},T) & = & \displaystyle \sum _{n=1}^{\infty }\displaystyle \sum _{m=0}^{n}{\varepsilon }^{n}{v}_{n,m}^{\perp }({\bf{X}},T){{\rm{e}}}^{{\rm{i}}m({x}_{\parallel }-{c}_{\parallel }t)}+{\rm{c}}.{\rm{c}}.,\end{array}\end{eqnarray} \tag{ 24 }$$

where V₀ denotes the homogeneous stationary solution $[{\bf{v}}=({v}_{\parallel },{v}_{\perp })=({V}_{0},0)]$ and ${\rm{c}}.{\rm{c}}.$ the complex conjugate. Similarly, we also expand the local Lagrange multiplier q and traveling speed ${c}^{\parallel }$ (see equations (E.2) and (E.3) in appendix E). Inserting all expansions (equations (24), (E.3) and (E.2)) into the dynamic equation (17) and using the incompressibility constraint ${\rm{\nabla }}\cdot {\bf{v}}=0$, one obtains solvability conditions in different modes ${{\rm{e}}}^{{\rm{i}}m({x}_{\parallel }-{c}_{\parallel }t)}$ and different orders of ε that all have to be satisfied. In ${ \mathcal O }({\varepsilon }^{3})$ we obtain a dynamic equation for the leading order amplitude ${v}_{1,1}^{\perp }$, which we will denote as A in what follows for the sake of brevity. All other contributions in the expansion either vanish, contribute in orders higher than ${ \mathcal O }({\varepsilon }^{3})$ or can be expressed as functions of A, i.e. are slaved to the dominating mode. For further technical details of the weakly nonlinear analysis, see appendix E. The final amplitude equation for the dominating mode is

$$\begin{eqnarray}&&{\partial }_{t}A+{v}_{{\rm{g}}}^{\parallel }{\partial }_{{x}_{\parallel }}A={\sigma }_{\mathrm{Re}}^{\parallel }A-g| A{| }^{2}A+{D}_{\parallel }{\partial }_{{x}_{\parallel }}^{2}A+{D}_{\perp }{\partial }_{{x}_{\perp }}^{2}A,\end{eqnarray} \tag{ 25 }$$

where we already scaled back to the fast scales.

We find that the coefficient g > 0 (see equation (E.12)), thus, the bifurcation at ${\tilde{B}}_{0}={\tilde{B}}_{0}^{\parallel }$ is supercritical. As for the vortex lattice at ${\tilde{B}}_{0}=0$, the obtained amplitude equation (25) corresponds to a real-valued Ginzburg–Landau equation [80] and describes a relaxation process to a uniform amplitude, given by the homogeneous stationary solution ${A}^{\mathrm{sat}}=\sqrt{{\sigma }_{\mathrm{Re}}^{\parallel }/g}$. However, there are several features not present in equations (22): 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}^{\parallel }$ and the group velocity ${v}_{{\rm{g}}}^{\parallel }=\lambda {V}_{0}$, which is equal to the linear traveling speed ${c}_{0}^{\parallel }$. In addition to the result from linear stability analysis, the traveling speed ${c}^{\parallel }$ contains contributions of higher order in ε (see equation (E.9) in appendix E). Up to second order, we have

$$\begin{eqnarray}&&{c}^{\parallel }=\lambda ({V}_{0}+{D}_{0}| A{| }^{2}).\end{eqnarray} \tag{ 26 }$$

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}_{\parallel }$ and ${D}_{\perp }$, respectively,

$$\begin{eqnarray}&&{D}_{\parallel }=4,\qquad {D}_{\perp }=2{{bV}}_{0}^{2}+\displaystyle \frac{1}{2}{\tilde{B}}_{0}{V}_{0}.\end{eqnarray} \tag{ 27 }$$

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.

5.2. 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 $\langle {\bf{v}}\rangle$ and can be written as $\lambda \langle {\bf{v}}\rangle \cdot {\rm{\nabla }}{\bf{v}}$ (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 ${\tilde{B}}_{0}=0$ and traveling stripes for ${\tilde{B}}_{0}^{\perp }\lt {\tilde{B}}_{0}\lt {\tilde{B}}_{0}^{\parallel }$ (see figures 6(a) and (b), respectively) Interestingly, when solving the dynamical equation (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. ${\tilde{B}}_{0}\gt 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 ${\tilde{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, $({\rm{\nabla }}\times {\bf{v}}{)}_{z}^{\max }=2{A}_{\parallel }+2{A}_{\perp }$) is plotted over the external field strength. Numerical results, obtained by solving the full dynamical equation for v (equation (17)), are denoted by red dots. The lines denote the stationary solutions of the amplitude equations (22) and (25), respectively. In detail, the blue line indicates the amplitude of the vortex lattice, valid at ${\tilde{B}}_{0}=0$, and the green line denotes the amplitude of the stripe pattern. Note that the stationary solutions of the amplitude equations are not solutions of the dynamical equation for v (equation (17)) in a strict mathematical sense, but rather approximate solutions that are obtained by expansions around the bifurcations at a = 0 and ${\tilde{B}}_{0}={\tilde{B}}_{0}^{\parallel }$, respectively. In this sense, the dashed green line does not indicate an unstable solution of the dynamical equation (17), but rather denotes the solution of the amplitude equation (25) for the stripe pattern, however, in a region (${\tilde{B}}_{0}\lt {\tilde{B}}_{0}^{* }$) where this equation is no longer valid (because perpendicular modes start to grow).

Zoom In Zoom Out Reset image size

Interestingly, the transition field strength ${\tilde{B}}_{0}^{* }$ is not given by ${\tilde{B}}_{0}^{\perp }$, 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 ${\tilde{B}}_{0}^{\perp }$. The transition field strength ${\tilde{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.

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}_{\parallel }$ and ${D}_{\perp }$. 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}_{\parallel }/{D}_{\perp }$ 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.

We thank the Deutsche Forschungsgemeinschaft for financial support through GRK 1558 and SFB 910 (projects B2 and B5), HE 5995/3 and BA 1222/7.

A detailed derivation of the continuum equations in the absence of an external field, equations (1)–(6), is given in [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 ${{\bf{X}}}^{\sigma }$ and the orientation ${{\bf{N}}}^{\sigma }$, respectively,

$$\begin{eqnarray}&&{\dot{{\bf{X}}}}^{\sigma }={v}_{0}{{\bf{N}}}^{\sigma }+{\bf{u}}({{\bf{X}}}^{\sigma })+\sqrt{2D}{{\boldsymbol{\xi }}}^{\sigma },\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}&&{\dot{{\bf{N}}}}^{\sigma }={\boldsymbol{\Omega }}({{\bf{X}}}^{\sigma })\cdot {{\bf{N}}}^{\sigma }+{\boldsymbol{\Pi }}\cdot \left[{a}_{0}{\boldsymbol{\Sigma }}({{\bf{X}}}^{\sigma })\cdot {{\bf{N}}}^{\sigma }-{{\rm{\nabla }}}_{{{\bf{N}}}^{\sigma }}{\rm{\Phi }}+\sqrt{1/\tau }{{\boldsymbol{\eta }}}^{\sigma }\right],\end{eqnarray} \tag{ A.2 }$$

where the functions ${{\boldsymbol{\xi }}}^{\sigma }$ and ${{\boldsymbol{\eta }}}^{\sigma }$ 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₀. Additionally, the swimmer is transported by advection with the averaged surrounding flow field ${\bf{u}}({\bf{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 ${\boldsymbol{\Omega }}=\tfrac{1}{2}\left[{\left({\rm{\nabla }}{\bf{u}}\right)}^{{\rm{T}}}-({\rm{\nabla }}{\bf{u}})\right]$ and deformation rate ${\boldsymbol{\Sigma }}=\tfrac{1}{2}\left[{\left({\rm{\nabla }}{\bf{u}}\right)}^{{\rm{T}}}+({\rm{\nabla }}{\bf{u}})\right]$, respectively. In analogy to Jeffrey's theory for oscillatory tumbling motion of elongated particles in flow [83–85], the shape parameter a₀ is given as a function of the aspect ratio, defined by length ℓand diameter d of the swimmers,

$$\begin{eqnarray}&&{a}_{0}=\displaystyle \frac{{\left({\ell }/d\right)}^{2}-1}{{\left({\ell }/d\right)}^{2}+1}.\end{eqnarray} \tag{ A.3 }$$

Further, the projector ${\boldsymbol{\Pi }}={\bf{I}}-{{\bf{N}}}^{\sigma }{{\bf{N}}}^{\sigma }$ (where ${\bf{I}}$ is the unit matrix) is introduced to conserve the length of the unit vector ${{\bf{N}}}^{\sigma }$. 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,

$$\begin{eqnarray}&&{\rm{\Phi }}=\displaystyle \sum _{\mu ,\nu }{{\rm{\Phi }}}_{\mathrm{int}}({{\bf{N}}}^{\mu },{{\bf{N}}}^{\nu },{r}_{\mu ,\nu })+\displaystyle \sum _{\mu }{{\rm{\Phi }}}_{\mathrm{ext}}({{\bf{N}}}^{\mu }).\end{eqnarray} \tag{ A.4 }$$

The pair potential ${{\rm{\Phi }}}_{\mathrm{int}}({{\bf{N}}}^{\mu },{{\bf{N}}}^{\nu },{r}_{\mu ,\nu })$ describes an activity-driven polar alignment of two swimmers with distance ${r}_{\mu ,\nu }=| {{\bf{X}}}_{\mu }-{{\bf{X}}}_{\nu }| \geqslant r$,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{int}}({{\bf{N}}}^{\mu },{{\bf{N}}}^{\nu },{r}_{\mu ,\nu })=-\displaystyle \frac{{\gamma }_{0}{v}_{0}}{2}{{\bf{N}}}^{\mu }\cdot {{\bf{N}}}^{\nu }{\rm{\Theta }}(r-{r}_{\mu ,\nu }),\end{eqnarray} \tag{ A.5 }$$

where ${\gamma }_{0}$ is the strength of the interaction, and ${\rm{\Theta }}(x)=1$ for $x\geqslant 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

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{ext}}({{\bf{N}}}^{\mu })=-{B}_{0}\displaystyle \frac{{v}_{0}}{{\ell }}{\bf{h}}\cdot {{\bf{N}}}^{\mu },\end{eqnarray} \tag{ A.6 }$$

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 relative to the active time scale ℓ/v₀ defines the dimensionless external field strength B₀. The dimension of the full prefactor ${B}_{0}{v}_{0}/{\ell }$ appearing in equation (A.6) then is an inverse time due to the scaling of the Langevin equations (A.1) and (A.2).

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 ${ \mathcal P }({\bf{x}},{\bf{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 $\rho ({\bf{x}},t)=\rho$ and employ a mean-field approximation to treat the two-particle correlations stemming from conservative interactions. We then project on to orientational moments $\overline{{\bf{n}}}$, $\overline{{\bf{nn}}}$, ... of ${ \mathcal P }({\bf{x}},{\bf{n}},t)$, which are directly connected to the polar order parameter ${\bf{P}}$ and the nematic order parameter ${\bf{Q}}$ via

$$\begin{eqnarray}&&{\bf{P}}=\overline{{\bf{n}}},\qquad {\bf{Q}}=\overline{{\bf{nn}}}-{\bf{I}}/3.\end{eqnarray} \tag{ A.7 }$$

Following our previous studies [18, 19] we apply a closure relation for the nematic order parameter,

$$\begin{eqnarray}&&{\bf{Q}}=q({\bf{PP}}-({\bf{P}}\cdot {\bf{P}}){\bf{I}}/3)+{\lambda }_{{\rm{K}}}{\boldsymbol{\Sigma }},\end{eqnarray} \tag{ A.8 }$$

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) that includes a coupling to the Stokes equation (6). The two coefficients c_I and c_F, not already specified in the main text, are given by

$$\begin{eqnarray}&&{c}_{{\rm{I}}}=\displaystyle \frac{8}{9}\tau \rho {\gamma }_{0}{v}_{0}{r}^{3},\qquad {c}_{{\rm{F}}}=\displaystyle \frac{{f}_{0}\rho {{\ell }}^{2}}{10\mu {v}_{0}},\end{eqnarray} \tag{ A.9 }$$

where f₀ denotes the strength of the force dipole exerted by every swimmer on the surrounding fluid, and μ is the effective viscosity of the suspension [19].

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.

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 $\delta \hat{{\bf{u}}}$, $\delta \hat{p}$ and $\delta \hat{q}$, thus reducing the set of variables to just the perturbation of the polar order parameter, $\delta \hat{{\bf{P}}}$.

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

$$\begin{eqnarray}&&-| {\bf{k}}{| }^{2}\delta \hat{{\bf{u}}}={c}_{{\rm{F}}}\left({\rm{i}}6{c}_{{\rm{I}}}{k}_{\parallel }{P}_{0}\delta \hat{{\bf{P}}}-| {\bf{k}}{| }^{2}\delta \hat{{\bf{P}}}+\displaystyle \frac{1}{28}| {\bf{k}}{| }^{4}\delta \hat{{\bf{P}}}\right)+{\rm{i}}{\bf{k}}\delta \hat{p}.\end{eqnarray} \tag{ C.1 }$$

We multiply equation (C.1) by the wavevector k and employ the incompressibility conditions which, for the perturbed system, yield ${\bf{k}}\cdot \delta \hat{{\bf{u}}}=0$ and ${\bf{k}}\cdot \delta \hat{{\bf{P}}}=0$. From this, we find that the pressure perturbation amplitude must satisfy $\delta \hat{p}=0$. Equation (C.1) then yields the perturbation $\delta \hat{{\bf{u}}}$ as a function of the perturbation of the polar order parameter,

$$\begin{eqnarray}&&\delta \hat{{\bf{u}}}={c}_{{\rm{F}}}\left(-{\rm{i}}\displaystyle \frac{6{c}_{{\rm{I}}}{k}_{\parallel }{P}_{0}}{| {\bf{k}}{| }^{2}}+1-\displaystyle \frac{1}{28}| {\bf{k}}{| }^{2}\right)\delta \hat{{\bf{P}}}.\end{eqnarray} \tag{ C.2 }$$

We now consider the equation of motion for the polar order parameter P, equation (8). Performing the linearization yields

$$\begin{eqnarray}\begin{array}{rcl}\sigma \delta {\hat{P}}_{\parallel } & = & -{\rm{i}}{\lambda }_{0}{k}_{\parallel }{P}_{0}\delta {\hat{P}}_{\parallel }-\alpha \delta {\hat{P}}_{\parallel }-3\beta {P}_{0}^{2}\delta {\hat{P}}_{\parallel }-{{\rm{\Gamma }}}_{2}| {\bf{k}}{| }^{2}\delta {\hat{P}}_{\parallel }-{{\rm{\Gamma }}}_{4}| {\bf{k}}{| }^{4}\delta {\hat{P}}_{\parallel }\\ & & -\displaystyle \frac{4}{5}{c}_{{\rm{I}}}{B}_{0}{P}_{0}\delta {\hat{P}}_{\parallel }-{\rm{i}}\displaystyle \frac{2}{15}{a}_{0}{P}_{{\rm{r}}}{B}_{0}{k}_{\parallel }\delta {\hat{u}}_{\parallel }-{\rm{i}}{k}_{\parallel }\delta \hat{q}+{\rm{i}}\kappa {k}_{\parallel }\delta {\hat{u}}_{\parallel }{P}_{0},\\ \sigma \delta {\hat{P}}_{\perp } & = & -{\rm{i}}{\lambda }_{0}{k}_{\parallel }{P}_{0}\delta {\hat{P}}_{\perp }-\alpha \delta {\hat{P}}_{\perp }-\beta {P}_{0}^{2}\delta {\hat{P}}_{\perp }-{{\rm{\Gamma }}}_{2}| {\bf{k}}{| }^{2}\delta {\hat{P}}_{\perp }-{{\rm{\Gamma }}}_{4}| {\bf{k}}{| }^{4}\delta {\hat{P}}_{\perp }\\ & & -\displaystyle \frac{3}{5}{c}_{{\rm{I}}}{B}_{0}{P}_{0}\delta {\hat{P}}_{\perp }-{\rm{i}}\displaystyle \frac{1}{15}{a}_{0}{P}_{{\rm{r}}}{B}_{0}({k}_{\parallel }\delta {\hat{u}}_{\perp }+{k}_{\perp }\delta {\hat{u}}_{\parallel })-{\rm{i}}{k}_{\perp }\delta \hat{q}\\ & & +{\rm{i}}\displaystyle \frac{1}{2}\left({k}_{\parallel }\delta {\hat{u}}_{\perp }-{k}_{\perp }\delta {\hat{u}}_{\parallel }\right){P}_{0}+{\rm{i}}\displaystyle \frac{1}{2}\kappa \left({k}_{\parallel }\delta {\hat{u}}_{\perp }+{k}_{\perp }\delta {\hat{u}}_{\parallel }\right){P}_{0}.\end{array}\end{eqnarray} \tag{ C.3 }$$

We can eliminate the dependence of equation (C.3) on the velocity perturbation $\delta \hat{{\bf{u}}}$ by inserting equation (C.2). The reduced system can then be written as

$$\begin{eqnarray}&&\sigma \delta \hat{{\bf{P}}}=-{\rm{i}}{\bf{k}}\delta \hat{q}+{\bf{M}}({\bf{k}})\cdot \delta \hat{{\bf{P}}},\end{eqnarray} \tag{ C.4 }$$

where the components of the Jacobian matrix M(k) are given by

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\parallel \parallel }({\bf{k}}) & = & -\alpha -3\beta {P}_{0}^{2}-{{\rm{\Gamma }}}_{2}| {\bf{k}}{| }^{2}-{{\rm{\Gamma }}}_{4}| {\bf{k}}{| }^{4}-{\rm{i}}{\lambda }_{0}{k}_{\parallel }{P}_{0}-\displaystyle \frac{4}{5}{B}_{0}{c}_{{\rm{I}}}{P}_{0}\\ & & -\displaystyle \frac{4}{5}{P}_{{\rm{r}}}{c}_{{\rm{F}}}{c}_{{\rm{I}}}{a}_{0}{B}_{0}\displaystyle \frac{{k}_{\parallel }^{2}}{| {\bf{k}}{| }^{2}}{P}_{0}-{\rm{i}}\displaystyle \frac{2}{15}{a}_{0}{B}_{0}{c}_{{\rm{F}}}{P}_{{\rm{r}}}{k}_{\parallel }\left(1-| {\bf{k}}{| }^{2}/28\right)\\ & & +6\kappa {c}_{{\rm{I}}}{c}_{{\rm{F}}}\displaystyle \frac{{k}_{\parallel }^{2}}{| {\bf{k}}{| }^{2}}{P}_{0}^{2}+{\rm{i}}\kappa {c}_{{\rm{F}}}{k}_{\parallel }\left(1-| {\bf{k}}{| }^{2}/28\right){P}_{0},\end{array}\end{eqnarray} \tag{ C.5 }$$

$$\begin{eqnarray}&&{M}_{\parallel \perp }({\bf{k}})=\ 0,\end{eqnarray} \tag{ C.6 }$$

$$\begin{eqnarray}&&{M}_{\perp \parallel }({\bf{k}})=\ 0,\end{eqnarray} \tag{ C.7 }$$

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\perp \perp }({\bf{k}}) & = & -\alpha -\beta {P}_{0}^{2}-{{\rm{\Gamma }}}_{2}| {\bf{k}}{| }^{2}-{{\rm{\Gamma }}}_{4}| {\bf{k}}{| }^{4}-{\rm{i}}{\lambda }_{0}{k}_{\parallel }{P}_{0}-\displaystyle \frac{3}{5}{B}_{0}{c}_{{\rm{I}}}{P}_{0}\\ & & -\displaystyle \frac{2}{5}{P}_{{\rm{r}}}{c}_{{\rm{F}}}{c}_{{\rm{I}}}{a}_{0}{B}_{0}\displaystyle \frac{{k}_{\parallel }^{2}-{k}_{\perp }^{2}}{| {\bf{k}}{| }^{2}}{P}_{0}\\ & & -{\rm{i}}\displaystyle \frac{1}{15}{a}_{0}{B}_{0}{c}_{{\rm{F}}}{P}_{{\rm{r}}}\displaystyle \frac{{k}_{\parallel }^{2}-{k}_{\perp }^{2}}{{k}_{\parallel }}\left(1-| {\bf{k}}{| }^{2}/28\right)\\ & & +\displaystyle \frac{3\kappa {c}_{{\rm{I}}}{c}_{{\rm{F}}}}{| {\bf{k}}{| }^{2}}\left[{k}_{\parallel }^{2}(1+\kappa )+{k}_{\perp }^{2}(1-\kappa )\right]{P}_{0}^{2}\\ & & +\displaystyle \frac{{\rm{i}}\kappa {c}_{{\rm{F}}}}{2{k}_{\parallel }}\left[{k}_{\parallel }^{2}(1+\kappa )+{k}_{\perp }^{2}(1-\kappa )\right]\left(1-| {\bf{k}}{| }^{2}/28\right){P}_{0}.\end{array}\end{eqnarray} \tag{ C.8 }$$

Equation (C.4) still involves the perturbation $\delta \hat{q}$ of the Lagrange multiplier. Utilizing again the incompressibility condition for the polar order parameter field P, yielding ${\bf{k}}\cdot \delta \hat{{\bf{P}}}=0$, and equation (C.4), we find that the perturbation $\delta \hat{q}$ must satisfy

$$\begin{eqnarray}&&\delta \hat{q}=-{\rm{i}}\displaystyle \frac{{\bf{k}}}{| {\bf{k}}{| }^{2}}\cdot {\bf{M}}\cdot \delta \hat{{\bf{P}}}.\end{eqnarray} \tag{ C.9 }$$

Inserting equation (C.9) into (C.4) we can identify the projector ${\boldsymbol{\Pi }}({\bf{k}})={\bf{I}}-{\bf{kk}}/| {\bf{k}}{| }^{2}$. 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

$$\begin{eqnarray}&&\sigma ({\bf{k}})=\displaystyle \frac{{k}_{\perp }^{2}}{| {\bf{k}}{| }^{2}}{M}_{\parallel \parallel }+\displaystyle \frac{{k}_{\parallel }^{2}}{| {\bf{k}}{| }^{2}}{M}_{\perp \perp }-\displaystyle \frac{{k}_{\parallel }{k}_{\perp }}{| {\bf{k}}{| }^{2}}\left({M}_{\parallel \perp }+{M}_{\perp \parallel }\right).\end{eqnarray} \tag{ C.10 }$$

For the explicit form of σ(k), see equations (15) and (16) in the main text.

Four dimensionless coefficients remain in the reduced and rescaled model equation (17). These are related to the coefficients of the full model given by equations (1), (2), (4) and (9) via

$$\begin{eqnarray}&&\lambda =\sqrt{\displaystyle \frac{5}{3{c}_{{\rm{I}}}}}\sqrt{\displaystyle \frac{-2{{\rm{\Gamma }}}_{4}}{{{\rm{\Gamma }}}_{2}}}\displaystyle \frac{1}{{{\rm{\Gamma }}}_{2}}{\lambda }_{0},\end{eqnarray} \tag{ D.1 }$$

$$\begin{eqnarray}&&a=1-\displaystyle \frac{4{{\rm{\Gamma }}}_{4}}{{{\rm{\Gamma }}}_{2}^{2}}\alpha ,\end{eqnarray} \tag{ D.2 }$$

$$\begin{eqnarray}&&b=\displaystyle \frac{4{{\rm{\Gamma }}}_{4}}{{{\rm{\Gamma }}}_{2}^{2}}\displaystyle \frac{5}{3{c}_{{\rm{I}}}}\beta ,\end{eqnarray} \tag{ D.3 }$$

$$\begin{eqnarray}&&{\tilde{B}}_{0}=\displaystyle \frac{2}{3}\displaystyle \frac{4{{\rm{\Gamma }}}_{4}}{{{\rm{\Gamma }}}_{2}^{2}}\sqrt{\displaystyle \frac{3{c}_{{\rm{I}}}}{5}}{B}_{0}.\end{eqnarray} \tag{ D.4 }$$

Here we present technical details for the weakly nonlinear analysis for the case ${\tilde{B}}_{0}\gt 0$. The analysis for the case ${\tilde{B}}_{0}=0$ is analogous, but less complex due to the overall rotational symmetry. In fact, the only non-standard feature in the field-free case is the coupling between the two leading-order amplitudes, see equation (22).

As introduced in section 5.1.2, the parameter ε denotes the distance to the bifurcation occurring at ${\tilde{B}}_{0}={\tilde{B}}_{0}^{\parallel }$ (when starting from large fields, where the system is homogeneous). At ${\tilde{B}}_{0}={\tilde{B}}_{0}^{\parallel }$, the instability sets in, that is, perturbations with a parallel wavevector start to grow and we observe a transition to a stripe pattern. Using the definitions of the long time and space scales, $T={\varepsilon }^{2}t$ and ${\bf{X}}=\varepsilon ({\bf{x}}-{{\bf{v}}}_{{\rm{g}}}t)$, the derivatives in the dynamic equation (17) for the effective velocity ${\bf{v}}$ are replaced by

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}\quad & \to & \quad {\partial }_{t}-\varepsilon {v}_{{\rm{g}}}^{\parallel }{\partial }_{{X}_{\parallel }}-\varepsilon {v}_{{\rm{g}}}^{\perp }{\partial }_{{X}_{\perp }}+{\varepsilon }^{2}{\partial }_{T},\\ {\partial }_{{x}_{\parallel }}\quad & \to & \quad {\partial }_{{x}_{\parallel }}+\varepsilon {\partial }_{{X}_{\parallel }},\qquad \quad {\partial }_{{x}_{\perp }}\to \quad {\partial }_{{x}_{\perp }}+\varepsilon {\partial }_{{X}_{\perp }},\\ {{\rm{\nabla }}}^{2}\quad & \to & \quad {\partial }_{{x}_{\parallel }}^{2}+{\partial }_{{x}_{\perp }}^{2}+\varepsilon \left(2{\partial }_{{x}_{\parallel }}{\partial }_{{X}_{\parallel }}+2{\partial }_{{x}_{\perp }}{\partial }_{{X}_{\perp }}\right)+{\varepsilon }^{2}\left({\partial }_{{X}_{\parallel }}^{2}+{\partial }_{{X}_{\perp }}^{2}\right),\\ {{\rm{\nabla }}}^{4}\quad & \to & \quad {\partial }_{{x}_{\parallel }}^{4}+{\partial }_{{x}_{\perp }}^{4}+2{\partial }_{{x}_{\parallel }}^{2}{\partial }_{{x}_{\perp }}^{2}\\ & & \qquad +\varepsilon \left(4{\partial }_{{x}_{\parallel }}^{3}{\partial }_{{X}_{\parallel }}+4{\partial }_{{x}_{\perp }}^{3}{\partial }_{{X}_{\perp }}+4{\partial }_{{x}_{\parallel }}^{2}{\partial }_{{x}_{\perp }}{\partial }_{{X}_{\perp }}+4{\partial }_{{x}_{\perp }}^{2}{\partial }_{{x}_{\parallel }}{\partial }_{{X}_{\parallel }}\right)\\ & & \qquad +{\varepsilon }^{2}\left(6{\partial }_{{x}_{\parallel }}^{2}{\partial }_{{X}_{\parallel }}^{2}+6{\partial }_{{x}_{\perp }}^{2}{\partial }_{{X}_{\perp }}^{2}+2{\partial }_{{x}_{\parallel }}^{2}{\partial }_{{X}_{\perp }}^{2}+2{\partial }_{{x}_{\perp }}^{2}{\partial }_{{X}_{\parallel }}^{2}+...\right.\\ & & \left.\qquad ...+8{\partial }_{{x}_{\parallel }}{\partial }_{{x}_{\perp }}{\partial }_{{X}_{\parallel }}{\partial }_{{X}_{\perp }}\right)+{ \mathcal O }({\varepsilon }^{3}).\end{array}\end{eqnarray} \tag{ E.1 }$$

Note that for the case ${\tilde{B}}_{0}=0$ the group velocity vanishes, ${v}_{{\rm{g}}}^{\parallel }={v}_{{\rm{g}}}^{\perp }=0$. 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 ${\tilde{B}}_{0}^{\parallel }$, only perturbations with a parallel wavevector grow. Therefore, we expand the effective velocity field ${\bf{v}}=({v}_{\parallel },{v}_{\perp })$ 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,

$$\begin{eqnarray}&&q({\bf{x}},t,{\bf{X}},T)={q}_{0}+\displaystyle \sum _{n=1}^{\infty }\displaystyle \sum _{m=0}^{n}{\varepsilon }^{n}{q}_{n,m}({\bf{X}},T){{\rm{e}}}^{{\rm{i}}m({x}_{\parallel }-{c}_{\parallel }t)}+{\rm{c}}.{\rm{c}}.\end{eqnarray} \tag{ E.2 }$$

Finally, the traveling speed ${c}^{\parallel }$ is also expanded in orders of ε,

$$\begin{eqnarray}&&{c}^{\parallel }={c}_{0}^{\parallel }+\displaystyle \sum _{n=1}^{\infty }{\varepsilon }^{n}{c}_{n}^{\parallel },\end{eqnarray} \tag{ E.3 }$$

where the zeroth component is proportional to the homogeneous stationary solution V₀, i.e. ${c}_{0}^{\parallel }=\lambda {V}_{0}$ (see equation (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, ${ \mathcal O }({\varepsilon }^{0})$, we recover equation (18) which determines the homogeneous stationary solution V₀.
• In first order, ${ \mathcal O }({\varepsilon }^{1})$, we find that the contributions ${v}_{1,0}^{\parallel }$, ${v}_{1,1}^{\parallel }$, ${v}_{1,0}^{\perp }$, ${q}_{\mathrm{1,0}}$, ${q}_{\mathrm{1,1}}$ and ${c}_{1}^{\parallel }$ must vanish to satisfy simultaneously the dynamical equation (17) and the incompressibility condition ${\rm{\nabla }}\cdot {\bf{v}}=0$. For ${v}_{1,1}^{\perp }$ we obtain
  $$\begin{eqnarray}&&0=(a-{{bV}}_{0}^{2}-\displaystyle \frac{3}{2}{\tilde{B}}_{0}{V}_{0}){v}_{1,1}^{\perp }={\varepsilon }^{2}{v}_{1,1}^{\perp },\end{eqnarray} \tag{ E.4 }$$
  where we used the definition of ε according to equation (23). Thus, we find that ${v}_{1,1}^{\perp }$ actually contributes in ${ \mathcal O }({\varepsilon }^{3})$ to the amplitude equation obtained below.
• In second order, ${ \mathcal O }({\varepsilon }^{2})$, using the incompressibility constraint, we find the relation
  $$\begin{eqnarray}&&{v}_{2,1}^{\parallel }={\rm{i}}{\partial }_{{X}_{\perp }}{v}_{1,1}^{\perp }.\end{eqnarray} \tag{ E.5 }$$
  As a solvability condition for equation (17) we obtain for the zeroth mode (m = 0)
  $$\begin{eqnarray}&&{v}_{2,0}^{\parallel }={D}_{0}| {v}_{1,1}^{\perp }{| }^{2},\end{eqnarray} \tag{ E.6 }$$
  where
  $$\begin{eqnarray}&&{D}_{0}=\displaystyle \frac{-2{{bV}}_{0}+{\tilde{B}}_{0}}{1+2{{bV}}_{0}^{2}+\tfrac{1}{2}{\tilde{B}}_{0}{V}_{0}}.\end{eqnarray} \tag{ E.7 }$$
  Matching all terms containing the first mode (m = 1) yields
  $$\begin{eqnarray}&&{q}_{\mathrm{2,1}}={\rm{i}}{D}_{\perp }{v}_{2,1}^{\parallel }=-{D}_{\perp }{\partial }_{{X}_{\perp }}{v}_{1,1}^{\perp },\end{eqnarray} \tag{ E.8 }$$
  where we inserted equation (E.5). The coefficient ${D}_{\perp }$ is given in equation (27) in the main text. Equations (E.5)–(E.8) show that the second order amplitudes ${v}_{2,0}^{\parallel }$, ${v}_{2,1}^{\parallel }$ and ${q}_{\mathrm{2,1}}$ are 'slaved' to the first order amplitude ${v}_{1,1}^{\perp }$. Further, we obtain the second order contribution to the traveling speed
  $$\begin{eqnarray}&&{c}_{2}^{\parallel }=\lambda {D}_{0}| {v}_{1,1}^{\perp }{| }^{2}.\end{eqnarray} \tag{ E.9 }$$
  All other higher order modes ($n\gt 1$), i.e. ${v}_{2,2}^{\parallel }$, ${v}_{2,0}^{\perp }$, ${v}_{2,1}^{\perp }$, ${v}_{2,2}^{\perp }$, ${q}_{\mathrm{2,0}}$, ${q}_{\mathrm{2,2}}$, ... either violate the incompressibility constraint or the solvability conditions obtained from equation (17), or they contribute in orders higher than ${ \mathcal O }({\varepsilon }^{3})$.
• In third order, ${ \mathcal O }({\varepsilon }^{3})$, we finally obtain a dynamic equation for the amplitude ${v}_{1,1}^{\perp }$ by matching terms containing the first mode (m = 1),
  $$\begin{eqnarray}\begin{array}{rcl}{\partial }_{T}{v}_{1,1}^{\perp } & = & +{v}_{1,1}^{\perp }-3b| {v}_{1,1}^{\perp }{| }^{2}{v}_{1,1}^{\perp }-2{{bV}}_{0}{v}_{2,0}^{\parallel }{v}_{1,1}^{\perp }-\displaystyle \frac{3}{2}{B}_{0}{v}_{2,0}^{\parallel }{v}_{1,1}^{\perp }\\ & & +4{\partial }_{{X}_{\parallel }}^{2}{v}_{1,1}^{\perp }-{\partial }_{{X}_{\perp }}{q}_{2,1}^{\parallel }.\end{array}\end{eqnarray} \tag{ E.10 }$$
  Inserting equations (E.6) and (E.8) into (E.10) and replacing ${v}_{1,1}^{\perp }\to A$ yields
  $$\begin{eqnarray}&&{\partial }_{T}A=A-g| A{| }^{2}A+{D}_{\parallel }{\partial }_{{X}_{\parallel }}^{2}A+{D}_{\perp }{\partial }_{{X}_{\perp }}^{2}A.\end{eqnarray} \tag{ E.11 }$$
  The diffusion coefficients ${D}_{\parallel }$ and ${D}_{\perp }$ are given in equation (27) in the main text. Further, the coefficient of the cubic term follows as
  $$\begin{eqnarray}&&g=3b+2{{bV}}_{0}{D}_{0}+\displaystyle \frac{3}{2}{\tilde{B}}_{0}{D}_{0}.\end{eqnarray} \tag{ E.12 }$$

Finally, scaling back to the fast time and length scales, t and ${\bf{x}}$, yields the amplitude equation (25) in the main text.

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).

Zoom In Zoom Out Reset image size

The emerging pattern near the bifurcation at ${\tilde{B}}_{0}^{\parallel }$ 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=\sqrt{{\sigma }_{\mathrm{Re}}^{\parallel }/g}$ of the amplitude equation (25) (and not to the homogeneous stationary solution V₀). The full stripe pattern solution on the level of the collective microswimmer velocity ${\bf{v}}$ is given as

$$\begin{eqnarray}\begin{array}{rcl}{v}_{\parallel } & = & {V}_{0}+{D}_{0}| A{| }^{2},\\ {v}_{\perp } & = & A{{\rm{e}}}^{{\rm{i}}({x}_{\parallel }-{c}_{\parallel }t)}+{A}^{* }{{\rm{e}}}^{-{\rm{i}}({x}_{\parallel }-{c}_{\parallel }t)},\end{array}\end{eqnarray} \tag{ G.1 }$$

where A^* denotes the complex conjugate of A. The form of the perturbation we consider is a perpendicular mode with critical wavenumber ${k}_{{\rm{c}}}=1$, i.e.

$$\begin{eqnarray}\begin{array}{rcl}{v}_{\parallel } & = & {V}_{0}+{D}_{0}| A{| }^{2}+\delta {\hat{v}}_{\parallel }{{\rm{e}}}^{{\tilde{\sigma }}_{\perp }t+{\rm{i}}{x}_{\perp }},\\ {v}_{\perp } & = & A{{\rm{e}}}^{{\rm{i}}({x}_{\parallel }-{c}_{\parallel }t)}+{A}^{* }{{\rm{e}}}^{-{\rm{i}}({x}_{\parallel }-{c}_{\parallel }t)}+\delta {\hat{v}}_{\perp }{{\rm{e}}}^{{\tilde{\sigma }}_{\perp }t+{\rm{i}}{x}_{\perp }}.\end{array}\end{eqnarray} \tag{ G.2 }$$

Utilizing the incompressibility condition ${\rm{\nabla }}\cdot {\bf{v}}=0$, we find $\delta {\hat{v}}_{\perp }=0$. We insert equation (G.2) into the dynamic equation for ${\bf{v}}$ (equation (17)) and obtain, after linearization, the growth rate of perpendicular modes with critical wavenumber k_c = 1,

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{\widetilde{\sigma }}_{\perp } & = & a-b(3{V}_{0}^{2}+3{D}_{0}^{2}| A{| }^{4}+6{P}_{0}{D}_{0}| A{| }^{2}+2| {A}^{2}| )\\ & & -2{\widetilde{B}}_{0}({V}_{0}+{D}_{0}| A{| }^{2}).\end{array}\end{array}\end{eqnarray} \tag{ G.3 }$$

The critical field strength ${\tilde{B}}_{0}^{* }$ defining the transition from the stripe pattern to the vortex lattice is then obtained by setting ${\tilde{\sigma }}_{\perp }=0$. Note that ${\tilde{\sigma }}_{\perp }$ depends on the amplitude of the parallel modes, A. This coupling is the reason why the field strength at the transition, ${\tilde{B}}_{0}^{* }$, is smaller than ${\tilde{B}}_{0}^{\perp }$.