Emergent pattern formation of active magnetic suspensions in an external field

We study collective self-organization of weakly magnetic active suspensions in a uniform external field by analyzing a mesoscopic continuum model that we have recently developed. Our model is based on a Smoluchowski equation for a particle probability density function in an alignment field coupled to a mean-field description of the flow arising from the activity and the alignment torque. Performing linear stability analysis of the Smoluchowski equation and the resulting orientational moment equations combined with non-linear 3D simulations, we provide a comprehensive picture of instability patterns as a function of strengths of activity and magnetic field. For sufficiently high activity and moderate magnetic field strengths, the competition between the activity-induced flow and external magnetic torque renders a homogeneous polar steady state unstable. As a result, four distinct dynamical patterns of collective motion emerge. The instability patterns for pushers include traveling sheets governed by bend-twist instabilities and dynamical aggregates. For pullers, finite-sized and system spanning pillar-like concentrated regions predominated by splay deformations emerge which migrate in the field direction. Notably, at very strong magnetic fields, we observe a reentrant hydrodynamic stability of the polar steady state.


I. INTRODUCTION
Self-propelled systems such as birds, fire ants and bacteria exhibit fascinating patterns of collective motion.Unraveling the physical principles governing collective self-organization of such autonomous systems have attracted tremendous attention in recent years.The efforts to understand the collective effects in self-propelled systems have led to emergence of the interdisciplinary field of active matter, see for example [1][2][3].Active matter is a fundamentally non-equilibrium class of materials which consist of particles transforming the ambient energy to some form of mechanical motion at the individual level.Many studies have focused on elucidating the influence of interparticle interactions on the collective behavior of active systems.It is found that the interplay between self-propulsion alone with simple short-ranged interactions in minimal models such as active Brownian particles with steric interactions [4][5][6] or Vicsek model with alignment interactions [1,7] leads to a rich phase behavior.Novel patterns of collective dynamics like dynamical clusters and traveling stripes have been identified [1,6] which have no counterparts in equilibrium systems.
The control of collective dynamics of microswimmers via an external field offers a promising route for high-tech applications such as micro-scale cargo transport, targeted drug delivery, and microfluidic devices [22][23][24][25].For instance, external magnetic field has been employed to control the rheological properties of magnetic swimmers [26,27].The collective dynamics of microswimmers in an external field is nonetheless poorly understood.Specifically, the effect of interplay between long-range hydrodynamic interactions and external fields on the pattern formation, with exception of few cases [27][28][29], has been little explored.
To make further progress in this direction, we focus on the large-scale collective dynamics of weakly magnetic microswimmers in a uniform magnetic field.We employ the kinetic theory framework [30] that allows us to overcome the size limitations of costly particle-based simulations and to capture the large-scale patterns of active suspensions over length scales much larger than the particle size.We provide an in-depth analysis of a kinetic continuum model that we have recently developed for dilute suspension of spherical microswimmers in an alignment field [31].Although our focus is on weakly magnetic swimmers, the model is in principle also applicable to bottom-heavy microswimmers in a gravitational field.
The kinetic model couples the Smoluchowski equation for arXiv:2006.03352v1[cond-mat.soft]5 Jun 2020 probability density function of fairly dilute active spherical suspensions in an alignment field to mean far field hydrodynamic interactions mainly generated by the swimmers motion.The hydrodynamic interactions are incorporated using the leading order flow field of a force-free microswimmer that is described by a force dipole.It decays as 1/r 2 where r is the distance from the swimmer.Independent of the details of motility mechanism, e.g.flagellar propulsion or surface distortions, the majority of microswimmers can be divided according to the their far field flow into pusher (extensile) and puller (contractile) swimmers, respectively.A pusher swimmer uses its tail to push fluid outward along its swimming axis whereas a puller swimmer employs its front appendages to pull the fluid towards its body in the direction of swimming.These two types of swimmers produce qualitatively different hydrodynamic flows and hence are expected to produce distinct spatio-temporal patterns.
We study the dynamics of both puller and pusher swimmers in a magnetic field by combining linear stability analysis and full numerical solution of 3D non-linear kinetic equations.Our linear stability analysis consists of investigating the stability of the probability density function of polar steady state as well as that of its orientational moments described by uniform density and polarization fields.Combining the two approaches we obtain complementary insights into the nature of instabilities.At low magnetic fields, a homogeneous weakly polarized state is stable, akin to an isotropic suspension of spherical swimmers.However, for sufficiently high activity strengths and moderately strong magnetic fields, a homogeneous polar state becomes unstable for both pushers and pullers.As we vary magnetic field and activity strengths, distinct spatio-temporal patterns emerge.At moderate field and activity strengths, pushers are driven by bend-twist hydrodynamic instabilities and form traveling bands perpendicular to the magnetic field.At stronger activity and field strengths, the density-driven hydrodynamic instabilities predominate pusher suspensions leading to formation of dynamical aggregates.Pullers at moderate field and activity strengths form system spanning pillars parallel to the field which are predominated by splay deformations.However, at stronger field and activity strengths, they form finite-sized pillar-like concentrated regions.Interestingly for very strong magnetic fields a homogenous polar state becomes stable again.Hence, we observe a re-entrant hydrodynamic stability; a hallmark of competition between alignment and hydrodynamic torques.
The remainder of this article is organized as follows.In section II, we discuss the main ingredients of the kinetic model for a dilute suspension of polar active particles in an alignment field.In section III, we analyze the linear stability of homogenous polar steady state to plane-wave perturbations for active polar suspensions aligned by an external field using a spectral method.Then, we calculate the stability diagram as a function of strengths of activity and magnetic field.In section IV, we first derive equations of motion for the orientational moments, density, polarization and nematic fields, using suitable closure approximations.Then, we analyze the linear stability of moment equations.In section V, we focus on numerical solution of the Smoluchowski equation coupled to the Stokes flow to explore the non-linear dynamics.We first outline our simulation method based on stochastic sampling method.Next, we investigate the emergent spatio-temporal pattern formation varying strengths of activity and magnetic field.We particularly discuss the distinguishing features of patterns observed at different instability regimes.Finally, our main conclusions and a discussion on comparison of linear stability analysis and non-linear dynamics solution can be found in section VI.
We consider a dilute suspension of N spherical magnetic microswimmers with a hydrodynamic radius a immersed in a fluid of volume V at a number density ρ m = N V .We assume that the self-propulsion is generated by a force-free mechanism of hydrodynamic origin such that its far field flow, averaged over swimmer's beat cycle, is well represented by that of a point-force dipole with an effective dipolar strength S eff [32][33][34].S eff depends on the geometrical parameters of the model swimmer [35][36][37][38], for instance on the body size a and the flagellum length [38].The translational and rotational friction coefficients of the swimmer are given by ξ r and ξ t .Each swimmer carries a weak magnetic dipole moment µ = µn along its body axis specified by the unit orientation vector n ≡ n and has a self-propulsion velocity U 0 n as depicted schematically in Fig. 1.The suspension is exposed to a uniform magnetic field B that exerts an alignment torque on each swimmer.We assume that µ is sufficiently small such that the dipole-dipole magnetic interactions at average inter-particle distance d int 3a are negligible relative to the thermal energy scale and no instabilities occur due to magnetic interactions.Therefore, for volume fractions Φ m 0.15 the dynamics of the system is governed by the interplay between the hydrodynamic interactions and the field-induced alignment torque.

B. Conservation equation: the Smoluchowski equation
For sufficiently low ρ m , the mean-field configuration of an ensemble of swimmers at a time T can be described by a single-particle distribution function Ψ(X, n, T ), i.e., the degrees of freedom of other particles have been traced out by integration.The function is normalized as As such Ψ(X, n, T )/N describes the probability density of finding a particle with the center of mass position X and the orientation vector n at time T .Therefore, a uniform and isotropic state can be described by the constant distribution function Ψ = ρ m /4π.The kinetic model for hydrodynamically interacting swimmers in an external field [31] is based on an evolution equation for the distribution function Ψ(X, n, T ) coupled to an equation for the mean-field fluid velocity U.The Smoluchowski equation for hydrodynamically interacting active particles carrying a weak magnetic dipole moment in an external field is given by in which ∇ • n ≡ (1 − nn) • ∇ n with dyadic product defined as (nn) i j = n i n j denotes the angular gradient operator and v x and v n are the translational and rotational flux velocities resulting from a swimmer's drift.D = D t ∇ 2 + D r ∇ • n 2 is the diffusion operator in which D t and D r describe the effective long-time translational D t and rotational diffusion coefficients, respectively.The diffusion coefficients can result from thermal or biological fluctuations, e.g., due to tumbling of bacteria in the case of rotational diffusion.The translational flux velocity includes the drift contributions from the self-propulsion U 0 n and an advection due to the local flow field U.The rotational flux velocity v n ≡ ṅ is modeled as in which P ⊥ n = 1 − nn, describes the projection operator to the space orthogonal to the orientation vector and W = 1 2 (∇U − (∇U) ), with (∇U) i j = ∂ i U j , is the vorticity tensor.The flux velocity v n includes the rotational drift contributions resulting from the torque due to the magnetic field and vorticity of the local flow.Using the relation between the angular velocity ω and rate of change of orientaitona vector ṅ = ω ×n, the first term P ⊥ n • µ ξ r B on the right-hand side is obtained from the balance between the magnetic torque τ = µn × B and the frictional hydrodynamic torque −ξ r ω in the overdamped and low Reynold's number limits.The second term models the interaction of a spherical swimmer with the local flow vorticity based on the the second Faxen's law [39].
From the distribution function, we define the local density field ρ(X, T ), polarization field p(X, T ), and the nematic order parameter field Q(X, T ), as the symmetric and traceless parts of the zeroth, first, and second order orientational moments of Ψ(X, n, T ) with respect to n, respectively, These moments will be used throughout the paper in the following sections.

