Orientation mixing in active suspensions

A BSTRACT . We study a popular kinetic model introduced by Saintillan and Shelley for the dynamics of suspensions of active elongated particles where the particles are described by a distribution in space and orientation. The uniform distribution of particles is the stationary state of incoherence which is known to exhibit a phase transition. We perform an extensive study of the linearised evolution around the incoherent state. We show (i) in the non-diffusive regime corresponding to spectral (neutral) stability that the suspensions experience a mixing phenomenon similar to Landau damping and we provide optimal pointwise in time decay rates in weak topology. Further, we show (ii) in the case of small rotational diffusion ν that the mixing estimates persist up to time scale ν − 1 / 2 until the exponential decay at enhanced dissipation rate ν 1 / 2 takes over. The interesting feature is that the usual velocity variable of kinetic models is replaced by an orientation variable on the sphere. The associated orientation mixing leads to limited algebraic decay for macroscopic quantities. For the proof, we start with a general pointwise decay result for Volterra equations that may be of independent interest. While, in the non-diffusive case, explicit formulas on the sphere allow to conclude the desired decay, much more work is required in the diffusive case: here we prove mixing estimates for the advection-diffusion equation on the sphere by combining an optimized hypocoercive approach with the vector ﬁeld method. One main point in this context is to identify good commuting vector ﬁelds for the advection-diffusion operator on the sphere. Our results in this direction may be useful to other models in collective dynamics, where an orientation variable is involved.

exponential convergence for analytic data, to polynomial convergence for Sobolev data.When diffusion is added, leading to the small scales created by the filamentation allow for an acceleration of the diffusive process; this leads to an enhanced dissipation on time scales ν −1/3 shorter than the usual ν −1 .Identifying these phenomena in more complex kinetic equations, either of transport or weakly dissipative type, linearly and nonlinearly, has attracted a lot of attention.A main example is the analysis of Landau damping in the Vlasov-Poisson equation of plasma physics, with pioneering works [26,30].More recent works were dedicated to weakly dissipative cases, trying to exhibit both Landau damping and enhanced diffusion at the same time [2,5].
Another field where similar mixing phenomena are identified is the one of collective dynamics.A famous example in this direction is the Kuramoto model, a system of ODEs describing the dynamics of coupled oscillators, known to exhibit phase transition, from incoherence to synchronization.At the mathematical level, the mean-field limit of this system of ODEs leads to a kinetic equation, with unknown f (t, θ, ω) describing the fraction of oscillators with phase θ ∈ T and natural frequency ω ∈ R. Depending on a bifurcation parameter quantifying the intensity of the coupling, solutions converge weakly either to the uniform distribution (incoherence), or to the so-called partially locked states modeling synchronization (and containing Dirac masses in θ).It turns out that this convergence is again due to a mixing process, as noticed for the first time in [37].Such phenomenon is now fully confirmed mathematically, both linearly and nonlinearly [7,15,20,4,17,16].
In this paper, we will study mixing properties of a popular model describing a dilute suspension of elongated active particles in a Stokes flow.
1.1.The model.In [36], see also [35] for a review, D. Saintillan and M. Shelley have introduced a system of PDEs to describe the dynamics of a dilute suspension of elongated active particles in a viscous Stokes flow.The word "active" refers to the fact that they convert chemical energy into mechanical work.A typical example are bacteria, which are able to swim and create stress through the use of their flagellas.
To write down the model, the first step is to consider a collection of N elongated particles (ellipsoids) immersed in a Stokes flow.They are described by the position of their center of mass x i ∈ T 3 L = (R/(LZ)) 3 and their director p i ∈ S 2 , 1 ≤ i ≤ N .The choice of a torus of size L for the spacial domain instead of a real container is for mathematical convenience: it introduces a typical length scale L without problems related to the boundary of the domain.Neglecting interaction between the particles, one can use the work of Jeffery about a single passive ellipsoidal particle in a Stokes flow, see [24,38].If the typical size of the particle is very small compared to L, the torque experienced by the particle i is approximated by P p ⊥ i (γE(u) + W (u)) | x i p i , where P p ⊥ i = I − p i ⊗ p i is the projection orthogonal to p i , while E(u) and W (u) are the symmetric and skew-symmetric parts of ∇u.The parameter γ is related to the aspect ratio of the ellipsoid, with γ = 1 in the limit case of a rod.Note that in the work of Jeffery, u refers to the unperturbed velocity field, that is the one forgetting the particles.It is thus defined everywhere on T 3 L , especially in x i .The dynamics is then given by ẋi = u(x i ) + U 0 p i , ṗi = P p ⊥ i γE(u) + W (u) | x i p i + possible Brownian noise.
The first equation describes the velocity of the particles as a sum of two contributions: the first one corresponds to advection by the velocity field u of the flow.The second one corresponds to swimming along the director p i , at constant speed U 0 > 0. The second equation describes the rotation of the particle, which is known in the absence of noise to experience periodic trajectories (Jeffery's orbits).As far as we understand, these orbits are not really observed in experiments, and it is more accurate to perturb them by adding additional rotational Brownian motion, of small amplitude ν ≪ 1.On the other hand, we neglect Brownian motion in x, that is translational diffusion.As pointed out in [36, Page 7], it is not expected to play a big role in the linear stability properties of the model.We stress though that it could be included in the analysis below, and all the results would then hold independently of the strength of this diffusion.FIGURE 1. Illustration of pusher and puller active particles moving to the right.Depending on configuration, the force for propulsion is differently applied on the fluid, which yields a different induced flow around the active particle [35].
Besides the swimming velocity U 0 p i , another feature of active suspensions that needs to be retained is the stress that they create on the fluid.Saintillan and Shelley depart here from the work of Jeffery by including this effect at the level of the Stokes equation on u.The extra stress due to particle i can be thought as a dipole: the sum of two opposite and close point forces, along the direction of p i : typically, for bacteria, one point force is centered at the body of the particles, while the other is centered at the flagella, see Fig. 1 distinguishing between two kinds of bacteria: pushers and pullers.One ends up with where this appearance of the divergence operator in front of a Dirac mass is typical of a dipole.Depending on the orientation of the dipole, the sign of α N can be positive (pullers) or negative (pushers).We refer to [6] for a more rigorous derivation.
Eventually, with the right scaling of α N , and performing a (formal) mean-field limit, one obtains the Saintillan-Shelley model, with two unknowns: • Ψ = Ψ(t, x, p), the distribution of particles in space and orientation, • u = u(t, x), the velocity field of the fluid.It reads Ψ(t, x, p) p ⊗ p dp, where we remind that 3) The first equation on Ψ reflects transport along the stream lines associated to inspired by the dynamics of the particles mentioned above.It also contains some diffusion term ν∆ p Ψ (ν ≪ 1), where ∆ p refers to the Laplace-Beltrami operator associated to possible Brownian rotational noise.The other equations describe the Stokes flow, which is forced by the extra stress div σ, with σ given by the marginal in p of αΨp ⊗ p.
Let us notice that this model is very close to the Doi model [18] which describes passive suspensions of rodlike polymers.In this setting, U 0 = 0 (passive particles), and α > 0 is proportional to the Boltzmann constant and the temperature: the tensor σ models viscoelastic stress.While local existence and uniqueness of smooth solutions of both the Doi and the Saintillan-Shelley models can be obtained by standard methods, the question of global in time well-posedness is harder.In the case of the Doi model, global well-posedness was proved in [32] for ν > 0, and in [9] for ν ≥ 0. Actually a look at the proof of [9] shows that it is still valid for α of any sign.However, for non-zero U 0 , as far as we know, global well-posedness is unknown (except for the addition of diffusion in x, see [6]).It seems to be an interesting open problem, especially when ν = 0: indeed, one can check that the force field in the equation for P p ⊥ [(γE(u) + W (u))p] has the same regularity in x as σ.This is to be compared with Vlasov-Poisson equation, where the force field has one degree of regularity more than the density ρ.In the present paper, we will not contribute to this well-posedness theory, but will rather investigate qualitative properties at the linear level.