C. Mean-field flow
The flow field U in Eqs. ( 3) and ( 4) may result from an imposed external flow or from hydrodynamic interactions.In this work, we consider the case that there is no external flow and U solely represents the self-generated flow due to motion of swimmers.In the limit of vanishing Reynolds number, applicable to microswimmers, the fluid reacts in good approximation instantaneously to changes in the particle configuration.The mean-field flow U[Ψ] resulting from the hydrodynamic interactions between the swimmers is well captured by the incompressible Stokes equation in which P and η denote the isotropic pressure and the viscosity of the suspending fluid and ∇ • Σ = ∂ X i Σ i j ê j .The meanfield stress Σ[Ψ] depends on the instantaneous suspension configuration encoded by Ψ.In the case of microswimmers, it can be decomposed into the sum of several contributions, arising from the self-propulsion, Brownian rotations, resistance to stretching and compression by the local flow field and steric and magnetic torques.For dilute suspensions of spherical microswimmers, we neglect stresses arising from Brownian rotations (can be incorporated into active stress by modifying the prefactor), inextensibility of the particles and steric torques because of their small contributions.We only consider active stress Σ a [Ψ], generated by the self-propulsion of swimmers, [33,34] and a magnetic stress Σ m [Ψ], caused by reorientation of swimmers in the external field.Hence, the stress in our model is given by Σ In a dilute suspension, for which the average the ratio of inter-particle distance to the swimmer size is large, the active stress of a force-free microswimmer Σ a [Ψ] can be modeled as that of a point-force dipole -the leading order non-zero singularity of the Stokes flow [34,40,41].The active stress of a suspension of dipolar microswimmers is proportional to the nematic order tensor field [30,42] as defined by Eq. ( 7): It can be interpreted as a superposition of stress contributions of all possible swimmer orientations at position X.The strength of the active stress is determined by the amplitude Σ a = −ρ m S eff .The sign of Σ a determines the nature of the swimmers, being a puller Σ a > 0 or a pusher Σ a < 0. The torque due to external field M B = µn × B leads to rotation of swimmers that in turn exerts a rotational stress on the fluid while dragging the surrounding fluid layers.This results in an antisymmetric stress contribution of the form in which B = B/B, Σ m ≡ ρ m µB, and ε is the Levi-Cevita symbol.This stress contribution is identical to that of passive magnetic suspensions.Note that the symmetric part of the magnetic stress is zero for spherical particles [43,44].

D. Non-dimensionalization
To facilitate the analysis of the model, we render the equations dimensionless, using the following characteristic velocity, length, and time scales: u c = U 0 , t c = 1/D r and x c = U 0 /D r .Note that our choice of characteristic time and length scales are different from our previous work [31].We rescale distribution function with the number density such that ψ(x, n,t) ≡ Ψ(xx c , n,tt c )/ρ m , is dimensionless and ψ/v represents a probability density normalized to unity: where v = V /x 3 c .The form of Smoluchowski equation for ψ(x, n,t) remains unchanged where the gradient operator ∇ ≡ ∂ /∂ x i êi is now with respect to the reduced coordinates.The dimensionless spatial and rotational flux-velocities reduce to in which α e = µB ξ r D r defines the alignment parameter.Likewise, the dimensionless diffusion operator simplifies to where d t = D t D r /U 2 0 is the reduced translational diffusion coefficient.The equation for the flow-field transforms into with the dimensionless stress tensor σ given by As such, two additional independent dimensionless parameters, the active stress amplitude σ a = Σ a D r η and the external field-induced stress amplitude σ m = Σ m D r η appear in our model.

III. LINEAR STABILITY ANALYSIS OF HOMOGENEOUS POLAR STEADY STATE
The set of equations ( 13) and ( 17) forms a closed system that can be solved for the evolution of the distribution function ψ and the flow field u in the suspension.However, it is not presently feasible to solve these coupled equations analytically.Therefore, we resort to the linear stability analysis that provides us with some degree of predictive insight into the dynamics of the equations with respect to a suitable base state.This kind of analysis allows us to divide the parameter space into a stable region described by the base state and an unstable region with yet unknown dynamics departing from the base state.Furthermore, the linear stability analysis offers some valuable insight into the dynamics at the onset of instability.

A. Homogeneous and steady solution as a base state
The external field breaks the rotational symmetry of the system but preserves translational invariance.Thus, we first seek for spatially-uniform ∂ x ψ 0 = 0 and steady ∂ t ψ 0 = 0 solutions of the Smoluchowski equation (13).Solutions of the form ψ 0 (α e , n) will serve as base states for the linear stability analysis.For a homogenous steady state, all spatial and time derivatives in Eq. ( 13) vanish.The same holds for the flow field as ∇ • σ = 0 ⇒ u 0 = 0.Only the rotational flux velocity terms remain.Hence, ψ 0 (α e , n) can be obtained by setting the total rotational flux velocity including both drift and diffusive contributions to zero, i. e., Solving this equation yields with the normalization S 2 dn ψ 0 (α e , n) = 1.This steady-state is identical to that of passive magnetic dipoles in an external field [43].For passive systems at the thermal equilibrium, the Einstein-Stokes-Debye relation D r = k B T ξ r [39] holds and the alignment parameter becomes α e = µB k B T which is equal to ratio between magnetic and thermal energy scales.More generally, it describes the ratio between two characteristic reorientation time scales α e = τ e /τ r .The time τ r = 1 D r represents the average decorrelation time of a diffusive particle from its initial orientation and τ e = ξ r mB is a measure of the typical alignment time of a non-diffusive dipole with the external field.The competition between alignment (order) and the randomization of orientation (disorder) determines the degree of alignment quantified by the mean polarization p 0 .It is the given by the magnitude of the polarization vector p 0 : The function p 0 (α e ) is identical to the well-known Langevin function appearing in the context of paramagnetism or forceextension relation of a freely jointed chain [45].
Assuming a magnetic field parallel to the z-axis, i.e.B = Bẑ, without loss of generality, the homogeneous polar state with axial symmetry takes the simple form of ψ 0 (α e , n) = ψ 0 (α e , θ ), where θ denotes the angle between the orientation vector and the magnetic field and it coincides with the polar angle in spherical coordinates for the orientation n(θ , φ ) = (sin θ cos φ , sin θ sin φ , cos φ ).The angular dependency of the homogeneous polar steady state for different values of α e is shown in Fig. 2 (a).A strong external magnetic field (large α e ) results in a focused angular distribution around the magnetic field axis corresponding to θ = 0 and thus a large mean polarization p 0 .The functional dependency of polarization magnitude on the alignment parameter α e is plotted in Fig. 2 (b).
The mean polarization continuously increases with increasing α e ∝ B. It asymptotically approaches a perfectly aligned state with p 0 = 1 in the limit of very large α e described by lim α e →∞ ψ 0 = δ (n − B).In the other extreme of very low magnetic field strengths, fluctuations will increasingly decorrelate the orientation of a swimmer, leading to a flat profile in the angular distribution, i. e., lim α e →0 ψ 0 = 1/4π, which corresponds to an isotropic suspension with p 0 = 0. We now proceed to analyze the linear stability of the homogeneous polar steady state presented in Sec.III A. We consider a small disturbance of the distribution function ψ with respect to ψ 0 .
where |ε| 1 and |ψ p | ∼ O(1).Likewise, the flow-field of the suspending medium is perturbed by, u = u 0 + εu p , in which u p is the flow-field caused by the perurbation ψ p .The corresponding flow field of the steady state is u 0 = 0, because all the spatial derivatives on the right-hand side of Stokes equation ( 17) vanish for ψ 0 (n).
After neglecting terms of O(ε 2 ) in the governing equations, we obtain the following linearized evolution equation for ψ p where the double contraction . . is defined as (ab . .C) = a i b j C i j .In our derivation, we have used antisymmetric property and tracelessness of the vorticity tensor W and the following identities: which hold for any arbitrary tensor A and vector a.The flow field resulting from the perturbation u p satisfies the same momentum equation as u, but forced by the linearized stress tensor given by where the time-dependence of the stress tensor stems from that of ψ p (x, n,t).
To progress further, we Fourier-transform the linearized Smoluchowski equation, Eq. ( 22) where the Fourier transform of ψ p is defined as ψ F p = dxψ p e ik•x and we use the factorization ansatz ψ F p (k, n,t) = ψ(k, n)e λ (k)t .This ansatz decomposes the contribution of Fourier mode of the perturbation ψ F p to a time-independent amplitude ψ(k, n), also known as mode shape, and an exponential growth factor with a complex growth rate given by λ (k).This ansatz, which arises from linearity of perturbation equations, implies that the Fourier transform of the stress tensor due to perturbation given by Eq. ( 23) can be written as σ F p = σp (k)e λ (k)t , where Consequently, the Fourier transform of the flow field u p can be obtained as is the Fourier transform of the Oseen tensor and k = k −1 k is the normalized wavevector.From above, we can see that the Fourier transform of the flow field can also be decomposed as u F p (k) = ũ(k)e λ (k)t where its amplitude is explicitly given by ũ After some algebraic manipulation, the governing equation for ψ(k, n) transforms into an eigenvalue problem of the form in which L represents a linear differentio-integro-operator and ψ(k, n) is the associated eigenvector encoding the form of the orientational perturbation for a given k.The explicit form of L is given by Based on the form of Eq. ( 28), we note that the stability is governed by four dimensionless parameters d t , α e , σ a and σ e .To determine the growth rate λ (k) and hence stability of the active suspension in external field for a given set of the parameters, we need to solve the eigenvalue problem defined by equation (27).We discuss our methodology for this problem in the following subsection.

Spectral method for solving the eigenvalue problem
The above analysis shows that it is sufficient to consider plane wave perturbations of the form: to investigate the linear stability of the steady state.Here, Re λ determines growth rate and Im λ gives the frequency of a travelling wave with wavevector k.To solve the eigenvalue problem of ( 27), we employ a spectral method where we also expand the orientational dependency of the eigenfunction ψ(k, n) as well as ψ 0 in the basis of spherical harmonics.We choose a spherical coordinate system in which B is aligned with the polar axis.Denoting the polar and azimuthal angles by θ ∈ [0, π] and φ ∈ [0, 2π), respectively, we have n = (sin θ cos φ , sin θ sin φ , cos θ ) In this coordinate system, the spherical harmonic function of degree l and order m = −l, . . ., l is defined as where P m l (cos θ ) is the associated Legendre polynomial.The spherical harmonics satisfy the orthogonality condition: where the scalar product is defined by with the star operator • * representing the complex conjugatation.These functions form a complete basis on the unit sphere, on which we expand the mode shape ψ(k, n) as where is the coefficient corresponding to spherical harmonics Y m l .After substituting Eq. ( 33) into Eq.( 27) and applying the orthogonality condition Eq. ( 31), the eigenvalue problem for the mode shape ψ(k, n) reduces into an algebraic eigenvalue problem for the vector | ψ whose components are given by the harmonic amplitudes ψ h l : . Expanding the operator L defined by Eq. ( 28) on the spherical harmonics basis generates terms which are products of two spherical harmonics.The product can in general be written as the following linear combination of spherical harmonics included in common software packages and computer algebra applications such as Mathematica.The tensor L in Eq. ( 34) is of infinite size, hampering further analytical progress.We solve the algebraic eigenvalue problem by truncating the sum at sufficiently large j = j max such that the convergence of the dominant eigenvalues and eigenvectors are ensured.The number of angular modes that have to be included for convergence depends on α e partly because of the growing number of modes needed to accurately represent the steady state for large α e .Truncating the coefficient tensor L introduces an error in the calculation of the eigensystem.However, the error gets progressively smaller and has rapid convergence when adding further modes.The value of the largest growth rate Re λ max and as a function of j max is plotted in Fig. 3 for two different values of α e .For α e = 4, we find that j max = 5 ( 55 angular modes) is sufficient to obtain a good convergence whereas for α e = 10, at least j max = 10 ( 210 angular modes) is required for a reasonable convergence.

C. Linear stability of homogeneous polar steady state
As discussed earlier, the linear stability of the homogeneous steady state ψ 0 in equation ( 20) depends on four dimensionless parameters d t , α e ∝ µB, σ a and σ e .Additionally, the eigenvalue problem defined by Eqs. ( 27) and (28) and thus the stability of the steady state depends on the direction of the wavevector k with respect to the field direction as the external field breaks the rotational symmetry.However, the system still holds an axial symmetry around the B axis.Hence, the direction of wavevector can be characterized by a single angle between the magnetic field and the wavevector Θ B = cos −1 ( k • B).For a given solvent viscosity and density of active particles, the experimentally tuneable parameters are  20) as a function of the dimensionless active stress σ a and alignment parameter α e ∝ B; while setting σ m = 0.4 α e and d t = 3 × 10 −6 .The borderline of neutral stability (red dash-dotted line) is calculated by finding Re λ max (k) = 0.The dashed amber lines correspond to the cases where Re λ max (k || ) = 0 for pushers and Re λ max (k ⊥ ) = 0 for pullers.The solid blue lines represent Re λ max (k) = 0 based on the linear stability analysis of density and polarization fields from truncated moment equations for wavevectors parallel and perpendicular to the magnetic field.On the right side, the corresponding polarization p 0 of the steady state ψ 0 (α e ) is plotted.For the points marked by crosses, the behavior of growth rate and pattern formation are further discussed in the paper.
the strengths of activity and magnetic field.Therefore, we construct a stability diagram as a function of α e ∝ µB and σ a .We set d t = 3 × 10 −6 ; chosen to be comparable to the parameter ranges relevant for magnetotactic bacteria [20] and vary the magnetic stress concomitantly with α e as σ m = 0.4 α e .
For a given set of parameters, ψ 0 is unstable if the maximum growth rate is positive for at least one mode parametrized by (k, Θ B ). Based on the results of linear stability analysis, we divide the parameter space spanned by (σ a , α) into stable and unstable regimes with respect to the steady state ψ 0 [31].In the stable regime, the system evolves towards the steady state ψ 0 and becomes stationary.In the unstable regime, even small fluctuations make the system depart from ψ 0 towards a non-trivial dynamics.A line of neutral stability, i.e., Re λ max = 0 divides the two regimes.In stability diagram of Fig. 4, the red dashed-dotted lines represent the lines of neutral stability for pushers σ a < 0 and pullers σ a > 0. On the right panel, the mean polarization p 0 (α e ) given by Eq. ( 21) is plotted, highlighting the dependency of the steady state ψ 0 on α e .The steady state ψ 0 is stable for either of small activity |σ a | 20 or a low external magnetic field α e 0.5.In the case of small σ a , hydrodynamic interactions are too weak to destabilize the steady state.For a small α e , the polarization p 0 (α e 0.5) 0.3 is rather weak and our system akin to an isotropic suspension of spherical swimmers remains stable.For a sufficiently large active stress σ a 20 and a moderate external magnetic field strength, the homogeneous polar steady state becomes unstable.In this regime, the combined effect of sufficiently strong hydrodynamic interactions ∝ σ a and orientation fluctuations drive the system away from a uniformly aligned state.Interestingly, by further in- creasing the external magnetic field strength, the steady state becomes stable again, and we observe a reentrant hydrodynamic stability.Reentrant stability at strong external fields is a consequence of magnetic torque overcoming the hydrodynamic torque.For given active stress amplitude σ a , the active force f a = σ a ∇ • Q and its resulting hydrodynamic torque has an upper bound that can be overcome by the alignment torque for sufficiently strong magnetic fields.As a result, the steady state becomes stable again.
Analyzing the nature of instability in the unstable regions, we recognize four distinct types of instability.To demonstrate their distinct nature, four representative points, corresponding to each type are picked out from the stability diagram.These points are marked by crosses in Fig. 4 and correspond to (σ a , α e ) = (−30, 4), (30,4), (−40, 20) and (40,20) in the parameter space.Fig. 5 shows the real and imaginary parts of the complex growth-rate with the largest real part, the socalled maximum growth rate (Re λ max ) and its oscillation frequency Im λ max , as a function of k = |k| at various perturbation angles Θ B for each of the points.We note that the maximum growth rate strongly depends on the direction of the perturbation wavevector, i. e., Θ B .For puller and pusher swimmers with equal strength of activity |σ a | = 30 and magnetic field α e = 4, long-wavelength perturbations k → 0 dominate the instabilities and destabilize the homogeneous polar state ψ 0 .However, for pushers σ a = −30, fluctuations in the direction of magnetic field grow fastest, whereas for pullers σ a = 30 both perturbation directions parallel and perpendicular to B predominate the system.Therefore, we expect dis-tinct patterns of instabilities for pushers and pullers as confirmed by our non-linear simulations presented in Section V. Notably, in both cases Im λ max (k → 0) → 0 which implies that the large wavelength fluctuations grow monotonically with time.Interestingly, for pushers and pullers with stronger activity σ a = ∓40 and much larger magnetic field α e = 20 but still in the unstable regime, the wavenumber corresponding to the maximum growth-rate k max ≈ 2.5 is finite and it occurs at Θ B ≈ 30 • for pushers and Θ B ≈ 60 • for puller, featuring clearly different instability regimes.In the case (d), Im λ (k max ) of the wavevector with the maximum growth-rate is non-zero pointing to the oscillatory behavior of the predominant growth mode.These examples represent the four distinct types of instabilities observed: parallel and perpendicular orientational instabilities for pushers and pullers at moderate external field strengths, and more complex perturbation structures at higher field strength featuring a finite characteristic wavenumber for the largest growth rate at an activitydependent angle intermediate between parallel and perpendicular directions.