Asymptotic stability of incoherence: main results. A class of equilibria of interest is given by
of which the isotropic suspension Ψ iso = 1/4π is a particular case.By analogy with the Kuramoto model, we will call Ψ iso the incoherent state, as it reflects zero alignment between the orientation of the particles.
The linearized system around (Ψ iso , u S = 0) reads Since E(u) is symmetric, W (u) is skew-symmetric and u is divergence-free, we have Hence the equations become (1.9) A partial analysis of system (1.9) is performed in [36,23] • in the case ν = 0, looking for unstable eigenmodes, the authors manage to calculate an explicit dispersion relation.They show that in the case of suspensions of pullers (α > 0), no unstable eigenvalue exists.As pointed out by the authors, this is consistent with the fact that for the full nonlinear system (1.9), for any α > 0, ν ≥ 0, the relative entropy of any solution ψ with respect to ψ iso decays.However, the situation changes completely for pushers (α < 0).In this case, unstable eigenvalues exist at low enough x frequencies.In other words, the length L of the torus is a bifurcation parameter: there is a critical value L c , computed numerically in [36], such that for L < L c , there exists no unstable eigenvalue, while for L > L c , there exists at least one.Above this threshold, the incoherent state loses its stability, some partial alignment of the ellipsoidal particles is observed numerically, while the corresponding velocity field develops patterns that are favourable to mixing of passive scalars advected by the flow.See the nice recent work [31] for more on the nature of the bifurcation process.
• in the case of small ν > 0, no explicit dispersion relation is available.But numerical computations in [36], confirm the picture given at ν = 0, with some threshold close to L c .This numerical observation will be confirmed rigorously here.
Our focus in the present paper is on what we call the incoherent regime, that is L < L c , both in the case ν = 0 and ν > 0. Again, some theoretical observations are already contained in [36], see also [23].For ν = 0, the absence of unstable eigenvalue for the linearized operator in (1.9) does not imply automatically linear stability, due to the possible unstable continuous spectrum.Numerical simulations in [23] show high frequency oscillations but no instability.Moreover, some decay at rate (kt) −2 , where k is the x-frequency of the perturbation, is seen for some moment of the solution ψ with respect to p. Stability of the incoherent state is conjectured on the basis of these simulations.Regarding ν > 0, stability is confirmed by simulations, but no clear rate of convergence with respect to ν is given.
Our aim is to clarify most of these aspects, through a detailed mathematical study of the linearized equation (1.9).Since x ∈ T 3 L , we take the Fourier transform of ψ and similarly for u k , the Fourier transform of u.The Fourier transform of the equation is then Note that for k = 0, the equation reduces to the simple heat equation ∂ t ψ 0 − ν∆ p ψ 0 = 0 so that we restrict to k ̸ = 0.Moreover, through the change of variables the system becomes where k ∈ S 2 , ε = ±1 with ε = 1 for pullers, ε = −1 for pushers.(1.14)In this rescaled setting, our main results are the following: THEOREM 1 (Inviscid case, pointwise decay through mixing).Let ν = 0, There is an absolute constant Γ c ∈ (0, +∞], with Γ c = +∞ for ε = 1, and THEOREM 2 (Diffusive case, mixing persists for small ν up to time ν −1/2 ).Let ν ∈ (0, 1), ) After the time ν −1/2 , we have the mixing estimates with an additional log factor until the enhanced dissipation takes over.THEOREM 3 (Diffusive case, enhanced dissipation at timescale ν −1/2 ).Under the assumptions of Theorem 2, there exist M, η, ν 0 > 0 depending on Γ, and C = C(δ, Γ) such that for 0 < ν < ν 0 and t ≥ ν −1/2 it holds that ) A few remarks are in order.
Remark 1.1.On the torus the largest length scale is L, so that in the unscaled equation the original wavenumber |k| is at least 2π/L.This leads to the upper bound Hence the stability condition on Γ can be translated into a maximal size of the torus.
Remark 1.2.Theorem 1 is the expression of an inviscid damping: it leads to a O(t −2 ) decay for the velocity field u k .Notice that u k involves ψ k through an average in p: this average allows to take advantage of the mixing phenomenon.As a byproduct, one has a O(t −1 ) decay of ψ k itself in weak topology (negative Sobolev space).The difference in the rate of decay is very much related to the special structure of u k in terms of ψ k .We stress that these polynomial rates cannot be improved, even taking analytic initial data.This is a major difference with Landau damping, and is related to the fact that the orientation variable p ∈ S 2 replaces the velocity variable v ∈ R 3 .
Remark 1.3.In our proof of Theorem 1, the stability constant Γ c that we exhibit when ε = −1 is sharp: it means that for Γ > Γ c there exist eigenmodes with u k growing.Concretely, condition Γ < Γ c is equivalent to the fact that some dispersion relation has no root in the unstable half-plane: Our rigorous analysis of this dispersion relation, including the identification of Γ c , is inspired by the work of Penrose [33] on the stability of plasmas.It substitutes for the numerical analysis carried in [36].
The fact that ν −1/2 is a natural time threshold for our problem is confirmed by Theorem 3: solutions of the equation experience an exponential decay after this typical threshold.This exponential decay before the natural diffusive timescale ν −1 reflects the well-known phenomenon of enhanced dissipation [3,10,12,13], evoked in the introduction, although one can notice again a strong qualitative difference between variable p ∈ S 2 and variable v ∈ R 3 .In the latter case, the enhanced dissipation would hold with typical time scale ν −1/3 .Here, the typical time scale is much longer, which creates strong mathematical difficulties in showing mixing up to this time scale, and enhanced diffusion afterwards.
For the transition time until the enhanced dissipation takes over, we show the persistence of the mixing estimates with additional log factors.In advection-diffusion problems, a diffusion-limited mixing behavior is often observed [29].That is, inviscid mixing does not persist as t → ∞, rather the ratio of the Ḣ−1 to the L 2 norm converges to a positive constant.This corresponds to the existence of a characteristic filamentation length scale, often referred to as the Batchelor scale in the physics literature.
Remark 1.5.While completing our manuscript, we noticed the release of the interesting preprint [1], about the same Saintillan-Shelley model (except for the introduction of an additional diffusion in variable x).
Although paper [1] and ours share a few common features, we believe that they are different enough to provide distinct and valuable insights into the stability of active suspensions.
Regarding mixing, which is the focus of the present paper, [1] only contains a weaker version of our Theorem 1, showing that under condition (1.22), a L 2 in time stability estimate Our extensive discussion of (1.22), see Remark 1.3, as well as our optimal pointwise in time decay estimates are not covered.We stress that these pointwise estimates are based on the general Theorem 5 on Volterra equations, which is of independent interest.More importantly, [1] does not contain any analogue of our Theorem 2, which is the heart of our paper, and requires completely different arguments from the inviscid case.
On the other hand, [1] contains the derivation of Taylor dispersion estimates when x ∈ R d , and two nonlinear stability results (diffusion in x is crucial there).First, in the case of pullers (ε = 1), using an approach à la Guo [22], the authors prove nonlinear stability of the incoherent state Ψ iso , both for x ∈ T d and x ∈ R d (but without enhanced dissipation).Second, in the case of pushers, they prove nonlinear stability of the incoherent state with enhanced diffusion, under the stringent assumption Γ = o(ν 1/2 ), which allows to treat the model as a perturbation of the advection-diffusion equation.With regards to these last two results, and although we restrict here to a linear analysis, our Theorem 3 is a significant progress.
1.3.Key ideas.The evolution (1.13) can be split as where is the advection-diffusion operator, and u k is the low-dimensional linear map from the kinetic distribution to the macroscopic velocity field given by Finally, L2,k is the vector field independent of time defined by Here, we have used that Applying the linear map u k then yields the Volterra equation with the source (1.29) and the matrix kernel K k defined by (1.30) In the proof of Theorems 1, 2 and 3 the key point is to work on the Volterra equation (1.28) and to show that u k decays like t −2 at infinity.Once this is obtained, to go back to the evolution of ψ k and estimate its decay in negative Sobolev norms is direct.Regarding the decay of u k , there are two main steps: i) To show that the source term U k and the kernel K k decay like in the corresponding theorem.
ii) To show that under appropriate conditions on Γ, this decay passes to the solution u k of the Volterra equation.Let us stress that the mathematical questions raised by the second step are not specific to our model, and can be asked for any Volterra equation.They will be examined in Section 2. Obviously, before relating the decay rate of the solution to the decay rates of the source and the kernel, a prerequisite is the stability of the solution.Using the classical theory of Volterra equation, one knows that it is equivalent to the spectral condition ) where L is the Laplace transform.Under this condition, we achieve Step ii) by proving Theorem 5 in Section 2. We provide a short proof, that uses the notion of Volterra kernel of type L ∞ and the underlying structure of Banach algebra of this class of kernels As regards Step i), there is a main difference between ν = 0 and ν > 0. In the inviscid case ν = 0, considered in Section 3, the evolution e tL 1,k can be solved explicitly: the kernel and sources are given by Fourier transforms on the sphere, whose decay properties are well-known, and analyzed through stationary phase arguments.This implies that U k and K k are O(t −2 ).Moreover, the spectral stability condition can be fully understood through an analysis à la Penrose, see Section 3, completing the proof of Theorem 1.
For the diffusive case ν > 0, we cannot solve the evolution e tL 1,k explicitly and the main effort in proving Theorem 2 and 3 is to obtain optimal mixing estimates on the advection-diffusion equation.This is performed in Section 4. In terms of the advection-diffusion evolution e tL 1 , the first result is an enhanced dissipation result, which can be proved by adapting hypocoercive methods: PROPOSITION 1.6.There exist constants C 0 , ε 0 , ν 0 > 0 with the following property: for every ν 0 > ν > 0 and ψ in the advection-diffusion evolution is bounded for all t ≥ 0 as (1.32) The main novelty is that the inviscid mixing estimates persists until the enhanced dissipation takes over: PROPOSITION 1.7.There exist ν 0 > 0 and M ≥ 2 such that for c, δ > 0 there exists a constant C with the following bounds: For any ν ≤ ν 0 , any scalar function F and t ≤ ν −1/2 it holds that and and The key idea is the use of the vector field method.This method, introduced for the analysis of decay properties of wave equations [25], was used recently in the context of Vlasov type equations on R d [5].It allowed to exhibit mixing when 0 < ν ≪ 1, despite the loss of explicit representation formula.Let us just mention at this stage that the key point is to construct good vector fields, meaning that they commute to the advection-diffusion operator ik • p − ν∆ p .In the Euclidean setting (1.1), the natural choice J = ∇ v + itk works well over the time scale ν −1/3 , after which enhanced dissipation dominates.On the sphere this has no good analogue as the obvious generalization J = ∇ p + i∇(p • k)t, when applied to a solution ψ of the advection-diffusion equation, creates commutators of the form νt∇ψ, with a time integral that can not be controlled over the large time scale ν −1/2 .We overcome this difficulty by combining two ideas.We first introduce a better vector field, Roughly, it allows to replace the bad commutator by one that behaves like νt(1 − k • p)∇ψ, hence vanishing near ±k.Then, we adapt Villani's hypocoercivity method [39], using additional weights in time.The key point of this adaptation is that, besides proving enhanced diffusion at time scale ν −1/2 , see Remark 1.4, it provides extra decay information for some quantities vanishing near the poles ±k of the sphere.This allows to control the commutator and apply the vector field method.All details will be found in Section 4.
Using the decay estimate, the proof of Theorem 2 is achieved in Section 5.1.Regarding the enhanced dissipation estimates of Theorem 3, a key point is to understand the behaviour of the diffusive Volterra kernel past time ν −1/2 , and notably to check the spectral condition (1.31) for small ν.In the Euclidean case, see for instance [5], this can be proved perturbatively from the inviscid case, using the decay of the diffusive kernel uniformly in ν.In our spherical setting, this is where we need the second part of Proposition 1.7.This analysis and the proof of Theorem 3 are performed in Section 5.2.