IV. LINEAR STABILITY OF ORIENTATIONAL MOMENT EQUATIONS
In this section, we take an alternative approach for investigating the linear stability of the homogeneous polar steady state.Instead of expanding the orientational part of the single particle distribution ψ in terms of spherical harmonics, as done in equation (32), it equivalently can be expanded in terms of dyadic products of the orientation vector n [48,49], where the l-fold dyadic product is denoted by and l denotes the l-fold contraction of two tensors, for tensors A ∈ T n 0 (R 3 ) and B ∈ T 0 m (R 3 ) using Einstein summation convention.As shown in reference [48], the coefficients M ( j, x, n) are proportional to the orientational expectation values of the symmetric and traceless (irreducible) part of n ⊗ j , where the orientational expectation value is defined as • n = S 2 dn ψ •.They are called the orientational moments of ψ.Especially, the zeroth, first and second moments coincide with the density, the polarization, and the nematic fields defined in equations (5-7).Hence, the distribution ψ expanded in terms of the orientational moments reads Truncating the moment expansion allows us to manipulate the resulting terms algebraically and to find approximate analytical expressions for the linear stability analysis.Evolution equations for each of the orientational moments can be directly derived by taking moments of the Smoluchowski equation (13).The dynamical moment equations for the first three moments are presented in the subsequent subsection.

A. Equations of moments
The time evolution of the density field ρ(x,t) = 1 n can be derived by tracing out the angular dependency in equation (13).Using the identity S 2 dn ∆ • n ψ = 0, its time evolution is given by where D t ≡ ∂ t + u • ∇ represents the material derivative.This represents a convection-diffusion type of equation with a source term which originates from the divergence of polarization field.Likewise, the time evolution of the polarization field p(x,t) = n n can be obtained by taking the first moment of equation ( 13) and using the following identities. and in which v t represents any tangential vector field on a sphere fulfilling v t • n = 0, and in our context, it is given by the rotational drift velocity v n defined by equation ( 15)).Consequently, the evolution equation for the polarization field is given by which is again a convection-diffusion type of equation with a more complex source term including contributions from density gradient, polarization and the divergence of the nematic tensor field and terms arising from the interaction of the active particles with the local flow vorticity and the magnetic field.Lastly, we obtain the time evolution of the nematic tensor field Q(x,t) = nn − 1/3 n by integrating equation ( 13) with S 2 dn nn • and using the following identities: The moment equation for Q is eventually given by in which we have used the antisymmetric property of the vorticity tensor W to further simplify the equation.The dynamics of nematic tensor is directly affected by source terms steming from the divergence of polarization.Similar to nematodynamics equation of active nematics [50], the dynamics of Q is strongly coupled to the flow velocity through the advection and vorticity terms.Moreover, additional terms appear due to coupling to the external field.
As can be seen from the equations of moments Eq. ( 38), Eq. ( 41), and Eq. ( 45), they constitute a hierarchy of equations where each moment equation depends on higher moments.In order to proceed further, we break this hierarchy by introducing the following closure relations which are compatible with a polar steady state All higher order moments are neglected.These closure approximation is sometimes referred to the Hand-closure [51].We will see in the next subsection that it generates perturbations that are structurally consistent with the results of linear stability analysis of steady distribution function ψ 0 in Section III B 1.

B. Stability of moments
After establishing the moment equations of the system, reduced to the density and polarization field, we proceed with their respective linear stability analysis employing the above closure relations.The homogeneous steady state solution of moment equations, Eqs.(38) and ( 41) is given by (ρ 0 , p 0 ) where which become equivalent to 1 ψ 0 n and n ψ 0 n = α e coth α e −1 α e (see Eq. ( 21)) for sufficiently small α e .Applying the closure relations Eq. ( 46) yields Q 0 = p 0 p 0 − 1 3 1.Next, we linearly disturb the steady state by small perturbations of the form in which 0 < ε 1 and u 0 = 0, and substituting them into equations (38) and (41) the dynamics of the linearized perturbations in O(ε) can be formulated as Analogous to Section III B 1, we make an eigenmode ansatz for the perturbations of the moments as below: Here, the flow field perturbation ũ is mainly driven by Q as in which k = k −1 k and Q ≡ p 0 p + pp 0 .Here, we have neglected the contribution of magnetic stress due to its small effect.Substituting the eigenmode ansatz for perturbations back into the density and polarization equations, they transform into where again W = i 2 ( ũk − k ũ).To analyze these equations, we consider perturbations of polarization which are perpendicular to the external field B = Bẑ and polarization p 0 ≡ p 0 ẑ.Without loss of generality, we assume the perturbation to be in the x-direction i.e., p ≡ p(k)x.This assumption reduces the problem to 2D and moreover, p • p 0 = 0.It also implies that small linear perturbations perpendicular to p 0 practically influence the orientation but not the magnitude of the polarization and allows us to investigate hydrodynamically induced orientational instabilities.Moreover, it naturally generates a nematic perturbation of the form pp 0 + p 0 p which is compatible with the closure approximation given by Eq. (46).Physically, a polarization perturbation of this form corresponds to bend (for pushers) and splay (for pullers) deformations of the polarization field [52], to be discussed in the following Section.Setting p • p 0 = 0 and thus focusing on ρ and p⊥ , the equations ( 60) and ( 61) can be written in the reduced form of with the operator 6.The growth rates Re λ against the alignment parameter α e of a linear perturbation of the polarization field, equation ( 64).This illustrates the re-entrant stability upon increase of magnetic field.The blue line is the eigenvalue of a mode crossing over from negative to positive and back to negative growth rate.It is constructed from the eigenvalues λ 1,2 defined in equation ( 63).
Here, k = k sin(Θ B )ê x + k cos(Θ B )ê z was used, where Θ B denotes the angle of the direction of the perturbation with respect to the external field axis pointing in the z-direction, Θ B = ∠(k, B).The eigenvalues can be found analytically [53].Solving this eigenvalue problem, we find that the largest eigenvalue can be found at k → 0, in agreement with the findings of Section III C for the cases (a) and (b) in Fig. 5. Therefore, for stability regimes (a) and (b), it is sufficient to restrict the analysis to parallel and perpendicular modes of truncated moment equations to assess the stability of a uniform polarization field.The k → 0 eigenvalues read with a = −4 − 2α e p 0 − p 2 0 σ a cos 2Θ B .The largest eigenvalue is given by λ max = max(0, a/2) (and the smallest by min(0, a/2)).However, simulations show that the perturbation associated with the λ = 0 can be considered stable.On the other hand, the eigenvalue corresponding to the mode that becomes unstable, i.e., the one having a change of sign in its growth rate, can be constructed by combining the non-zero parts of eigenvalues into one, yielding The dependence of the non-zero eigenvalue λ on the alignment parameter α e is plotted in Fig. 6.This result demonstrates that the simplified approach is sufficient to recover the re-entrant stability obtained earlier based on the full linear stability analysis of the steady state.
The line of neutral stability can be found by solving λ = 0 for σ a : It is shown in Fig. 4 in direct comparison to the result of linear stability analysis of the steady state distribution function ψ 0 .The line of neutral stability based on the stability of the first two moments (blue solid line) nearly coincides with the results of the parallel and perpendicular perturbations obtained from the linear stability analysis of ψ 0 discussed in Section III C (red dashed line).For larger external field strengths, the full analysis reveals additional unstable modes that do not fall into the same scheme.For the regions between the red dash-dotted lines and dashed amber lines in Fig. 4), a wider spectrum of orientational modes contribute to the instability that are not captured by the closure approximation of equation (46).Additional moments would be required to obtain a complete description.

V. NONLINEAR DYNAMICS SIMULATIONS
The linear stability analysis predicts the stability of a given steady state and provides us with a qualitative insight into the dynamics near it.However, as the system departs from the initial steady state, non-linearities prevail the dynamics and the linearized equations fail to describe the dynamics correctly.Therefore, it is necessary to investigate solution of the full non-linear equation (13).Below, we first outline our methodology for solving the full non-linear Smoluchowski equations coupled to the Stokes flow.Then, we discuss the pattern formation emerging from the long-time dynamics of active magnetic swimmers in the external field.