Volterra equations
In this section we study a general Volterra equation: with unknown u : R + → C n , and data K : We are interested in the global in time solvability of this equation, and in accurate polynomial decay estimates for solution u under assumptions of polynomial decay on K and v.The classical path to study this equation, detailed in [21], is to construct the so-called resolvent of the equation, that is a matrix-valued R : (2.2) Note that the convolution product does not commute when n > 1, so that in this case, the resolvent satisfies two distinct equations.It is then straightforward to check that if K and v are integrable, and if there exists an integrable solution R of (2.2), then an integrable solution u of (2.1) is given by u = v − R ⋆ v and is unique.
For the construction of the resolvent, a vague idea of proof is the following.Assuming there is a solution R to (2.2), and extending both R and K by zero on R − , the relation (2.2) holds now for all t on R, replacing the convolution on R + by the usual convolution on R. Taking the Fourier transform yields in particular This suggests to impose the condition and to define This is however too quick.We remind that, in order to obtain (2.3), we have extended our hypothetical solution R by zero on R − .Hence, we must not only verify that the formula (2.5) makes sense, but also that it gives a function that vanishes on R − .Condition (2.3) is not sufficient for that.Still, under the stronger condition where LK(z) = +∞ 0 e −zt K(t) dt is the Laplace transform of the kernel, everything goes nicely.This is the content of following theorem.THEOREM 4 ([21, Chapter 2]).Assume that K ∈ L 1 (R + , M n (C)), and that the spectral condition (2.6) is satisfied.Then, there exists a unique solution Once global existence of an L 1 solution has been obtained, a natural question is its decay, depending on the decay of the kernel and data.While many works, for instance on Landau damping, treat this kind of question, most results fall into one of the following two types.Either they use weighted spaces and use the stronger assumption [15,20].Or they establish L p in time estimates without loss for p = 1, 2. Indeed, while weighted L 1 estimates behave well with respect to convolution formulas such as (2.2), L 2 estimates may be established using (2.3) and Plancherel formula, see [1].
Note that this is a significant loss of information, and a look at (2.1) shows that we may expect much more: if K, v are O((1+t) −α ) for α > 1, one may hope that u is O((1+t) −α ) as well because this property is stable by convolution.Our main result exactly shows this: THEOREM 5. Assume that (2.6) holds.Let α > 1, and assume that K satisfies for α > 1 and a constant C K .Then, for any v, the solution u of (2.1) satisfies Remark 2.1.Faou, Horsin and Rousset [19, Corollary 3.3] have a related result for specific weights.Our proof is different and easily applies for general weights, notably weight ln(2 + t)(1 + t) −α , that will be used later.
Our proof relies on the analysis of Volterra equations of non-convolution type carried in [21,Chapter 9].These equations take the form where J is a subinterval of R + .An important notion developed there is the one of Volterra kernel of type L p , that we introduce here only for p = ∞.
Definition 2.2.A Volterra kernel on J is a measurable mapping k : J × J → M n (C), such that k(t, s) = 0 for all s, t ∈ J with s > t.
On J we can define a generalized convolution product and one can directly verify that the set of Volterra kernels of type L ∞ on J, equipped with (+, ⋆), is a Banach algebra for the norm |||•||| ∞,J .Moreover, one can show that the space L ∞ (J, C n ) is a left Banach module over it, through (kv As for classical Volterra equations, one has a notion of resolvent: Definition 2.4.Given a Volterra kernel k on J, a resolvent of k on J is another Volterra kernel satisfying As in the convolution case, the resolvent determines the solution.
given by Using the standard von Neumann series for perturbations of the resolvent map in the Banach algebra, we can control the resolvent around a known resolvent.PROPOSITION 2.6 ([21, Chapter 9, Theorem 3.9]).If k = k 1 + k 2 is the sum of two Volterra kernels of type L ∞ on J, if k 1 has a resolvent r 1 of type L ∞ on J, and if then k has a resolvent of type L ∞ on J.
PROOF OF THEOREM 5.For ϵ > 0 (to be chosen later sufficiently small), we consider

By (2.1), they satisfy
with the Volterra kernel By Lemma 2.5 it suffices to show that k is of type L ∞ on R + and has a resolvent of type L ∞ on R + .We decompose k as which can be made arbitrary small by choosing λ small enough.This λ being fixed, we consider now t > λ/ϵ and split the integral at δt for 0 < δ < 1.First, which can be made arbitrary small by choosing some δ < 1 close to 1.Then, for fixed λ, δ we find which again can be made arbitrary small by choosing ϵ sufficiently small.□

Isotropic suspensions in the inviscid case
As the incoherent state is rotational invariant, we can always choose a coordinate system so that k ∈ S 2 equals the coordinate vector e := (0, 0, 1) t .More precisely, if ψ k [ψ in ] is the solution of (1.13) with initial data ψ in , one has for any rotation matrix R: so that it is enough to assume k = e to prove Theorem 1 or Theorem 2. For easier readability we then also drop the explicit dependence in the index k.The system (1.13) reduces to where ).The point is to obtain the decay estimates, first for u, then for ψ itself.We consider in this section the case ν = 0.