A. Numerical simulation method
Our methodology consists of a hybrid stochastic particle based sampling method for obtaining Ψ(x, n,t) in the Smoluchowski equation with periodic boundary conditions.It is based on integrating coupled translational and rotational Langevin equations, which are the counterpart of the Smoluchowski equation.They are given by: where v x and v n are defined by equations ( 3) and ( 4).Γ and Λ represent the stochastic force and torques with the following statistical properties: Within our theory, direct inter-particle dependencies are replaced by mean-field interactions.Consequently, once the mean-field stress profile and the resulting flow is computed from the distribution function, different initial value problems for a given particle can be simulated independently of each other.This realization is the basis of our numerical method that we dub "Stochastic Sampling" method.Fig. 7 summarizes the flow diagram of our method.A detailed description of the methodology can be found in reference [54].
To solve these equations numerically, we employ the Euler forward integration scheme based on the Itô interpretation of noise.For every time step, we integrate the corresponding Langevin stochastic differential equations for the FIG. 7. Flow diagram illustrating the algorithm of Stochastic Sampling method used to solve full non-linear Smoluchowski equation coupled to the mean-field Stokes flow.Given a probability density function (PDF) Ψ, the mean flow field can be calculated by evaluating the stress profile based on Eq. (18).Then, the PDF can be sampled to obtain particles which are integrated using the interaction field.Employing the updated sample configuration the PDF can be estimated by a Kernel Density Estimation (KDE).By iterating the process, the Smoluchowski equation describing the dynamics of the PDF can be integrated in time.
positions and orientations of a large number (10 6 − 10 8 ) of independent and randomly initialized test particles.These sample configurations provide us with sufficient statistics to infer a realization of the distribution function ψ(x, n,t) by using a kernel density estimation method.From the estimated distribution function, we compute the stress profile in the fluid.We use a spectral method based on the decomposition of u into the Fourier modes for solving the Stokes equation.In the Fourier domain, the flow field is obtained as p is the Fourier transform of the stress tensor.Given the stress profile, the flow field in terms of its Fourier modes on a periodic lattice is obtained.Then, it is Fourier transformed back to the real space.Eventually, u(x,t) is fed back into the next integration time step for the Langevin equations.
In the reported numerical simulations, we use a grid of 100 lattice points with box dimensions of 5 x c for each of the spatial coordinates, and 24 and 16 points for the spherical polar and azimuthal orientational coordinates θ and φ in n(θ , φ ).This choice of box dimension ensures that the initial perturbation spans both unstable and stable modes for all the four instability regimes presented in Fig. 5.The simulations, conducted in a box size of V sim = (5 x c ) 3 , are initialized with the homogeneous polar steady state ψ 0 , given by equation (20) for different system parameters.As in the previous sections, the translational diffusion is fixed to d t = 3 × 10 −6 , the alignment stress magnitude is set to σ m = 0.4α e , while the external magnetic field and the active stress magnitude are varied through the alignment parameter α e and the active stress amplitude σ a .

B. Pattern formation and nature of instabilities
Starting from a spatially homogeneous polar state ψ 0 (α e ) given by Eq. ( 20), we evolve ψ and u for each state point characterized by the (α e , σ a ) pair by employing the method outlined above.For the unstable state points, ψ departs from ψ 0 significantly whereas for the stable points the system converges towards ψ 0 even starting from an initial uniform isotropic state.The snapshots of Fig. 8 depict the time evolution of density field, projected to 2D by averaging along the y-axis, for (α e , σ a ) values corresponding to the points marked by crosses in Fig. 4. We observe a general trend that a uniform density profile becomes unstable towards density fluctuations.Over time, small-scale fluctuations disappear and the the field profiles become smoother owing to diffusion.Only predominant fluctuations at wavelengths of the order of the box size persist.As a consequence, smooth non-uniform density, polarization and flow fields develop.At long times, the configuration of the active suspensions is not steady but constantly fluctuates in time.The distribution of swimmer orientation appears to converge towards a dynamical steady state which depends on α e and σ a , leading to a constant average polarization in time [31].In contrast, the density fields exhibit distinct spatial patterns for different instability regimes which keep evolving and reorganizing.Using the results from the non-linear dynamics simulations, we assess the validity of the phase diagram of Fig. 4 predicted by the linear stability of ψ 0 (α e ).Fig. 9 presents an overview of density field projections into the x-z-plane, at a late time t=1800, after the instability has already established itself, for different values of the active stress σ a and the alignment parameter α e .To compare with the linear stability analysis predictions, we have plotted the lines of neutral stability (red lines), for which the largest growth rate is zero, i.e., Re λ max = 0. Additionally, the dashed amber lines depict the borderlines beyond which instability is governed by parallel and perpendicular perturbations for negative (pusher) and positive (puller) σ a , respectively.Comparing the predictions of the linear stability analysis with results of simulations for different activity and magnetic field strengths, find excellent agreement.ψ 0 (α e ) is stable for the (σ a , α e ) values where density profile remains homogenous, whereas ψ evolves towards an inhomogeneous time-dependent density profile for the regions that are unstable according to the linear stability analysis.
For moderate external magnetic field strengths ∝ α e and moderate activities ∝ σ a , corresponding to the unstable regions beyond the amber lines, patterns with the most distinct characteristics appear.The higher the external field, the patterns display finer structures (of higher spatial frequency) in the density profile suggesting a the predominance of the characteristic length scale associated with the external field varying inversely with magnetic field strength.Indeed, based on dimensional analysis, one can identify a length scale e ∝ D t η/ρ m µB ≡ D t η/Σ m .In these regions, the basic characteristics seem to be conserved, with a band-like structure for pushers and a pillar-like structure for pullers.In the intermediate regions between the red and amber lines, the perturbative mode structure of the instabilities differs from the region beyond the amber line, as discussed in section III C and particularly in Fig. 4. For pushers, we observe band-like density structures which are not perpendicular to the magnetic field.Pullers in the intermediate region mainly produce a similar density profile as in the region beyond amber line but with finer density structures.Notably, at higher magnetic fields, i.e., larger α e , we observe hatch-like patterns of density with finite width, which are not parallel to the magnetic field.In general, as expected, the instabilities near the line of neutral stability are rather weak and the resulting dynamics only depart very slowly from the steady state.Consistent with the linear stability analysis prediction, we observe distinct patterns for pushers and pullers.More insight into instabilities and their underlying mechanism can be gained by investigating the polarization and self-generated flow fields.In the following, we discuss the prominent features of the long-time density, polarization and flow fields of each of the representative points discussed in Fig. 5 and Fig. 8 corresponding to four distinct instability regimes.