Volterra equation on u.
As indicated in Section 1.3, we can rewrite the evolution for u as a Volterra equation We recognize a Fourier transform over the sphere, with well-known decay properties, quantified in Lemma 3.1.Let M ∈ N, δ > 0 and F ∈ H 2M +1+δ (S 2 ) be a function over S 2 .For a unit vector e ∈ S 2 define the integral e i e•p t F (p) dp, t ≥ 0.
Then there exist complex numbers c m,± , 1 ≤ m ≤ M , and and For completeness, we shall provide the proof of this lemma, directly inspired from the lecture notes [8], at the end of this subsection.In the special case of our source U and kernel K, the integrand F of the lemma contains the projection P k ⊥ , that is vanishing at p = ±k so that Lemma 3.1 with M = 1 implies the bounds Hence we find by Theorem 5 that under the spectral condition the first bound (1.15) on the decay of u holds.To complete the proof of Theorem 1, it remains to check when assumption (3.7) is satisfied, and finally to analyse the decay of ψ itself.This will be done in the following Subsections 3.2 and 3.3.
PROOF OF LEMMA 3.1.We introduce a smooth partition of unity χ 0 , χ ± , with χ 0 + χ + + χ − = 1 over S 2 , where χ 0 is supported away from p = ±e and χ ± is supported in a neighborhood of ±e.We decompose By the rotational symmetry, one can introduce coordinates such that e is along the z-axis.Then any point p ∈ S 2 can be parametrized as ).The surface measure on the sphere is dσ = dγ dz so that The key point is that F 0 is compactly supported in (−1, 1) so that its extension by zero to R, still denoted We find in particular that It remains to treat I + (I − can be handled in the same way).Let We find where j is the Jacobian from the change of variable.
The key point is that the phase 1 − |x| 2 has a non-degenerate critical point at x = 0, with Hessian matrix at zero being −2I.By Morse lemma, there exists a smooth diffeomorphism ψ from a neighborhood U of 0 in R 2 to B(0, η), for η small enough, so that 1 − |ψ(y)| 2 = 1 − |y| 2 .By taking the support of χ + sufficiently small, we can perform another change of variables and arrive at where F + is the product of F 0 • φ • ψ with a smooth function, compactly supported in U coming from the Jacobian.Extending F + by zero outside U , we end up with with the Fourier transform F+ of F + and where the last line comes from Plancherel identity.
We can then perform a Taylor expansion Setting we obtain that One can then notice that The result of the lemma follows directly.□