Instability regime (a): traveling bands
The snapshots in Fig. 8 (a) present pattern formation for pushers with dimensionless active stress σ a = −30 and alignment parameter α e = 4.They correspond to the panels (a.1) and (a.2) of Fig. 5, where the linear stability analysis predicts the prevalence of the long wavelength instabilities of wavevectors parallel to the magnetic field.At late stages of the simulations, we observe density modulations in the direction parallel to B, confirming the development of such long wavelength instabilities.The swimmers concentrate in bands perpendicular to the magnetic field spanning the whole transverse dimension of the box while traveling collectively in the field direction.To better understand the origin of these hydrodynamic instabilities, we investigate the corresponding polarization and flow fields.We first focus on a 2D slice of the sample.Fig. 10(a) displays the polarization field superimposed by the density field at a late stage t = 1800.In Fig. 10 (b), the corresponding flow and vorticity fields are presented.We note that concomitant with density modulations, the polarization and flow fields also become non-uniform.The polarization streamlines begin to deflect from straight lines forming bend-like deformations.Such distortions can be understood in terms of bend instability of polarization field similar to those observed in liquid crystals.
Bend fluctuationsconsist of small polarization perturbations which are perpendicular to p 0 while their magnitude is modulated in the direction parallel to p 0 B ≡ ẑ, i. e., p p = p⊥ exp(ik z z) [52,55].Such distortions increase the density in volumes of negative divergence ∇ • p < 0, as expected based on time evolution of density given by Eq. (38).As a result, pushers form dense layers perpendicular to B that migrate parallel to the magnetic field.Bend-like distortions also generate a position-dependent active stress ∝ Q that results in a net active force density in the fluid given by f a ≈ σ a ∇ • pp − 1 3 1 .The active force density leads to alternating flow layers perpendicular to the magnetic field as can be observed from Fig. 10(b).The ensuing vorticity field, encoded by background color, is also modulated in a similar fashion.According to the Faxen's second law [39], the spherical microswimmers is affected by the hydrodynamically induced torque M HD ∝ 1 2 ∇ × u due to the flow vorticity, which rotates the swimmers further away from the magnetic field axis.Thus, the self-generated flow amplifies the bend distortions and renders a uniform homogenous polar phase unstable.This self-amplification would lead to a highly unstable feedback loop if the external torque M B ∝ α e p × B would not eventually counterbalance the hydrodynamic torque.The competi- tion between the alignment and hydrodynamic torques continues until they almost balance each other, hindering further growth of instabilities.As a result, fairly stable patterns at dynamic equilibrium are established.
The emerging picture from a 2D slice of instability snapshot, provides the ground for discussion of 3D patterns.The 3D visualization of the density field shown in Fig. 11(a) is consistent with the picture drawn from a 2D slice.It clearly shows that the pushers concentrate in band-like structures perpendicular to the magnetic field that migrate in the field direction B ≡ ẑ.Now, if we plot the variation of the perpendicular  component of the polarization field averaged in the x-y plane along the z-axis i.e. p ⊥ (z) x,y = (p x (z), p y (z)) x,y , shown in Fig. 12, we observe a helical-like evolution of p ⊥ (z) x,y .For clarity, each of the perpendicular components of polarization, p x (z) x,y and p y (z) x,y are also shown in the p x -z and p y -z planes by red and blue lines, respectively.We note that each of the perpendicular components exhibit a bend-like instability.Therefore, the resultant p ⊥ (z) x,y can be interpreted as a superposition of phase shifted bend deformations of the orientation.The observed behavior is reminiscent of the bendtwist phase predicted for passive chiral or bent-shape liquid crystalline mesogens [56,57] and observed experimentally in achiral molecules [58].Moreover, in hydrodynamic theory of vectorially ordered suspensions of self-propelled particles it was predicted that the coupling between polar or- der and self-generated flow vorticity can lead to formation of bend-twist waves [59].In a bend-twist phase, the polarization orientation vector draws an oblique helicoid, maintaining a constant oblique angle 0 < θ 0 < π/2 with the helix axis z: p = (sin θ 0 cos φ , sin θ 0 sin φ , cos θ 0 ), in which the azimuthal angle varies as φ = 2πz/ p with p being the pitch of the helicoid.
To evaluate if the observed helicoidal pattern is associated with a bend-twist instability, we have extracted the values of the mean polar angle θ p (z) = θ ( p/ρ) x,y and the azimuthal angle of the mean polarization φ p (z) = φ ( p/ρ) x,y averaged in the x-y plane.We find that the mean polarization angle is almost independent of z, θ p (z) ≈ 0.475π, whereas φ p shown in Fig. 13(a) varies nearly linearly with z, apart from the regions close to the box boundaries.These results confirm that polar pushers in an alignment field are prevailed by a bend-twist instability.Moreover, we have also calculated the azimuthal angle of the flow field averaged in the x-y plane φ u = φ ( u x,y ), in Fig. 13(a).φ u also varies almost linearly with z, but shows a clear π/2 phase shift relative to φ p .For comparison, we have also plotted the variation of density averaged over the x-y plane along the magnetic field axis in Fig. 13(b), which clearly shows the modulation of density as a result of band formation.

Instability regime (b): moving pillars
The time series snapshots for pullers at the same alignment parameter α e = 4 and activity strength σ a = 30 are presented in Fig. 8 (b).At late stages, pullers tend to form dynamic pillar-like structures parallel to the external field axis.A 3D rendering of density field is shown in Fig. 11(b).Pullers exhibit a more complex density patterns relative to pushers with the same activity and magnetic field strength.This can be understood in light of the linear stability analysis of homogeneous polar steady state, presented in panels (b1) and (b2) of Fig. 5, which predicts the predominance of long wavelength instabilities with wavevectors both parallel and perpendicular to the magnetic field.Although, our non-linear dynamics simulations display some density undulations in the direction parallel to the magnetic field, we observe that pattern formation is primarily prevailed by the perpendicular perturbations as predicted by the linear stability of the moment equations, Eqs.(38) and (41).
Next, we examine the corresponding polarization and flow fields on a 2D slice of the sample in the y-z plane.Fig. 14(a displays the polarization field superimposed by the density field at a late stage t = 1800 and the corresponding flow and vorticity fields are presented in Fig. 14 (b).We observe that the polarization streamlines significantly deviate from straight lines.The distortions of polarization field can be understood in terms of splay instabilities in the language of liquid crystals [55], see also the appendix for an idealized description of bend and splay distortions.
Splay deformations, similar to bend fluctuations, consist of small polarization perturbations which are perpendicular to p 0 B ≡ ẑ but in this case their magnitude is modulated in directions perpendicular to it, i. e., p p = p⊥ exp(ik ⊥ • x ⊥ ), where k ⊥ ≡ (k x , k y , 0) and x ⊥ ≡ (x, y, 0) [52,55].Again based on Eq. ( 38) for the density moment, splay fluctuations increase the density in volumes of negative divergence ∇ • p < 0 and generate alternating pillar-like flow regions along B as shown in Fig. 14(a); see also the appendix for illustration of an idealized splay distortion.Splay distortions also generate a position-dependent active force density f ≈ σ a ∇ • pp − 1 3 1 in the fluid that result in alternating flow layers parallel and anti-parallel to the magnetic field, see Fig. 14(b).The vorticity field also become heterogenous and induces a hydrodynamic torque M HD ∝ 1 2 ∇ × u.This torque rotates the swimmers further away from the magnetic field axis and renders a uniform homogenous polar phase unstable.The more concentrated regions of pullers, where ∇•p < 0 coincide with regions carrying a flow anti-parallel to B and high self-generated flow vorticity.They result in a net convection anti-parallel to the magnetic field and reduce the mean transport speed [31].Furthermore, we have also shown the polarization field superimposed by the density field and the corresponding flow and vorticity fields of a x-y slice perpendicular to the magnetic field at z = 0.75 in Fig. 15.We note that the perpendicular component of polarization is largest in denser regions.Consistent with our picture from a y-z slice, the regions of large vorticity are correlated with the concentrated regions of swimmers.5, where the linear stability analysis of homogeneous polar steady state predicts predominance of perturbations with a finite wavelength.In both cases, the length scale associated with concentrated regions is finite, k max ≈ 2.5, and smaller than the box size as opposed to cases (a) and (b), where the maximum growth rate occurs at the limit k → 0. However, the angle of wavevector relative to the magnetic field Θ B , for which the the maximum growth rate occurs is different for pushers Θ max B ≈ 30 and pullers Θ max B ≈ 60.In both cases, we see fluctuating concentrated regions which on the average migrate in the direction of magnetic field suggesting that some kind of dynamical aggregates are formed.Consistent with the stability analysis prediction, the morphology of the aggregates are different for pushers and pullers with the identical activity and magnetic field strengths.To gain more insight into similarities and differences between pushers and pullers, we look into the polarization and flow fields in each case.
For pushers in the case (c), concentrated regions form bands with a finite length and a wide range of angles relative to the magnetic field.These patterns are distinct from those of pushers in the case (a) of Fig. 8, where bands are perpendicular to the magnetic field and expand the whole lateral dimension of the box; see Fig. 11(a).Looking into the polar- ization field in a y-z slice shown in Fig. 16(a), we observe very weak deviations from a uniform polarization, whereas the density is notably heterogenous.The generated flow field and its associated vorticity are shown in Fig. 16(b) and they are weaker than the flow and vorticity created in the case (a) presented in Fig. 10(b), which is predominated by the bendtwist instability.Examining the polarization field superimposed by density in a x-y slice perpendicular to the magnetic field shown in Fig. 17(a), we find that the deviations of polarization field from the B ≡ ẑ occur at concentrated regions.In other words, the finite wavelength instabilities are dominantly density driven and the polarization distortions merely stem from density perturbations.In the language of linear stability analysis, the predominant mode of perturbation is given by the mode shape ψ 0 0 (k max )e ik max •x+λt ∝ ρ p (k max , x,t), where k max corresponds to the wavevector with the largest growth rate shown in Fig. 5(c).Likewise, from Fig. 17(b), we note that the self-generated flow velocity component perpendicular to the magnetic field (u x , u y ) is rather weak and it only becomes considerable in concentrated regions.Unlike the case (a), the flow field has an appreciable component along the magnetic field (parallel or anti-parallel) as encoded by red and blue colors in Fig. 17(b).
For the pullers in the stability regime (d), the long-time density pattern resembles that of pullers with moderate strengths of the activity and magnetic field in case (b).However, concentrated regions consist of of finite-sized pillar-like aggregates in contrast to the case (b), where pillar-like dense regions span the whole box dimension in the field direction, ver- ifying the predominance of smaller wavelength density fluctuations.Moreover, finite-sized concentrated regions are on average not parallel to the magnetic field and have a wider orientation distribution.Inspecting the polarization field superimposed by the density field in a y-z slice perpendicular to the magnetic field shown in Fig. 18(a), we notice that polarization field displays some splay deformations.However, its distortions are weaker in comparison to the case (b).Likewise, the self-generated flow field is very similar to the case (b) and we observe a notable vorticity field in concentrated regions.Looking into the polarization and velocity fields in a x-y slice perpendicular to B ≡ ẑ, we find that similar to case (c), the deviations from a uniform polarization occurs at concentrated regions, which lead to a very heterogenous flow field as shown in Fig. 18(b).Despite the similarities of polarization and flow field with the case (b), the predominant mode of perturbation in the linear stability is the density mode similar to the case (c).It is given by ψ 0 0 (k max )e ik max •x+λt ∝ ρ p (k max , x,t), where k max corresponds to the wavevector with the largest growth rate shown in Fig. 5(d).

VI. DISCUSSION AND CONCLUDING REMARKS
We have presented a continuum kinetic model for active suspensions of weakly magnetic spherical particles in an external field.The model is based on first principles, namely, a conservation equation for the particle configuration distribution in an alignment field, coupled to the Stokes equation for the fluid flow which incorporates stress contributions steming from activity and alignment torque.It is applicable to moderately dilute suspensions of magnetotactic bacteria or artificial magnetic microswimmers with a small magnetic moment and focuses on the interplay between hydrodynamic interactions arising from self-generated flow and external alignment torque.We investigated the nature of hydrodynamic instabilities and emergent pattern formation by combining linear stability analysis and the full numerical solution of kinetic model equations.
We first performed a linear stability analysis of steady state solution of the model, which corresponds to a homogenous polar distribution function ψ 0 .The stability analysis of steady state as a function of activity and magnetic field strengths reveals that a uniformly polarized suspension becomes unstable for moderate magnetic field and sufficiently large activity strengths for both pushers and pullers.Based on the dispersion relation of the maximum growth rate, we have drawn a non-equilibrium phase diagram as presented in Fig. 4. We recognize four distinct instability regimes.For moderate activity and field strengths, the long wavelength instabilities dominate both pusher and puller suspensions.However, the nature of instabilities are different for the two types of swimmers.Pushers are dominated by wave perturbations parallel to the field, whereas pullers are unstable with respect to both parallel and perpendicular wave perturbations.For stronger activities and magnetic fields, the wavevector with the largest growth rate has a finite wavelength and its angle with the field differs for pushers and pullers with the same activity and field strengths.These instability regimes are driven by density fluctuations as opposed to long wavelength instabilities which are driven by the orientational fluctuations.Increasing the magnetic field strength further, the alignment torque is strong enough to overcome the hydrodynamically induced torque.As a consequence, the homogenous polar state becomes stable again and we observe a reentrant hydrodynamic stability.
Next, we obtained the dynamical equations for the first three orientational moments, i.e., density, polarization and nematic field, imposing suitable closure approximations.Moment equations, although less accurate, provide us with new insights into the nature of instabilities.As can be seen from the moment equations, Eqs. ( 38), ( 41) and ( 45), density, polarization and nematic fields are coupled to each other.This implies that any heterogeneity in one of them generates a heterogeneity in the other fields leading to a feedback loop until a new dynamical equilibrium is reached.Linear stability analysis of moment equations for uniform density and polarization fields predicts the predominance of long wavelength instabilities with wavevectors parallel to the alignment field for pushers and wavevectors perpendicular for pullers at moderate magnetic fields.Based on these results, we deduce that pushers are prevailed by bend deformations, whereas pullers are predominated by splay distortions.These findings are in agreement the linear stability analysis of the steady distribution function ψ 0 for a large region of stability diagram, where long wavelength instabilities prevail the system, although perturbation of ψ 0 , equivalent to considering the full hierarchy of moments, predicts the predominance of both long wavelength parallel and perpendicular perturbations for pullers.Moreover, the coarse-grained approximate moment equations do not capture the finite wavelength instability regimes at higher magnetic fields and activity strengths.
To evaluate the accuracy of predictions of the linear stability analysis, we investigated the numerical solution of kinetic model equations.Numerical simulations show very good agreement with predictions of linear stability analysis for the borderlines of instability.Although linearly unstable modes do not capture the full non-linear dynamics, many aspects of the dynamics observed in simulations can be understood in the light of the stability analysis.According to Fig. 4, for a large region in the parameter space the most unstable mode for pushers is parallel to the external field, whereas for pullers both parallel and perpendicular unstable modes dominate the system.Simulations show that indeed in a large part of the unstable region predominant modes of instability for pushers are bend-twist distortions of the polarization field.For the pullers, splay deformations prevail the pattern formation suggesting that the perpendicular perturbation is the predominant mode of deformation.As a consequence, we observe distinct patterns for the two kinds of swimmers: traveling bands perpendicular to the magnetic field for pushers and pillar-like concentrated regions parallel to the field for pullers.As discussed in our prior work [31], the deflections of polarization field lead to a reduction of the average polarization and mean transport speed.In the regions of stability diagram of Fig. 4, where the maximum growth rate occurs at finite wavelengths, we observe finite-sized concentrated regions suggesting formation of dynamical aggregates in external field.However, the morphology of these regions is different for pushers and pullers in agreement with predictions of linear stability analysis.
We conclude by pointing out a few limitations of the present model and future directions.Our results are obtained in the limit of negligible magnetic interactions and only consider the interaction of a single particle with a mean-field flow.This limits the validity of our model to moderately dilute suspensions to magnetic swimmers with weak dipole moments such as magnetotactic bacteria.Nevertheless, we believe that the present model captures most salient features of interplay between hydrodynamic interactions and external field in not so concentrated active suspensions.For instance, band formation observed for pushers are in agreement with experimental findings of magnetotactic bacteria at moderate field strengths B ∼ 3 mT [18].For synthetic magnetic microswimmers with larger magnetic dipole moments or dense suspensions, the magnetic dipolar interactions alone can lead to clustering instabilities [60] and the interplay between long-range magnetic and hydrodynamic interactions on development of instabilities merits further investigations.Moreover, the role of swimmerswimmer correlations [61], and near-field hydrodynamic interactions in more concentrated suspensions remains an open question.Finally, our results show that a sufficiently strong alignment field can overcome hydrodynamic instabilities calling for further exploration of controlling the collective behavior and transport of active matter in various external fields.pushers u = u ⊥ (x )x ⊥ and parallel to the field for pullers u = u (x ⊥ )x .To very good approximation, it can be described by u {⊥, } ({x , x ⊥ }) ≈ c 1 sin{x , x ⊥ } + c 2 sin 3{x , x ⊥ }, in which the braces {} evaluate to the first entry for pushers, the second entry for pullers, and c 1 , c 2 ∈ R are some numerical prefactors.
The polarization field and the associated flow and vorticity fields of pushers with bend deformations and pullers with splay distortions in the alignment field are shown in Fig. 20 and Fig. 21, respectively.The subplots (a) shows the streamlines of the polarization field p with its divergence color encoded in the background.The swimmers concentrate in the red areas where the divergence is negative.As a consequent, pushers form band-like regions perpendicular to the field, whereas pullers concentrate in pillar-like dense regions parallel to the field.The subplots (b) depict a vector plot of the self-generated flow fields as a result of the bend and splay deformations of the polarization fields of pushers and pullers, respectively.In both cases, we observe alternating flow layers modulated in the direction perpendicular to field.For the pushers, the flow velocities with alternating directions are perpendicular to the alignment field, while for pullers the the alternating flow velocities are parallel and anti-parallel to the field.In subplots (c), the flow fields' vorticity fields ∇ × u are shown, where the red color encodes a counter clockwise rotation and the blue color a clockwise rotation.We observe the alternating clockwise and anticlockwise vorticity fields are also modulated in directions parallel perpendicular to the field for pushers and pullers, respectively.According to the Faxen's law, the flow vorticity induces effectively a hydrodynamic torque M HD ∝ 1 2 ∇ × u, which rotates the particles away from the magnetic field direction and competes with the magnetic torque.A homogenous polar steady state becomes unstable due to these competing torques.
In the idealized case, both torques are easy to calculate and are plotted in arbitrary units in Fig. 22 for a bend-deformation (it is qualitatively the similar for splay-deformations).The mean, dimensionless external torque, given by M B = α e p × B nearly fully balances the M HD hindering further growth of bend deformation.As a result, a stable dynamical pattern is established.The competition between the two torques becomes apparent looking into the rotational drift velocity given by Eq. ( 15) which can be equivalently written as ṅ = α e n × B + 1 2 (∇ × u) × n. (A1) Under conditions that both terms compensate each other, the bracket vanishes and the orientation n does not change any more.

FIG. 1 .
FIG.1.A schematics of a spherical microswimmer with unit orientation vector n and swimming speed U 0 n, carrying a magnetic dipole moment µ.The dynamics of the swimmer is influenced by the local flow U and its vorticity Ω = ∇ × U, and an external magnetic flux density B.

FIG. 2 .
FIG. 2. (a) Angular distribution of the homogeneous steady state ψ 0 (θ ) for different values of the external alignment parameter α e = µB ξ r D r .(b) The total polarization p 0 of the steady state ψ 0 as a function of α e ∝ B.

B. Linearized equations and eigenvalue problem 1 .
Linear perturbation of the base state

40 FIG. 3 .
FIG. 3. (a) Magnitude of change of the corresponding eigenvectors (in the basis of spherical harmonics) with respect to a truncation order j max , diff( ψ, j max ) = ψ[ j max ] − ψ[ j max − 1] .The remaining parameters are fixed to σ m = 0.4α e , d t = 3 × 10 −6 .(b) Change of the largest eigenvalues diff(λ , j max ) = |λ [ j max ] − λ [ j max − 1]| as a function of the number of included modes j max on a logarithmic scale.(c) The largest growth rate Re λ max [ j max ] as a function of truncation order j max .

FIG. 4 .
FIG.4.Stability diagram of the steady state ψ 0 (α e ) given by Eq.(20) as a function of the dimensionless active stress σ a and alignment parameter α e ∝ B; while setting σ m = 0.4 α e and d t = 3 × 10 −6 .The borderline of neutral stability (red dash-dotted line) is calculated by finding Re λ max (k) = 0.The dashed amber lines correspond to the cases where Re λ max (k || ) = 0 for pushers and Re λ max (k ⊥ ) = 0 for pullers.The solid blue lines represent Re λ max (k) = 0 based on the linear stability analysis of density and polarization fields from truncated moment equations for wavevectors parallel and perpendicular to the magnetic field.On the right side, the corresponding polarization p 0 of the steady state ψ 0 (α e ) is plotted.For the points marked by crosses, the behavior of growth rate and pattern formation are further discussed in the paper.

FIG. 8 .
FIG. 8. Representative snapshots of density field projections averaged along the y-axis from 3D non-linear simulations at different time steps, as shown on the snapshots, for magnetic swimmers with activity strength and alignment parameter values (a) σ a = −30 and α e = 4 (b) σ a = 30 and α e = 4, (c) σ a = −40 and α e = 19 and (d) σ a = 40 and α e = 19.The colors encode the probability density integrated in the y-direction ρ(x, z) = ∆y∑ y ρ(x, y, z).In all simulations, the translational diffusion is fixed to d t = 3 × 10 −6 and the alignment stress is varied along the external field as σ m = 0.4α e .

FIG. 9 .
FIG. 9. Density projections for different (α e σ a ) values superposed by the lines of neutral stability from linear stability analysis.The red lines marks the onset of instability calculated by the linear stability analysis probing all possible orientations of wavevector.The amber dashed lines correspond to predictions for neutral stability (Re λ max (k) = 0) when considering only perturbations parallel k || (pushers) and perpendicular k ⊥ (pullers) to the external field.The translational diffusion coefficient is fixed to d t = 3 × 10 −6 and the alignment stress is varied with the external field as σ m = 0.4 α e .

8 FIG. 10
FIG. 10.(a) Streamlines of the polarization field (p y , p z ) of a slice in the y-z-plane at x = 0.75, demonstrating a characteristic bend fluctuations for a suspension of pushers with dimensionless active stress σ a = −30 and alignment parameter α e = 4 in a magnetic field pointing in z-direction.The color encoded local density ρ/V × 10 6 is shown in the background.(b) The corresponding flow field (u y , u z ) is represented as vector arrows where the length of the vector is weighted by its magnitude.The flow vorticity, responsible for the hydrodynamically induced particle rotation, is color encoded in the background.Red colors correspond to a vorticity vector pointing out of the plane (counter clockwise rotation).Lines of constant density are overlaid as contour lines on top to guide the eyes.

FIG. 11 .
FIG. 11. 3D volumetric rendering of the density field of a pusher (a) and puller (b) at t = 1800 corresponding to the last time-step presented in panels (a) and (b) in Fig. 8, respectively.

FIG. 12 .
FIG.12.3D representation of the radial component of the mean orientation p/ρ) x,y revolving around the magnetic field axis in a helical twist pattern.The red and blue lines at the wall represent the respective projections in the are plotted, the black lines mark the of the density bands.