Stability condition.
We now study the stability condition (3.7).This condition was already studied through explicit numerical computations in [36], but here we provide another angle through the argument principle which allows a complete solution.
We first compute the determinant.
Lemma 3.2.For λ ∈ C with Re λ ≥ 0 we have (3.12) PROOF.We remind that We again take the parametrization (3.8).We then deduce from (3.2) and (3.3) for Re λ > 0 that  (2) If ε = −1, the spectral condition is satisfied if and only if Γ < Γ c , where PROOF.For Re λ ≥ 0 define the analytic function F (λ) by By Lemma 3.2, there exists then an eigenmode if and only if F attains in the right half plane the value ε/Γ.
By the explicit expression, we also see that F (λ) → 0 as |λ| → ∞.Hence by the argument principle the attained values are exactly those values encircled by the curve b → F (ib).This curve is plotted in Fig. 2 but can also be understood analytically.
For |b| ≥ 1, the expression directly shows that F (b) is purely imaginary, while for |b| ≤ 1 we find that Re F ≤ 0. Hence it cannot encircle positive real numbers, which proves the first statement of the proposition.
In {|b| < 1} the curve crosses the real axis at b = 0 and b = ±b c .By the expression of F , we see that it indeed crosses the real axis at −3πb We apply Lemma 3.1 with M = 0, to obtain ) is an algebra.As regards I 2 , we write We bound the parenthesis in the first term using again Lemma 3.1 with M = 0: it follows where we have used (1.15) for the second inequality.Eventually, we use again (1.15) to get where we applied the Sobolev imbedding H 1+δ (S 2 ) → L ∞ (S 2 ).This concludes the proof of (1.16) and of Theorem 1.