FIG. 13 .
FIG. 13.(a) Mean azimuthal angles around the magnetic field axis for the orientation field φ p = φ ( p/ρ) x,y and the flow field φ u = φ ( u x,y ). and (b) mean density along the magnetic field axis marking the position of two bands perpendicular to the magnetic field (compare with FIG.8(a)) for a representative snapshot of pushers with a dimensionless active stress σ a = −30 and alignment parameter α e = 4 at t = 1800.The magnetic field points along the z-axis.

8 FIG. 14
FIG. 14.(a) Streamlines of the polarization field (p y , p z ) of a slice in the y-z-plane at x = 0.75, demonstrating a characteristic splay patterns for a suspension of pullers with dimensionless active stress σ a = 30 and alignment parameter α e = 4 in a magnetic field pointing z-direction.The color encoded local density ρ/V × 10 6 is shown in the background.(b): The corresponding flow field (u y , u z ) is represented as vector arrows, where the length of the vector is weighted by its magnitude.The flow vorticity, leading to hydrodynamically induced particle rotation, is color encoded in the background.Red colors correspond to a vorticity vector pointing out of the plane (counter clockwise rotation).Lines of constant density are overlaid as contour lines on top to guide the eyes.

FIG. 15
FIG. 15. (a): Streamlines of the polarization field (p x , p y ) of a slice in the x-y-plane at z = 0.75, where the size of arrows shows the relative magnitude of polarization for a suspension of pullers with dimensionless active stress σ a = 30 and alignment parameter α e = 4 for a magnetic field pointing in z-direction.The color encoded local density ρ/V × 10 6 is superimposed in the background.(b): The flow field in the plane (u x , u y ) is represented as black arrows, where the length of the vector is weighted by its magnitude.The z-component of flow field, u z is color encoded, where a red color denots a flow in the direction of the external magnetic field and a blue color opposed to it.Lines of constant density are overlaid as contour lines on top to guide the eyes.

3 .
Instability regimes (c) and (d): finite-sized concentrated regions The snapshots in Fig. 8 (c) and (d) show the evolution of density profiles of pushers with α e = 19 and |σ a | = 40.They correspond to the panels (c) and (d) of Fig.

8 FIG. 16
FIG. 16.(a) Streamlines of the polarization field (p y , p z ) of a slice in the y-z-plane at x = 0.75, for a suspension of pullers with dimensionless active stress σ a = −40 and alignment parameter α e = 19 in a magnetic field pointing in z-direction.The color encoded local density ρ/V × 10 6 is shown in the background.(b): The corresponding flow field (u y , u z ) is represented as vector arrows, where the length of the vector is weighted by its magnitude.The flow vorticity, leading to hydrodynamically induced particle rotation, is color encoded in the background.Red colors correspond to a vorticity vector pointing out of the plane (counter clockwise rotation).Lines of constant density are overlaid as contour lines on top to guide the eyes.

50 FIG. 17
FIG. 17. (a): Streamlines of the polarization field (p x , p y ) of a slice in the x-y-plane at z = 0.75, where the size of arrows shows the relative magnitude of polarization for a suspension of pushers with dimensionless active stress σ a = −40 and alignment parameter α e = 20 for a magnetic field pointing in z-direction.The color encoded local density ρ/V × 10 6 is superimposed in the background.(b): The flow field in the plane (u x , u y ) is represented as black arrows, where the length of the vector is weighted by its magnitude.The z-component of flow field, u z , is color encoded with a red color denoting a flow in the direction of the external magnetic field and a blue color opposed to it.Lines of constant density are overlaid as contour lines on top to guide the eyes.

8 FIG. 18
FIG. 18.(a) Streamlines of the polarization field (p y , p z ) of a slice in the y-z-plane at x = 0.75, exhibiting weak splay distortions for a suspension of pullers with dimensionless active stress σ a = 40 and alignment parameter α e = 19 in a magnetic field pointing in z-direction.The color encoded local density ρ/V × 10 6 is shown in the background.(b): The corresponding flow field (u y , u z ) is represented as vector arrows, where the length of the vector is weighted by its magnitude.The flow vorticity, leading to hydrodynamically induced particle rotation, is color encoded in the background.Red colors correspond to a vorticity vector pointing out of the plane (counter clockwise rotation).Lines of constant density are overlaid as contour lines on top to guide the eyes.

FIG. 19
FIG. 19.(a): Streamlines of the polarization field (p x , p y ) of a slice in the x-y-plane at z = 0.75, where the size of arrows shows the relative magnitude of polarization for a suspension of pullers with dimensionless active stress σ a = 40 and alignment parameter α e = 20 for a magnetic field pointing in z-direction.The color encoded local density ρ/V × 10 6 is superimposed in the background.(b): The flow field in the plane (u x , u y ) is represented as black arrows, where the length of the vector is weighted by its magnitude.The z-component of flow field, u z , is color encoded with a red color denoting a flow in the direction of the external magnetic field and a blue color opposed to it.Lines of constant density are overlaid as contour lines on top to guide the eyes.

FIG. 22 .
FIG.22.Illustration of the competition between a (virtual) hydrodynamic torque M HD ∝ 1 2 ∇ × u and an external magnetic torque M B ∝ α e n × Be for bend and splay instabilities in arbitrary units.The instabilities cause the flow vorticity to grow while rotating the particles, until the rotation is compensated by the external torque (amber).