Decay estimates for advection-diffusion on the sphere
In this central section of our paper we prepare for the results under small angular diffusion ν by proving Proposition 1.6 and Proposition 1.7.
Setting ψ(t) = e tL 1,k ψ in , we thus study decay estimates for solutions to where, from now on, we suppress the subscript on all differential operators, as it is understood that they are all with respect to the variable p on the sphere S 2 .
A first feature of (4.1) is that the interaction of transport and diffusion leads to an enhanced diffusion time scale O(1/ √ ν) that is much shorter than the heat equation one O(1/ν).This phenomenon is analysed in details in Section 4.1.Furthermore, we exhibit an auxiliary time ν −1/3 for the decay of ∥∇(p • e)ψ∥ L 2 , a quantity that vanishes near the pole.This additional time scale is coherent with classical results for enhanced dissipation in the Euclidean setting, where the absence of critical points also leads to typical time ν −1/3 .
In the following Section 4.2, we build a vector field J ν under the form such that, roughly: In the usual Euclidean setting, where e • p is replaced by e • v, v ∈ R 3 , the vector field J = ∇ v + ite is an easy and convenient choice: it commutes with the advection-diffusion operator, so that Jψ can be controlled easily for all times.In the case of the sphere, we do not know how to construct such a commuting vector field.Designing a J ν for which we can show properties (i) and (ii) is difficult, and relies on the refined estimates of Section 4.1.
By these vector fields we can obtain the first part of the mixing estimates of Proposition 1.7 up to time ν −1/2 in Section 4.3.For the estimates on the longer time scale, we adapt the estimates in Section 4.4 which then concludes the proof of Proposition 1.7.

Hypocoercive estimate for advection-diffusion.
In this subsection we will introduce the hypocoercive estimates which will yield the enhanced dissipation of Proposition 1.6.
Using the covariant derivatives ∇ as discussed in Appendix A, we find in our case of S 2 that ∇∆ψ = ∆∇ψ − ∇ψ where the Laplacian ∆ = tr(∇ 2 ) is the connection Laplacian (the correction comes from the Ricci curvature tensor which equals the metric on the sphere).Taking the covariant derivative of (4.1), we therefore get We explicitly compute ∆(∇(p • e)ψ) in Lemma A.1 and this concrete expression yields We further note that |p • e| ≤ 1 and |∇(p • e)| ≤ 1.In spherical coordinates, we find explicitly that |∇(p • e)| = sin θ which is uniformly lower-bounded away from the poles where it vanishes linearly.We now derive several energy identities for the solution to (4.1).In all what follows, ⟨•, •⟩ and ∥ • ∥ stand for the complex scalar product (with the conjugate on the second variable) and norm on L 2 (S 2 ) respectively.For simplicity, we will use the same notations for the scalar product and norm of L 2 vector fields, or more generally for L 2 sections of any tensor bundle over S 2 .A direct computation using the antisymmetry of the operator ip • e allows us to derive the L 2 balance  (4.7) For positive constants a, b, c to be chosen later independently of ν, define the energy functional Such a form of time-dependent functional takes inspiration from [28,14,40].Assuming that 2b 2 ≤ ac, the mixed term can be controlled by the squares so that It follows that (4.11) We will choose the constants according to the following lemma, which is crucial for the subsequent analysis and holds for the longer time-scale ν −1 .
PROOF OF PROPOSITION 1.6.Fix the constants a, b, c according to Lemma 4.1.Let λ > 0 to be fixed later independently of ν.Let ν 0 such that 4aν 0 cλ 2 ≤ 1.From the definition of E[ψ], we have, for all ν ≤ ν 0 , at time where the last line comes from Lemma 4.2, with σ := 4aν cλ 2 ≤ 1.By the evolution estimate of Lemma 4.1 we also find (as ∥ψ∥ is non-increasing), that at time Hence we find that at time t = λν −1/2 , Taking λ small enough (depending on a, b, c but independently of ν), the factor at the right-hand side is less than 1, which implies exponential decay with a rate proportional to ν 1/2 .□ 4.2.Hypocoercive estimate for vector fields.In this subsection we introduce the vector fields to show the mixing estimate Proposition 1.7 for ν > 0 when the semigroup e tL 1 cannot solved explicitly.This has already been used in a few instances for mixing estimates [5,11].In terms of L 1 , the evolution of ψ of (4.1) can be written as and the idea is to find a vector field J for which we control Jψ.A natural candidate for a vector field is J = ∇ + it∇(p • e), which commutes with the inviscid part of the equation.However, it does not commute well with the diffusion operator.Namely, we find that By Duhamel's formula, If we were to rely only on the straightforward (yet optimal on time scales O( 1)) bounds we would get The second term behaves like Cνt 3 , and therefore diverges for t ≫ ν −1/3 , which is a faster time scale than the enhanced dissipation one for this problem.
To overcome this issue, we introduce the viscosity-adapted vector field J ν of the form for scalar functions α = α ν , β = β ν .For the evolution, we then find In order to control the commutator error terms, we set and thus take For the time frame t ≤ ν −1/2 we see that α ∼ 1 and β ∼ t.By this choice, we find that now The gain from the adapted vector field J ν is that the right-hand side now vanishes at the north pole p = e.
As the enhanced dissipation estimate from Section 4.1 provides better decay properties for quantities that vanish at the poles ±e, this will provide a better control of the source term, and in turn a better control of X. Obviously, the right-hand side still does not vanish at the south pole p = −e, so that we need to localize the estimates away from this pole.Symmetrically, one could construct another vector field X for which the roles of the north and south poles would be reversed.
To obtain good control of X (away from the south pole), our starting point is Lemma 4.1 with δ = 1.After integration in time, using that ∥ψ(t)∥ ≤ ∥ψ in ∥ for all t, we get for all t * ≤ ν −1 that sup t∈(0,t * ) (4.19) From there, as a preliminary step, we deduce a few easy bounds on X and R: Lemma 4.3.There exists ν 0 > 0, such that for all 0 < ν ≤ ν 0 and all t * ≤ ν −1 it holds that Furthermore, for and any cutoff χ excluding the south pole −e we have In particular, PROOF.The bounds relative to X follow from the definition X = α∇ψ + iβ∇(p • e)ψ and from the hypocoercive estimate (4.19).For the control of R, we decompose and the estimate also follows from (4.19).□ We now state the key estimates of this paragraph, where we first focus on times up to ν −1/2 .In the first iteration, we use the nested cutoffs χ and χ ′ excluding the south pole −e.Lemma 4.4.There exists ν 0 > 0, such that for all 0 < ν ≤ ν 0 and cutoffs all χ, χ ′ : S 2 → [0, 1] such that χ ′ = 1 on the support of χ and such that −e does not belong to the support of χ ′ , one has for some C > 0 the estimate Before the proof, we show that iterating this estimates shows the estimate with the sharp rates.
Corollary 4.5.There exists ν 0 > 0, such that for all 0 < ν ≤ ν 0 and cutoff χ : S 2 → [0, 1] such that −e does not belong to the support of χ, one has for some C > 0 the estimate In a first step, applying Lemma 4.4 with χ ′ and χ ′′ yields in particular with the factor ν 1/2 the bound where the last inequality comes from Lemma 4.3.Hence we gain the control so that the corollary now follows by another application of Lemma 4.4.□ We now come to the proof of Lemma 4.4.
PROOF OF LEMMA 4.4.The basic idea is to reproduce the hypocoercive estimates of Section 4.1, where we have additional contributions from the cut-off χ and the remainder R. We rewrite (4.17) as We find 1 2 For the last term, we use the explicit expression (4.20) of R and integrate by parts to find It follows from Young's inequality that 1 2 where, for some absolute constant C, For the gradient, we find, using the commutation Hence, by Young's inequality, 1 2 where, for some absolute constant C, It follows by Young's inequality that where Finally, we find using the explicit calculation of Lemma A.1 that For the contribution of R, use the explicit expression of R and the definition of X to find that Hence we find 1 2 where From here, we can proceed exactly as in the proof of Lemma 4.1: for suitable positive constants a, b, c, ν 0 , for all ν ≤ ν 0 , an energy of the form and satisfies for times t ≤ ν −1 the inequality In particular, Regarding E 3 , we have by Young's inequality that for all κ > 0: where still using Lemma 4.3, we have for all In particular, Regarding E 4 , we estimate for all t ≤ ν −1 where, for all Here we used in the last step (4.19) and Lemma 4.3.In particular, this shows For κ small enough independently of ν, and for all t ≤ ν −1 , we get 31) The desired inequality follows by time integration from 0 to ν −1/2 , using the bounds on E 1 , E 2 , Ẽ3 , Ẽ4 .□ To obtain a good decay rate for the velocity field of our linearized model, we will need to control also the L ∞ norm of ψ near the pole.This is the purpose of the next lemma.Lemma 4.6.There exist absolute constants C > 0 and ν 0 > 0 such that for all ν ≤ ν 0 and t ≤ ν −1/2 it holds that Using the Leibniz rule for the covariant derivative, we have where ∆X still refers to the connection Laplacian of the vector field X and where |X| = g(X, X) is the usual norm induced by the metric on tensors.Back to (4.18), it follows that so that where the latter inequality holds for all t ≤ ν −1/2 .It is then easy to find absolute constants c 0 > 0, ν 0 > 0, a 0 , a 1 , a 2 , a 3 > 0, such that for all ν ≤ ν 0 , the function The lemma then follows from the maximum principle for the scalar heat flow on the sphere.□ The estimates on X = J ν ψ can roughly yield a decay like t −1 .For the sharp decay of t −2 for the velocity field, we need to iterate and thus also need a control of J ν X = α∇X + iβ∇(p • e) ⊗ X.From (4.17), we get where we recall R = 2iβν∇([p • e − 1]ψ).By Ricci's formula, see Appendix A, one has [∆, ∇]X = R∇X for some tensor R so that and eventually Using the previous estimates, we shall prove: Lemma 4.7.There exists ν 0 > 0, such that for all 0 < ν ≤ ν 0 and cutoffs χ : S 2 → [0, 1] such that −e does not belong to the support of χ, one can find C > 0 independent of ν satisfying for all t ≤ ν −1/2 ∥J ν X(t)χ∥ 2 ≤ C∥ψ in ∥ 2 H 2 .PROOF.Let χ, χ ′ such that χ ′ = 1 on the support of χ and such that −e does not belong to the support of χ ′ .Performing an L 2 estimate on (4.33), we find after standard integrations by parts that 1 2 We first notice that For the last term, we use that For the first part, we find by integration by parts that This right-hand side was already encountered above and we bound it as For the other component, we further compute that

One has
Re⟨R PROOF OF PROPOSITION 1.7 FOR TIMES UP TO ν −1/2 .We start with the proof of the first inequality.We split the integrand where χ is a smooth function which is 1 near e and 0 near −e.By symmetry consideration, it is enough to show the decay for S 2 ψ(t)F χ.We introduce a cutoff χ ϵ which is 1 in a ball of radius ϵ around the pole p = e, zero outside a ball of radius 2ε, and satisfies |∇χ ϵ | ≲ ϵ −1 .We write As regards I 2 , we write We introduce again spherical coordinates (θ, φ), with colatitude θ ∈ (0, π) and longitude φ ∈ (0, 2π), so that p = sin θ cos φ e x + sin θ sin φ e y + cos θ e, and ∇(p • e) = − sin θe θ , while the surface measure on the sphere is ds = sin θdθdφ.For the first term, we find that In particular with Corollary 4.5 this shows for times For the second term, we get In particular for times t ≤ ν −1/2 it shows that Taking ϵ = 1 √ t yields the first inequality.As regards the second inequality, we write Applying the first inequality to the integral in J 2 , we find for times t ≤ ν −1/2 that The term J 1 is similar, except that ψ is replaced by X, on which we have a weaker control.Using the same strategy as in the proof of the first inequality, we can write We treat K 1 as I 1 , except that we use Lemma 4.6 as a substitute to the L ∞ bound on ψ.We get We then treat K 2 as I 2 , resulting in As before we find For K 2,2 , we further decompose it as For H 1 we obtain the bound and in particular for times t ≤ ν −1/2 we get by Lemma 4.6 that For H 2 , we compute In particular for t ≤ ν −1/2 it shows by Lemma 4.7 that In particular, for t ≤ ν −1/2 we use Lemma 4.6 to find Taking ϵ = 1 √ t yields the second inequality.□ 4.4.Mixing estimates after time ν −1/2 .In this subsection, we prove the remaining estimates in Proposition 1.7, i.e. for times t ∈ [ν −1/2 , cν −1/2 | ln ν|].The general idea is the same and it has already been anticipated in the previous subsections.The main difference is to notice that we cannot approximate anymore α by 1 and β by t but that α and β are exponentially growing in this time scale.PROOF OF PROPOSITION 1.7 FOR TIMES AFTER ν −1/2 .Performing the proof of Lemma 4.4, we do not substitute α and β in the last step but rather take out a supremum.As an analogue of Corollary 4.5 this yields for all t * ≤ ν −1 the bound Regarding the L ∞ bound on X near the poles, we restart from (4.32) to deduce where we used that |β| 2 /(|α| 2 t 2 ) ≲ 1 for the second inequality.Mimicking the end of the proof of Lemma 4.6, we find that for all t * ≤ cν −1/2 | ln ν|, and using (4.39)-(4.40)-(4.41),we can adapt the proof of Section 4.3 to obtain the second part of Proposition 1.7.□

Proof of mixing and enhanced diffusion for the diffusive Saintillan-Shelley model
The good decay estimates for the semigroup e tL 1 of the previous subsection now allow us to conclude the results by studying the Volterra equation (1.28).
5.1.Proof of Theorem 2. We first prove Theorem 2. Coming back to the Volterra equation (1.28), we drop the subscript k and restore the superscript ν so that it now reads with We remind that Thanks to the identity we find that Then, Proposition 1.7 shows that there exists a constant C independent of ν such that and consider the corresponding modified Volterra equation Let K inv be the inviscid kernel considered in Section 3, that is the kernel for ν = 0. We have seen that K inv decays like t −2 .We now claim that Indeed, let κ > 0 be arbitrarily small.We fix some large T > 0 independent of ν, so that Note that the second condition can be achieved because the constant C in (5.3) is independent of ν.Now, for all ν such that ν −1/2 > T , we write Eventually, it is standard to show that, given any initial data in L 2 (S 2 ), and any finite time interval (0, T ), the solution ψ ν of the advection-diffusion equation )) to the solution ψ inv of ∂ t ψ inv + i(p • k)ψ inv = 0 (with the same initial data).
This implies in particular that det(I + LK ν ) cannot vanish in the unstable half plane.
Acting on a general tensor X, we find a similar expression with a successive application of the Riemann curvature tensor.As the Riemann curvature tensor is bounded on a sphere, we just note that [∆, ∇]X = R∇X for a bounded tensor R.
In the enhanced dissipation estimates, we use the commutator between the Laplacian and the tensor (p • e)X, for a fixed vector e.It is provided by PROOF.We compute the expression in the spherical coordinates (θ, φ) where we take e to be the north pole.Using the Leibniz rule, we find that ∆(∇ i 0 (p • e)X i 1 ,...,in ) = (∆∇ i 0 (p • e))X i 1 ,...,in ) + 2g mn (∇ m ∇ i 0 (p • e))∇ n X i 1 ,...,in + ∇ i 0 (p • e)∆X i 1 ,...,in .
Hence we have arrived at the claimed expression.□ Conflict of interest.The authors declare that they have no conflict of interest.

FF/ 2 F/ 2
Pusher: bacteria with flagella at the back, e.g.E. coli.Induced flow (B) Puller: bacteria with flagella at the side at the front, e.g. C. reinhardtii.

. 13 )PROPOSITION 3 . 4 .
This is an odd function, that vanishes at 0 and at ±b c where b c > 0 can be evaluated numerically to b c ≈ 0.62375.Moreover, the function is negative on (0, b c ), positive on (b c , 1).For the inviscid system, one has depending on ε(1) If ε = 1, then the spectral condition (3.7) is satisfied for all Γ ∈ R + .