# Hydrodynamic limit of multiscale viscoelastic models for rigid particle
suspensions

Mitia Duerinckx, Lucas Ertzbischoff,
Alexandre Girodroux-Lavigne, Richard M. Höfer  Université Libre de Bruxelles,
Département de Mathématique, 1050 Brussels, Belgium<EMAIL_ADDRESS>Department of Mathematics, Imperial College London, London, SW7 2AZ, United-
Kingdom<EMAIL_ADDRESS>University of Regensburg, Faculty of
Mathematics, Regensburg, Germany<EMAIL_ADDRESS>regensburg.de University of Regensburg, Faculty of Mathematics, Regensburg,
Germany<EMAIL_ADDRESS>
###### Abstract.

We study the multiscale viscoelastic Doi model for suspensions of Brownian
rigid rod-like particles, as well as its generalization by Saintillan and
Shelley for self-propelled particles. We consider the regime of a small
Weissenberg number, which corresponds to a fast rotational diffusion compared
to the fluid velocity gradient, and we analyze the resulting hydrodynamic
approximation. More precisely, we show the asymptotic validity of macroscopic
nonlinear viscoelastic models, in form of so-called ordered fluid models, as
an expansion in the Weissenberg number. The result holds for zero Reynolds
number in 3D and for arbitrary Reynolds number in 2D. Along the way, we
establish several new well-posedness and regularity results for nonlinear
fluid models, which may be of independent interest.

###### Contents

  1. 1 Introduction
  2. 2 Doi–Saintillan–Shelley theory
  3. 3 Ordered fluid models
  4. 4 Statement of main results
  5. 5 Small-$\operatorname{Wi}$ expansion of Doi–Saintillan–Shelley theory
  6. A Well-posedness of the Doi–Saintillan–Shelley system
  7. B Perturbative well-posedness of ordered fluid equations
  8. C Derivation of third-order fluid equations

##  1\. Introduction

###  1.1. General overview

Suspensions of rigid particles in Stokesian fluids are ubiquitous both in
nature and in applications and are known to typically exhibit non-Newtonian
behaviors. These systems can be described on different scales: either by
macroscopic non-Newtonian fluid models, or by so-called multiscale kinetic
models, or else on the microscopic scale as suspended particles moving in a
fluid flow. Let us briefly describe these different levels of physical
description:

  1. $\bullet$

Macroscopic non-Newtonian fluid models:
An incompressible fluid flow is generally described by the Navier–Stokes
equations,

$\left\\{\begin{array}[]{l}\rho_{\operatorname{fl}}(\partial_{t}+u\cdot\nabla)u-{\operatorname{div}}(\sigma)+\nabla
p\,=\,h,\\\\[2.84526pt] {\operatorname{div}}(u)=0,\end{array}\right.$ (1.1)

where $u:\mathbb{R}^{+}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ is the fluid
velocity field, where the source term $h$ accounts for internal forces, where
$\rho_{\operatorname{fl}}$ stands for the fluid density (constant, say), and
where $\sigma$ is the deviatoric fluid stress (a trace-free symmetric matrix
field). For Newtonian fluids, the stress is linear in the strain rate, that
is,

$\sigma=2\mu\operatorname{D}(u),$ (1.2)

where the strain rate $\operatorname{D}(u)=\frac{1}{2}(\nabla u+(\nabla
u)^{T})$ is defined as the symmetric gradient of the fluid velocity, and where
the viscosity $\mu\geq 0$ is a constant. In contrast, non-Newtonian fluids are
characterized by more complex constitutive laws relating the stress $\sigma$
to the strain rate $\operatorname{D}(u)$, describing various possible types of
nonlinear and memory effects. In _macroscopic_ non-Newtonian models, these
laws are assumed to take some explicit form and are usually fitted
phenomenologically to experimental rheological measurements. In Reiner–Rivlin
and in generalized Newtonian fluid models, the stress $\sigma$ is simply taken
to be a local function of $\operatorname{D}(u)$, meaning that only nonlinear
effects are retained. To further describe memory effects of non-Newtonian
fluids, such as viscoelastic properties, more realistic models rather relate
$\sigma$ and $\operatorname{D}(u)$ via integral or differential equations,
which aim to take into account the dependence of the stress on the fluid
deformation history. Such models that are frame invariant are generically
called _simple fluids_ , of which the celebrated Oldroyd–B and FENE–P models
are particular cases. In the fast relaxation limit, simple fluid models reduce
to the hierarchy of so-called ordered fluid models. We refer to Section 3 for
details.

  2. $\bullet$

Multiscale kinetic models for suspensions:
So-called multiscale or micro-macro models go one step away from the pure
phenomenology, towards a microscopically more accurate description of particle
suspensions. More precisely, the macroscopic fluid equation (1.1) gets coupled
via the stress $\sigma$ to a kinetic equation describing the evolution of the
suspended solid phase. The latter is modeled by a particle density function
$f(t,x,n)\in\mathbb{R}^{+}$ at time $t$, where $x$ is the position of
particles and $n$ is their ‘state’: for instance, $n\in\mathbb{R}^{d}$ may
describe the relative position of endpoints of elongated particles (hence
$n\in\mathbb{S}^{d-1}$ in case of rigid suspended particles as will be
considered in this work). The evolution of the particle density $f$ is then
modeled by a Fokker–Planck equation describing transport with the fluid and
diffusive effects. Finally, the coupling to the macroscopic fluid equation is
expressed through an explicit constitutive law
$\sigma=\sigma(f,\operatorname{D}(u))$, which is typically derived from formal
microscopic considerations in dilute regime. Such models describe how the
microscopic state of the particles adapts collectively to local fluid
deformations and how the macroscopic fluid flow gets itself effectively
impacted. Popular models include the kinetic FENE and Hookean dumbbell models
for dilute suspensions of flexible polymers, the so-called Doi model for
suspensions of Brownian rigid rod-like particles, and the
Doi–Saintillan–Shelley for corresponding active particles. We refer to Section
2 for details.

  3. $\bullet$

Microscopic models:
At the particle scale, we can formulate a fully detailed hydrodynamic model
describing the motion of suspended particles in the Stokesian background fluid
flow. It takes form of equation (1.1) restricted to the fluid domain, with
Newtonian constitutive law (1.2), coupled with Newton’s equations of motion
for the particles. We refer for instance to [HS23] for a discussion of this
complex dynamics.

From the modeling perspective, while microscopic models are certainly
impractical due to the huge number of particles in real-life systems, one can
argue that multiscale kinetic models are more satisfactory than macroscopic
fluid models as they retain some information from the fluid-particle coupling
and therefore better reveal mechanisms leading to non-Newtonian behavior. Yet,
macroscopic models appeal through their simpler description: in particular,
they are much more accessible for numerical simulations and proner to
comparison with experimental rheological measurements.

Most previous works on particle suspensions have aimed either to study
properties of macroscopic non-Newtonian fluid models or multiscale kinetic
models, or else to derive multiscale kinetic models rigorously from
microscopic particle dynamics, which has indeed attracted considerable
interest in recent years; see Section 1.3 below for references. In the present
contribution, we propose to fill the gap in the micro-macro understanding of
non-Newtonian effects of particle suspensions by further studying the
derivation of explicit macroscopic fluid models from multiscale kinetic models
in suitable regimes.

In some exceptional cases, the kinetic equation describing the state of the
particles in micro-macro models can be integrated out and leads to a closed
equation for the stress $\sigma$ in terms of the strain rate
$\operatorname{D}(u)$. This is for instance the case for the kinetic Hookean
dumbbell model, which is well-known to be formally equivalent to the Oldroyd–B
macroscopic fluid model (see e.g. the recent rigorous analysis in [DS23]). In
general, however, no exact macroscopic closure is available and we can only
hope for perturbative closures to be valid in suitable asymptotic regimes.

In fact, in the regime of weak non-Newtonian effects, a fairly large class of
non-Newtonian macroscopic fluid models is believed to be well approximated by
a special family of models, called _ordered fluid_ models. More precisely, the
latter are expected to be good approximations for all viscoelastic fluids in
the fading memory regime, that is, in the regime when the elastic time-
dependent effects due to suspended particles in the fluid have an inherent
relaxation timescale that is much shorter than the overall timescale of the
fluid flow. This ratio of timescales is the so-called Weissenberg number
$\operatorname{Wi}$. In other words, ordered fluid models arise formally as
expansions of any viscoelastic fluid model at small $\operatorname{Wi}$, which
is sometimes referred to as the retarded motion expansion. A first-order fluid
is a Newtonian fluid, a second-order fluid is a non-Newtonian fluid where
effects of order $O(\operatorname{Wi}^{2})$ are neglected, a third-order fluid
amounts to neglecting effects of order $O(\operatorname{Wi}^{3})$, etc.

In the present work, we focus on the so-called Doi model, which is a
multiscale kinetic model describing suspensions of Brownian rigid rod-like
particles, and we further consider its generalization by Saintillan and
Shelley for active (self-propelled) particles. We rigorously analyze the
hydrodynamic limit of these models in the small-$\operatorname{Wi}$ expansion,
which was extensively studied on a formal level in the physics literature in
the 1970s, see [HL72, Bre74], and we confirm the asymptotic validity of
ordered fluid models in this setting. As a prerequisite, the justification of
the asymptotic expansion requires a careful study both of the
Doi–Saintillan–Shelley model and of ordered fluid models: in particular, we
establish some new well-posedness and regularity results for these nonlinear
viscoelastic fluid models, which we believe are nontrivial and of independent
interest. Moreover, as the Doi–Saintillan–Shelley model can itself be derived
from a microscopic hydrodynamic model (at least formally), see Section 1.3,
our derivation of macroscopic ordered fluid models comes together with
explicit expressions for rheological parameters in terms of microscopic
characteristics of the underlying particle suspension.

###  1.2. Informal statement of main results

We start from the following dimensionless Doi–Saintillan–Shelley model
describing suspensions of (active or passive) Brownian rigid rod-like
particles in a Stokesian fluid,

$\left\\{\begin{array}[]{l}\operatorname{Re}\big{(}\partial_{t}u_{\varepsilon}+(u_{\varepsilon}\cdot\nabla)u_{\varepsilon}\big{)}-\Delta
u_{\varepsilon}+\nabla
p_{\varepsilon}=h+\frac{1}{\varepsilon}{\operatorname{div}}(\sigma_{1}[f_{\varepsilon}])+{\operatorname{div}}(\sigma_{2}[f_{\varepsilon},\nabla
u_{\varepsilon}]),\\\\[2.84526pt]
\partial_{t}f_{\varepsilon}+{\operatorname{div}}_{x}\big{(}(u_{\varepsilon}+U_{0}n)f_{\varepsilon}\big{)}+{\operatorname{div}}_{n}\big{(}\pi_{n}^{\bot}(\nabla
u_{\varepsilon})nf_{\varepsilon}\big{)}=\tfrac{1}{\operatorname{Pe}}\Delta_{x}f_{\varepsilon}+\tfrac{1}{\varepsilon}\Delta_{n}f_{\varepsilon},\\\\[2.84526pt]
{\operatorname{div}}(u_{\varepsilon})=0,\end{array}\right.$ (1.3)

in terms of the elastic and viscous stresses

$\displaystyle\sigma_{1}[f]$ $\displaystyle:=$
$\displaystyle\lambda\theta\int_{\mathbb{S}^{d-1}}\big{(}n\otimes
n-\tfrac{1}{d}\operatorname{Id}\big{)}f(\cdot,n)\,dn,$ (1.4)
$\displaystyle\sigma_{2}[f,\nabla u]$ $\displaystyle:=$
$\displaystyle\lambda\int_{\mathbb{S}^{d-1}}(n\otimes n)(\nabla u)(n\otimes
n)\,f\,dn,$ (1.5)

where $\varepsilon:=\operatorname{Wi}>0$ stands for the Weissenberg number,
$\operatorname{Re}\geq 0$ is the Reynolds number, $\operatorname{Pe}>0$ is the
so-called Péclet number, and where the source term $h$ is taken to be smooth
and accounts for internal forces. The state variable $n$ in the kinetic
Fokker–Planck equation describes the orientation of rigid particles on the
unit sphere, hence $n\in\mathbb{S}^{d-1}$: we write $\Delta_{n}$ and
${\operatorname{div}}_{n}$ for the Laplace–Beltrami operator and the
divergence on the sphere $\mathbb{S}^{d-1}$, and we also use the short-hand
notation $\pi_{n}^{\bot}:=\operatorname{Id}-n\otimes n$ for the orthogonal
projection onto $n^{\bot}$. The above system is introduced in detail in
Section 2, where in particular the constants
$U_{0},\lambda,\theta\in\mathbb{R}$ are further described. The Doi model for
passive suspensions is recovered for the special choice $U_{0}=0$ and
$\theta=6$. We set this model for simplicity in a finite box
$\mathbb{T}^{d}=[0,1)^{d}$ with periodic boundary conditions, and consider
space dimension $d=2$ or $3$.

Note that we take into account a non-vanishing spatial diffusion
$O(\frac{1}{\operatorname{Pe}})$ in (1.3), which differs from the setting
usually considered in applications where $\operatorname{Pe}=\infty$: a
nontrivial spatial diffusion $\operatorname{Pe}<\infty$ is actually needed in
the present work for technical well-posedness reasons. Due to this spatial
diffusion, the structure of ordered fluid equations needs to be slightly
adapted, leading to the nonstandard definition of Rivlin–Ericksen tensors in
(1.13) below; see Section 3 for details. This is reminiscent of the version of
the Oldroyd–B model with stress diffusion that is often considered both for
analytical and numerical studies [RT21].

Our main result is the following asymptotic validity of second-order fluid
models. We refer to Theorem 4.2 in Section 4 for a more detailed statement,
including the precise well-preparedness requirement for kinetic initial data,
as well as the explicit expression for the effective coefficients
$\mu,\nu,\gamma_{1},\gamma_{2}$ and their rheological interpretation. New
results on the well-posedness of the kinetic model (1.3) and of the second-
order fluid model (1.13) are postponed to Sections 2 and 3, respectively.

######  Theorem 1.1 (Informal statement of the main result).

Consider either the Stokes case $\operatorname{Re}=0$ with $d\leq 3$, or the
Navier–Stokes case $\operatorname{Re}=1$ with $d=2$. Given an initial particle
density $f_{\varepsilon}^{\circ}\in
C^{\infty}\cap\mathcal{P}(\mathbb{T}^{d}\times\mathbb{S}^{d-1})$ that is well-
prepared in a sense that will be clarified later (see Assumption 4.1), and
given also an initial fluid velocity $u^{\circ}\in
C^{\infty}(\mathbb{T}^{d})^{d}$ with ${\operatorname{div}}(u^{\circ})=0$ in
the Navier–Stokes case, consider a weak global solution
$(u_{\varepsilon},f_{\varepsilon})$ of the Cauchy problem for the
Doi–Saintillan–Shelley model (1.3). For all $T>0$, provided that

$\varepsilon\ll 1\qquad\text{and}\qquad\lambda\theta(1+\operatorname{Pe})\ll
1$ (1.6)

are small enough, the fluid velocity $u_{\varepsilon}$ satisfies

$\displaystyle\begin{array}[]{rlll}\|\nabla(u_{\varepsilon}-\bar{u}_{\varepsilon})\|_{\operatorname{L}^{2}(0,T;\operatorname{L}^{2}(\mathbb{T}^{d})^{d^{2}})}&\lesssim&\varepsilon^{2},&\quad\text{if\leavevmode\nobreak\
\leavevmode\nobreak\ $\operatorname{Re}=0$},\\\\[2.84526pt]
\|u_{\varepsilon}-\bar{u}_{\varepsilon}\|_{\operatorname{L}^{\infty}(0,T;\operatorname{L}^{2}(\mathbb{T}^{d})^{d})}+\|\nabla(u_{\varepsilon}-\bar{u}_{\varepsilon})\|_{\operatorname{L}^{2}(0,T;\operatorname{L}^{2}(\mathbb{T}^{d})^{\small{d^{2}}})}&\lesssim&\varepsilon^{2},&\quad\text{if\leavevmode\nobreak\
\leavevmode\nobreak\ $\operatorname{Re}=1$},\end{array}$

and the particles’ spatial density
$\rho_{\varepsilon}:=\fint_{\mathbb{S}^{d-1}}f_{\varepsilon}(\cdot,n)\,\mathrm{d}n$
satisfies

$\|\rho_{\varepsilon}-\bar{\rho}_{\varepsilon}\|_{\operatorname{L}^{\infty}(0,T;\operatorname{L}^{2}(\mathbb{T}^{d}))}+\|\nabla(\rho_{\varepsilon}-\bar{\rho}_{\varepsilon})\|_{\operatorname{L}^{2}(0,T;\operatorname{L}^{2}(\mathbb{T}^{d})^{d})}\,\lesssim\,\varepsilon^{2},$

where $(\bar{u}_{\varepsilon},\bar{\rho}_{\varepsilon})$ solves the following
(non-standard) second-order fluid equation

$\displaystyle\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+\bar{u}_{\varepsilon}\cdot\nabla)\bar{u}_{\varepsilon}-{\operatorname{div}}(\bar{\sigma}_{\varepsilon})+\nabla\bar{p}_{\varepsilon}\,=\,h+O(\varepsilon^{2}),\\\\[2.84526pt]
(\partial_{t}+\bar{u}_{\varepsilon}\cdot\nabla)\bar{\rho}_{\varepsilon}-(\tfrac{1}{\operatorname{Pe}}+\varepsilon\nu)\Delta\bar{\rho}_{\varepsilon}=O(\varepsilon^{2}),\\\\[2.84526pt]
\bar{\sigma}_{\varepsilon}=(1+\mu\bar{\rho}_{\varepsilon})A_{1}(\bar{u}_{\varepsilon})+\varepsilon\gamma_{1}A_{2}^{\prime}(\bar{u}_{\varepsilon},\bar{\rho}_{\varepsilon})+\varepsilon\gamma_{2}\bar{\rho}_{\varepsilon}A_{1}(\bar{u}_{\varepsilon})^{2},\\\\[2.84526pt]
{\operatorname{div}}(\bar{u}_{\varepsilon})=0,\end{array}\right.$ (1.12)

for some explicit coefficients $\mu,\nu,\gamma_{1},\gamma_{2}\in\mathbb{R}$,
where $A_{1}(u):=2\operatorname{D}(u)=(\nabla u)^{T}+\nabla u$ is the strain
rate and where $A_{2}^{\prime}$ is a (non-standard) inhomogeneous, diffusive
version of the second Rivlin-Ericksen tensor, defined by

$A_{2}^{\prime}(\rho,u)\,:=\,\big{(}\partial_{t}-{\tfrac{1}{\operatorname{Pe}}}\Delta+u\cdot\nabla\big{)}\big{(}\rho
A_{1}(u)\big{)}+\rho\big{(}(\nabla u)^{T}A_{1}(u)+A_{1}(u)(\nabla
u)\big{)}.\qed$ (1.13)

######  Remark 1.2.

A few comments are in order:

  1. (a)

The well-preparedness assumption for the initial kinetic data will be
clarified later on in Assumption 4.1. Informally, it amounts to assuming that
initial data are locally at equilibrium with respect to the dynamics of
orientations and are perturbatively compatible with the formal
$\varepsilon$-expansion. It allows to avoid initial boundary layers.

  2. (b)

The existence of global weak solutions for the Doi–Saintillan–Shelley model
(1.3), which is assumed to hold in the above statement, is indeed proven in
Section 2.2 below. We further establish a new weak-strong uniqueness principle
for this system, which implies some stability that is at the very heart of the
above result.

  3. (c)

Although the second-order fluid model (1.12) is well-known to be ill-posed
whenever $\gamma_{1}<0$ (which is indeed the case for relevant effective
parameters), we can define several well-posed notions of approximate solutions
that only satisfy (1.12) up to higher-order $O(\varepsilon^{2})$ errors. This
discussion is postponed to Section 3, where we present two approaches to fix
this issue: First, we introduce the notion of approximate hierarchical
solutions, which naturally appear as low-$\operatorname{Wi}$ expansions; see
Propositions 3.2 and 3.4. Second, by means of a Boussinesq-type perturbative
rearrangement, we additionally provide a reformulation of second-order fluid
equations in terms of a closed well-posed system; see Proposition 3.3.

  4. (d)

The explicit expression for the effective second-order fluid coefficients
$\mu,\nu,\gamma_{1},\gamma_{2}$ is postponed to Section 4, where we further
describe the rheological properties of the obtained macroscopic fluid model.
The expressions for the coefficients agree with those computed in [HL72,
Bre74] in the case of passive suspensions. Moreover, they qualitatively match
experimental data and formal predictions on active suspensions.

  5. (e)

A similar result could be obtained with the same approach when starting from
the co-rotational kinetic FENE model for elastic polymers (see e.g. [LM07] for
a review of this system). For conciseness, we do not repeat our analysis in
that setting and leave the adaptation to the reader. ∎

###  1.3. Previous results

We briefly review previous rigorous results related to the multiscale
description of particle suspensions and related systems.

#### Derivation of the Doi–Saintillan–Shelley model

The systematic theoretical study of the effective rheology of suspensions has
been initiated by Einstein [Ein06], who found that passive non-Brownian
spherical rigid particles effectively increase the fluid viscosity by
$\frac{5}{2}\phi\mu$, where $\phi$ is the volume fraction of the particles and
$\mu$ is the viscosity of the solvent. Jeffery [Jef22] studied the analogous
problem for ellipsoidal particles and found an increase of the viscosity
depending locally on orientations of the particles. By slender-body theory, in
the limit of very elongated particles, Jeffery’s viscous stress exactly takes
the form of $\sigma_{2}$ in (1.3); see e.g. [Bre74, KK13]. Starting in the
1930s, there is a vast literature in physics on the rheology of suspensions of
non-spherical rigid particles, see e.g. [Kuh32, Eis33, Pet38, Bur38], but this
early work was restricted to specific fluid flows like simple shear, and
Brownian effects were neglected. Brownian particles were first considered in
[Sim40, KK45, RK50, Sai51], by means of different models and justifications,
finally leading to the additional elastic stress $\sigma_{1}$ in (1.3). These
models were largely reviewed in [LH71, HL72, Bre74]: in particular, the
multiscale kinetic model (1.3) in the passive case ($U_{0}=0$, $\theta=6$)
then entered textbooks such as [DE88, Gra18] and became known in the
mathematical literature as the Doi model. The extension to active suspensions
has been proposed by Saintillan and Shelley [SS08]: by coarse-graining force
dipoles exerted by the particles on the fluid, they derived a further
contribution of the elastic stress $\sigma_{1}$ due to particles’ activity;
see also [HABK08, Hai+09, Sai10, PRB16, DVY19].

On the mathematical side, the derivation of the viscous stress $\sigma_{2}$
from microscopic models has received considerable attention in the last years.
When the fluid is modeled by the Stokes equations and the particle
distribution is given, the effective increase of the fluid viscosity is by now
well understood [HM12, NS20, GVH20, HW20, DG23, GV21, GVH21, DG21, GVM22,
Due22, DG23a]; see [DG22] for a review. In a similar setting, for active
suspensions, the elastic stress $\sigma_{1}$ has been derived in [GL22,
BDG22]. Yet, non-Newtonian effects originate from the retroaction of the fluid
on the particles, that is, from coupling the particle density to the fluid
flow: beyond the derivation of $\sigma_{1}$ and $\sigma_{2}$ in the static
setting, the derivation of kinetic equations for the particle density in the
time-dependent setting is of key interest. First steps in that direction have
been undertaken in [HS21, HMS24, Due23]. Finally, regarding the derivation of
the passive part of the elastic stress $\sigma_{1}$ for Brownian particles, we
refer to [HLM23], where the authors start from a simplified microscopic model
where the particle dynamics is given by Brownian motion and not coupled to the
fluid.

In a different direction, we also mention recent work [HS10, AO22, CZDGV22],
where the (linear) stability and mixing properties of the
Doi–Saintillan–Shelley model have been investigated (neglecting however the
viscous stress $\sigma_{2}$).

#### Macroscopic rheology and formal closures

Most of the above-mentioned works in the physics literature on the derivation
of the Doi model do not stop at multiscale models but also aim at macroscopic
non-Newtonian models, as well as at explicit calculations of stresses in
specific flows such as simple shear. In particular, they typically give
formulas for the shear viscosity in simple shear flow for very large or very
small Weissenberg number $\operatorname{Wi}$. Normal-stress differences (see
(3.6) below) have also been computed at small $\operatorname{Wi}$ by Giesekus
[Gie62], showing that the elastic stress $\sigma_{1}$ does not contribute to
the second normal-stress difference but that the viscous stress does. This was
in contradiction with Weissenberg’s original conjecture that all real-world
fluids must have vanishing second normal-stress difference, a conjecture that
was later falsified also through experiments. Hinch and Leal [LH71, HL72]
systematically computed expansions for the stress in simple shear flow both at
small and at large $\operatorname{Wi}$. They also noticed that their findings
at small $\operatorname{Wi}$ agree with second-order fluid models, but they
did not investigate whether this holds in more general fluid flows.

Although no exact macroscopic closure is available for the Doi model in the
general non-perturbative setting, formal approximate closures have been
studied for example in [DE88, Chapter 8.7]. A particular instance is the
Oldoyd–B model, which is well-known to be remarkably an exact macroscopic
closure for the kinetic Hookean dumbbell model; this formal connection was
made rigorous in [DS23] (in the stress-diffusive case). For the Doi and
Doi–Saintillan–Shelley models, an exact closure is not available, thus raising
the question whether macroscopic closures are at least asymptotically valid in
some scaling regimes: for instance, we will see that the validity of the
Oldoyd–B closure already fails at order $O(\operatorname{Wi})$ for the Doi
model in the small-$\operatorname{Wi}$ regime (see Remark 3.1). The situation
is very similar for kinetic FENE models for elastic polymers: no exact
macroscopic closure holds, but the so-called FENE–P model is still a very
popular formal approximate closure (see e.g. [BDJ80], where the name FENE–P is
attributed in reference to an earlier paper of Peterlin [Pet66]).

#### Small-$\operatorname{Wi}$ regime and hydrodynamic limits

At small $\operatorname{Wi}$, the Doi–Saintillan–Shelley system (1.3)
undergoes a strong diffusion in orientation (cf. factor
$\frac{1}{\varepsilon}=\frac{1}{\operatorname{Wi}}$ in front of $\Delta_{n}$):
to leading order, the orientations of the particles relax instantaneously to
the steady state, which corresponds to isotropic orientations. While this
leading order amounts to a trivial Newtonian behavior, next-order
$O(\operatorname{Wi})$ corrections encode nontrivial non-Newtonian effects
where the stress starts to depend on the local fluid deformation and its
history. The formal perturbative expansion in powers of $\operatorname{Wi}$ is
comparable to the Hilbert expansion method in the Boltzmann theory [Hil12,
Caf80, Gol05, SR09], where one looks for a solution of the Boltzmann equation
as a formal power series in terms of the Knudsen number $\mathrm{Kn}\ll 1$ and
where the leading-order approximation simply leads to compressible Euler
equations. From the closely related Chapman–Enskog asymptotics expansion, one
can (formally) obtain compressible Navier–Stokes equations as a
$O(\mathrm{Kn})$ correction to the compressible Euler system, up to
$O(\mathrm{Kn}^{2})$ errors (see e.g. [Gol05, Section 5.2] or [SR09, Section
2.2.2]). In a similar way, for the Doi model, second-order fluid equations are
obtained in this work as a $O(\operatorname{Wi})$ correction to the Stokes
equations, up to $O(\operatorname{Wi}^{2})$ errors.

Our results can be compared with corresponding results for the Doi–Onsager
model for liquid crystals, which indeed shares some similarities with the Doi
model that we consider in this work. In [EZ06, WZZ15], the macroscopic
Ericksen–Leslie system is derived from the Doi–Onsager model by means of a
Hilbert expansion. Note however that in that case the leading term in the
expansion already yields a non-trivial system, so that higher-order
corrections are not investigated in [EZ06, WZZ15]. We also mention recent
related work on hydrodynamic limits for alignment models [DM08, DFMA17,
DFMAT19, DFL22] and for flocking models [KMT15, KV15, FK18], as well as a
preliminary result for kinetic FENE models for elastic polymers [LPD02].

Note that hydrodynamic limits have been investigated also for various other
kinetic models for particle suspensions in different settings: for instance,
in the context of sedimentation for small inertial particles, let us mention
the inertialess limits studied in [Jab00, Höf18, HKM23, Ert23], as well as the
high-diffusion limit in velocity investigated in [GJV04, GJV04a, MV08, SY20].

###  1.4. Structure of the article

The article is split into five main sections, in addition to three appendices
containing proofs of secondary results:

  1. $\bullet$

In Section 2, we give a detailed account of the Doi–Saintillan–Shelley model:
we describe the non-dimensionalization leading to its form (1.3), and we state
a new well-posedness result for this system, see Proposition 2.1, which we
prove in Appendix A.

  2. $\bullet$

Section 3 provides an introduction to ordered fluid models, starting with
their basic non-Newtonian rheological properties. We then discuss the ill-
posedness of second-order fluid equations for the relevant range of
coefficients, and we present two approaches to fix this issue: we introduce
hierarchical solutions in Proposition 3.2 and we define a Boussinesq-type
perturbative rearrangement in Proposition 3.3. As we allow the suspended
particle density to be spatially inhomogeneous in general, and as we include a
non-vanishing spatial diffusion $O(\frac{1}{\operatorname{Pe}})$ in the model
for technical reasons, we actually derive ordered fluid models of the form
(1.12), which slightly differ from their standard version: we motivate and
introduce these nonstandard models in Section 3.4 and give the corresponding
definition of hierarchical solutions in Proposition 3.4. The proofs of
Propositions 3.2, 3.3 and 3.4 are postponed to Appendix B.

  3. $\bullet$

In Section 4, we provide a more detailed formulation of our main result on the
derivation of second-ordered fluid equations from the Doi–Saintillan–Shelley
theory at small Weissenberg number, see Theorem 4.2. We further comment on the
derivation of higher-order fluid models, for which details are postponed to
Appendix C.

  4. $\bullet$

In Section 5 we prove our main result, that is, the rigorous
$\varepsilon$-expansion of the solution $(u_{\varepsilon},f_{\varepsilon})$ of
the Doi–Saintillan–Shelley system (1.3).

### Notations

We summarize the main notations that we use in this work:

  1. $\bullet$

We denote by $C\geq 1$ any constant that only depends on the dimension $d$ and
possibly on other controlled quantities to be specified. We use the notation
$\lesssim$ for $\leq C\times$ up to such a multiplicative constant $C$. We
write $\ll$ (resp. $\gg$) for $\leq C\times$ (resp. $\geq C\times$) up to a
sufficiently large multiplicative constant $C$. When needed, we add subscripts
to indicate dependence on other parameters.

  2. $\bullet$

For a vector field $u$ and a matrix field $S$, we set $(\nabla
u)_{ij}:=\nabla_{j}u_{i}$, $S^{T}_{ij}:=S_{ji}$,
$\operatorname{D}(u):=\frac{1}{2}(\nabla u+(\nabla u)^{T})$, and
${\operatorname{div}}(S)_{i}:=\nabla_{j}S_{ij}$ (we systematically use
Einstein’s summation convention on repeated indices).

  3. $\bullet$

We denote by $\mathrm{d}n$ the (not normalized) Lebesgue measure on the
$(d-1)$-dimensional unit sphere $\mathbb{S}^{d-1}$, and we denote its area by
$\omega_{d}:=|\mathbb{S}_{d-1}|$. Differential operators with a subscript $n$
(such as $\mathrm{div}_{n}$ and $\Delta_{n}$) refer to differential operators
on $\mathbb{S}^{d-1}$, endowed with the natural Riemannian metric.

  4. $\bullet$

For $n\in\mathbb{S}^{d-1}$, we denote by
$\pi_{n}^{\bot}:=\operatorname{Id}-n\otimes n$ the orthogonal projection on
$n^{\bot}$.

  5. $\bullet$

We let $\langle
g\rangle(x):=\fint_{\mathbb{S}^{d-1}}g(x,n)\,\mathrm{d}n:=\frac{1}{\omega_{d}}\int_{\mathbb{S}^{d-1}}g(x,n)\,\mathrm{d}n$
be the angular averaging of a function $g$ on
$\mathbb{T}^{d}\times\mathbb{S}^{d-1}$. We also use the short-hand notation
$P_{1}^{\bot}g:=g-\langle g\rangle$.

  6. $\bullet$

We denote by $H^{k}(\mathbb{T}^{d})$ (resp.
$H^{k}(\mathbb{T}^{d}\times\mathbb{S}^{d-1})$) the standard
$\operatorname{L}^{2}$ Sobolev spaces for functions depending on
$x\in\mathbb{T}^{d}$ (resp. $(x,n)\in\mathbb{T}^{d}\times\mathbb{S}^{d-1}$),
and we use the notation $\|\cdot\|_{H^{k}_{x}}$ (resp.
$\|\cdot\|_{H^{k}_{x,n}}$) for the corresponding norms. For time-dependent
functions, given a Banach space $X$ and $t>0$, the norm of
$\operatorname{L}^{p}(0,t;X)$ is denoted by
$\|\cdot\|_{\operatorname{L}_{t}^{p}X}$.

  7. $\bullet$

The space of probability measures on $\mathbb{T}^{d}$ (resp. on
$\mathbb{T}^{d}\times\mathbb{S}^{d-1}$) is denoted by
$\mathcal{P}(\mathbb{T}^{d})$ (resp. by
$\mathcal{P}(\mathbb{T}^{d}\times\mathbb{S}^{d-1})$).

##  2\. Doi–Saintillan–Shelley theory

For simplicity, we focus in this work on systems in a finite box
$\mathbb{T}^{d}=[0,1)^{d}$ with periodic boundary conditions. In terms of the
fluid velocity field
$u\colon\mathbb{R}^{+}\times\mathbb{T}^{d}\to\mathbb{R}^{d}$ and of the
particles’ position and orientation density function
$f\colon\mathbb{R}^{+}\times\mathbb{T}^{d}\times\mathbb{S}^{d-1}\to\mathbb{R}^{+}$,
we consider the so-called Doi–Saintillan–Shelley kinetic model for a
suspension of very elongated, rigid, active particles in an incompressible
viscous fluid flow (see [DE88, Chapter 8] and [Sai18]),

$\left\\{\begin{array}[]{l}\rho_{\operatorname{fl}}(\partial_{t}+u\cdot\nabla)u-\mu_{\operatorname{fl}}\Delta
u+\nabla
p\,=\,h+{\operatorname{div}}(\sigma_{1}[f])+{\operatorname{div}}(\sigma_{2}[f,\nabla
u]),\\\\[5.69054pt]
\partial_{t}f+{\operatorname{div}}_{x}\big{(}(u+V_{0}n)f\big{)}+{\operatorname{div}}_{n}\big{(}\pi_{n}^{\bot}(\nabla
u)nf\big{)}\\\ \hskip
142.26378pt=D_{\operatorname{tr}}{\operatorname{div}}_{x}\big{(}(\operatorname{Id}+n\otimes
n)\nabla_{x}f\big{)}+D_{\operatorname{ro}}\Delta_{n}f,\\\\[5.69054pt]
{\operatorname{div}}(u)=0,\end{array}\right.$ (2.1)

where elastic and viscous stresses are respectively given by

$\displaystyle\sigma_{1}[f]$ $\displaystyle:=$
$\displaystyle\left(3k_{B}\Theta+\alpha\mu_{\operatorname{fl}}|V_{0}|\ell^{2}\right)\int_{\mathbb{S}^{d-1}}\big{(}n\otimes
n-\tfrac{1}{d}\operatorname{Id}\big{)}f(\cdot,n)\,\mathrm{d}n,$
$\displaystyle\sigma_{2}[f,\nabla u]$ $\displaystyle:=$
$\displaystyle\tfrac{1}{2}\zeta_{\operatorname{ro}}\int_{\mathbb{S}^{d-1}}(n\otimes
n):\operatorname{D}(u)(n\otimes n)\,f(\cdot,n)\,\mathrm{d}n,$ (2.2)

where $h$ is an internal force, and where the different parameters
$\rho_{\operatorname{fl}},\mu_{\operatorname{fl}},D_{\operatorname{tr}},D_{\operatorname{ro}},\zeta_{\operatorname{ro}},V_{0},\ell$
are all assumed to be constant. In this kinetic description,
$\rho_{\operatorname{fl}}$ stands for the fluid density,
$\mu_{\operatorname{fl}}$ for the solvent viscosity, $D_{\operatorname{tr}}$
and $D_{\operatorname{ro}}$ for the translational and rotational diffusion
coefficients, $\zeta_{\operatorname{ro}}$ for the rotational resistance
coefficient, $V_{0}$ for the self-propulsion speed of the particles, and
$\ell$ for the length of the particles. The fluid flow is coupled to the
kinetic equation for the particle density via the additional stresses
$\sigma_{1}$ (elastic stress) and $\sigma_{2}$ (viscous stress), which make
the fluid equations non-Newtonian. We briefly describe the structure and
physical origin of these contributions (see Section 1.3 for references):

  1. —

The viscous stress $\sigma_{2}$ arises from the rigidity of suspended
particles in the fluid flow: it is formally understood by homogenization of
the solid phase, viewed as inclusions with infinite shear viscosity in the
fluid. The above expression (2.2) is an approximation for very elongated
particles: for general axisymmetric particles, the viscous stress involves
additional terms depending on the precise shape of the particles (see e.g.
[HL72] for spheroids), but slender body theory indeed shows that it reduces to
the above form in the limit of very elongated particles (see e.g. [KK13,
Section 3.4]).

  2. —

The elastic stress $\sigma_{1}$ contains a passive and an active contribution.
The passive part, proportional to the the Boltzmann constant $k_{B}$ and to
the absolute temperature $\Theta$, is created by the random torques that are
responsible for the rotational Brownian motion of the particles. These torques
indeed create stresses due to the rigidity and anisotropy of the particles. On
the other hand, the active part arises directly from the swimming mechanism,
which is encoded in the parameter $\alpha\in\mathbb{R}$: a so-called puller
particle corresponds to $\alpha>0$, and a pusher particle corresponds to
$\alpha<0$.

For $V_{0}=0$, the model (2.1) reduces to the classical Doi model for passive
Brownian particles [DE88]. Finally, for very elongated particles in 3D, we
also note that the translational and rotational diffusion and resistance
coefficients are asymptotically given as follows (see [Dho96, Section 5.15]),

$\begin{array}[]{rllrrll}D_{\operatorname{ro}}&=&\tfrac{k_{B}\Theta}{\zeta_{\operatorname{ro}}}+D_{\operatorname{act}},&&\zeta_{\operatorname{ro}}&=&\tfrac{\pi\mu_{\operatorname{fl}}\ell^{3}}{3\log(\ell/a)},\\\\[5.69054pt]
D_{\operatorname{tr}}&=&\tfrac{k_{B}\Theta}{\zeta_{\operatorname{tr}}},&&\zeta_{\operatorname{tr}}&=&\tfrac{2\pi\mu_{\operatorname{fl}}\ell}{\log(\ell/a)},\end{array}$
(2.3)

where $a$ is the width of the particles and where $D_{\operatorname{act}}$ is
some possible active contribution to the rotational diffusion (tumbling).
Compared to the above model (2.1), we henceforth make two minor
simplifications:

  1. —

While the translational diffusion in (2.1) is proportional to
$\operatorname{Id}+n\otimes n$, hence is twice as strong in the direction of
particle orientation as in the orthogonal directions, we choose to neglect
this $O(1)$ difference and assume that the diffusion is isotropic. This choice
is for simplicity and does not change anything in the analysis of the model.

  2. —

It has been observed in the seminal work [BB72] that, for the example of E.
coli bacteria and related microswimmers, the contribution of thermal rotation
and active tumbling is of the same order, and we therefore set
$D_{\operatorname{act}}=0$ for simplicity.

###  2.1. Non-dimensionalization and relevant regimes

We non-dimensionalize the above model (2.1) (after the two above-described
minor simplifications), in terms of the typical speed $u_{0}$ of the fluid,
the typical macroscopic length scale $L$, and the typical number $N$ of
particles in a cube of side length $L$. We further rescale time according to
the time scale of the fluid flow, that is, $T=L/u_{0}$. More precisely, we
define

$\displaystyle\hat{u}(t,x):=\tfrac{1}{u_{0}}u(Tt,Lx),$
$\displaystyle\hat{f}(t,x):=\tfrac{L^{d}}{N}f(Tt,Lx),$
$\displaystyle\hat{h}(t,x):=\tfrac{L^{2}}{\mu u_{0}}h(Tt,Lx).$

This leads to the following dimensionless model, dropping the hats for
simplicity,

$\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+u\cdot\nabla)u-\Delta
u+\nabla
p\,=\,h+{\operatorname{div}}(\sigma_{1}[f])+{\operatorname{div}}(\sigma_{2}[f,\nabla
u]),\\\\[5.69054pt]
\partial_{t}f+{\operatorname{div}}_{x}\big{(}(u+U_{0}n)f\big{)}+{\operatorname{div}}_{n}\big{(}\pi_{n}^{\bot}(\nabla
u)nf\big{)}=\frac{1}{\operatorname{Pe}}\Delta_{x}f+\frac{1}{\operatorname{Wi}}\Delta_{n}f,\\\\[5.69054pt]
{\operatorname{div}}(u)=0,\end{array}\right.$

where the dimensionless counterparts of the additional stresses
$\sigma_{1},\sigma_{2}$ now take the form

$\displaystyle\sigma_{1}[f]$ $\displaystyle=$
$\displaystyle(6\tfrac{\lambda}{\operatorname{Wi}}+\gamma)\int_{\mathbb{S}^{d-1}}\big{(}n\otimes
n-\tfrac{1}{d}\operatorname{Id}\big{)}f(\cdot,n)\,\mathrm{d}n,$
$\displaystyle\sigma_{2}[f,\nabla u]$ $\displaystyle=$
$\displaystyle\lambda\int_{\mathbb{S}^{d-1}}(n\otimes n)(\nabla u)(n\otimes
n)\,f(\cdot,n)\,\mathrm{d}n.$

Here, $\operatorname{Re},\operatorname{Pe},\operatorname{Wi}>0$ stand for the
so-called Reynolds, Péclet, and Weissenberg numbers, $\lambda>0$ depends only
on the shape and number density of the particles, and $\gamma\in\mathbb{R}$
accounts for the active contribution to the stress. More precisely, these
parameters are given by

$\begin{array}[]{rllrrll}\operatorname{Re}&=&\tfrac{\rho_{\operatorname{fl}}u_{0}L}{\mu_{\operatorname{fl}}},&&\lambda&=&\tfrac{\zeta_{\operatorname{ro}}N}{2\mu_{\operatorname{fl}}L^{d}},\\\\[5.69054pt]
\operatorname{Pe}&=&\tfrac{u_{0}L\zeta_{\operatorname{tr}}}{k_{B}\Theta},&&\operatorname{Wi}&=&\tfrac{u_{0}\zeta_{\operatorname{ro}}}{k_{B}\Theta
L},\\\\[5.69054pt] \gamma&=&\alpha
N\tfrac{|U_{0}|\ell^{2}}{L^{d-1}},&&U_{0}&=&\tfrac{V_{0}}{u_{0}},\end{array}$

and we briefly comment on their range and interpretation:

  1. —

The Weissenberg number $\operatorname{Wi}$ is the ratio between convection and
relaxation time scales. For the kinetic viscoelastic models under
consideration, it coincides with the rotational Péclet number. For elongated
particles in 3D, as $\zeta_{\operatorname{ro}}$ is proportional to the cube of
the particle length, cf. (2.3), the regime when $\operatorname{Wi}$ is of
order $O(1)$ is very narrow, and lies under standard flow conditions at
particle lengths of around $10$ microns. In this work, we are interested in
the derivation of hydrodynamic approximations in case of very small
$\operatorname{Wi}\ll 1$: this means that we have applications in mind where
the particles have a length of a few microns and below, which is in particular
the case for many types of bacteria. For notational convenience, we rename the
Weissenberg number as

$\varepsilon\,:=\,\operatorname{Wi}\ll 1.$

  2. —

The (translational) Péclet number $\operatorname{Pe}$ is typically much larger
than its rotational counterpart $\operatorname{Wi}$ since111The same holds
with an additional logarithmic correction in 2D.

$\tfrac{\operatorname{Pe}}{\operatorname{Wi}}\,=\,\tfrac{\zeta_{\operatorname{tr}}L^{2}}{\zeta_{\operatorname{ro}}}\,=\,\tfrac{6L^{2}}{\ell^{2}}\,\gg\,1.$

In fact, from the application-oriented perspective, it makes sense to consider
$\operatorname{Pe}\gg 1$. However, due to well-posedness and stability issues,
our analysis crucially relies on keeping $\operatorname{Pe}$ not too large,
and we shall generally assume $\operatorname{Pe}\simeq 1$. We can actually
allow for $\operatorname{Pe}$ to slightly diverge, but more slowly than the
inverse of the volume fraction of the particles, cf. (1.6).

  3. —

The shape parameter $\lambda$ is typically quite small as it is proportional
to the volume fraction $NL^{-d}$ of the particles.

  4. —

The prefactor in the viscous stress $\sigma_{1}$ reads

$6\tfrac{\lambda}{\operatorname{Wi}}+\gamma\,=\,\tfrac{\lambda}{\operatorname{Wi}}\Big{(}6+2\alpha\tfrac{|V_{0}|\ell^{2}\mu_{\operatorname{fl}}}{k_{B}\Theta}\Big{)},$

where the term $\frac{|V_{0}|\ell^{2}\mu_{\operatorname{fl}}}{k_{B}\Theta}$ is
of order $1-10$ for typical microswimmers like E. coli bacteria [BB72]. We
shall set for abbreviation

$\displaystyle\theta\,:=\,6+2\alpha\tfrac{|V_{0}|\ell^{2}\mu_{\operatorname{fl}}}{k_{B}\Theta}.$
(2.4)

Note in particular that for passive particles $V_{0}=0$ we find $\theta=6$.

  5. —

The self-propulsion speed $V_{0}$ of active particles is typically around $10$
micron per second. This is so slow with respect to typical shear flows
considered in experimental settings that the ratio $U_{0}=\frac{V_{0}}{u_{0}}$
is typically tiny. However, extremely low shear rates leading to $u_{0}\sim
V_{0}$ can possibly also be achieved, as for instance in the experiments
reported by [Lóp+15]. For that reason, we choose not to neglect $U_{0}$ in the
equations and to keep track of its effects.

In conclusion, we are led to consider the following system,

$\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+u_{\varepsilon}\cdot\nabla)u_{\varepsilon}-\Delta
u_{\varepsilon}+\nabla
p_{\varepsilon}=h+\frac{1}{\varepsilon}{\operatorname{div}}(\sigma_{1}[f_{\varepsilon}])+{\operatorname{div}}(\sigma_{2}[f_{\varepsilon},\nabla
u_{\varepsilon}]),\\\\[5.69054pt]
\partial_{t}f_{\varepsilon}+{\operatorname{div}}_{x}\big{(}(u_{\varepsilon}+U_{0}n)f_{\varepsilon}\big{)}+{\operatorname{div}}_{n}\big{(}\pi_{n}^{\bot}(\nabla
u_{\varepsilon})nf_{\varepsilon}\big{)}=\tfrac{1}{\operatorname{Pe}}\Delta_{x}f_{\varepsilon}+\tfrac{1}{\varepsilon}\Delta_{n}f_{\varepsilon},\\\\[5.69054pt]
{\operatorname{div}}(u_{\varepsilon})=0,\\\\[5.69054pt]
\int_{\mathbb{T}^{d}}u_{\varepsilon}=0\leavevmode\nobreak\
\leavevmode\nobreak\ \text{if $\operatorname{Re}=0$},\\\\[5.69054pt]
\sigma_{1}[f]=\lambda\theta\int_{\mathbb{S}^{d-1}}\big{(}n\otimes
n-\tfrac{1}{d}\operatorname{Id}\big{)}f(\cdot,n)\,\mathrm{d}n,\\\\[5.69054pt]
\sigma_{2}[f,\nabla u]=\lambda\int_{\mathbb{S}^{d-1}}(n\otimes n)(\nabla
u)(n\otimes n)\,f(\cdot,n)\,\mathrm{d}n,\end{array}\right.$ (2.5)

which is complemented by initial conditions

$\left\\{\begin{array}[]{l}f_{\varepsilon}|_{t=0}=f^{\circ}_{\varepsilon},\\\\[2.84526pt]
u_{\varepsilon}|_{t=0}=u^{\circ}\leavevmode\nobreak\ \leavevmode\nobreak\
\text{if $\operatorname{Re}\neq 0$},\end{array}\right.$ (2.6)

and we consider the asymptotic limit $\varepsilon\downarrow 0$ in the regime
with $\lambda\leq 1$ and with $\lambda\theta\ll 1$ small enough. Regarding the
Reynolds number, we focus on $\operatorname{Re}\in\\{0,1\\}$ for simplicity.
Due to classical regularity issues for the Navier–Stokes equations, we
actually limit ourselves to the Stokes case $\operatorname{Re}=0$ for the 3D
model, but we also study the Navier–Stokes case $\operatorname{Re}=1$ in 2D.
Note that the 3D Navier–Stokes case is much more complicated and is not
discussed in this work as the well-posedness of the system is still open in
that case.

###  2.2. Well-posedness of the Doi–Saintillan–Shelley model

The system (2.5) is an instance of the more general class of
Fokker–Planck–Navier–Stokes systems, but we emphasize two main peculiarities:

  1. —

We include in (2.5) the contribution of the viscous stress $\sigma_{2}$, which
arises from the rigidity of underlying suspended particles on the microscale
and effectively modifies the solvent viscosity.

  2. —

We also include the effect of particle swimming via $U_{0}$. This creates
local changes in the spatial density
$\rho_{\varepsilon}=\fint_{\mathbb{S}^{d-1}}f_{\varepsilon}(\cdot,n)\,\mathrm{d}n$,
which no longer remains constant in general.

In contrast, most previous works have focused on the corresponding model
without viscous stress $\sigma_{2}$ and without particle swimming $U_{0}=0$.
In that simplified setting, for the 3D Stokes case, the existence of global
entropy solutions was proven in [OT08], and the global well-posedness of
smooth solutions was proven in [CS10] (without translational diffusion,
$\operatorname{Pe}=\infty$). In the Navier–Stokes case, corresponding well-
posedness results were obtained for instance in [CM08]. The model including
the viscous stress $\sigma_{2}[f_{\varepsilon},\nabla u_{\varepsilon}]$ but
without particle swimming $U_{0}=0$ was first studied in [LM07], where the
existence of global weak solutions was proven for the Navier–Stokes case in 2D
and 3D (without translational diffusion, $\operatorname{Pe}=\infty$). We also
refer to [La19] for the global well-posedness of smooth solutions in the 2D
Navier–Stokes case. Particle swimming $U_{0}\neq 0$ was first considered in
[CL13], where the authors studied the corresponding model without viscous
stress $\sigma_{2}$ and proved the global existence of weak entropy solutions
both for the Stokes and Navier–Stokes cases in 2D and 3D, as well as the
existence of energy solutions for the Stokes case and their uniqueness in 2D.

Building on similar ideas, we can actually establish the existence of global
energy solutions for the full model (2.5) in the 2D and 3D Stokes cases, as
well as in the 2D Navier–Stokes case, and we further obtain a weak-strong
uniqueness principle. To our knowledge, this is surprisingly the first result
where both the viscous stress and the swimming forces are included at the same
time. The proof is postponed to Appendix A.

######  Proposition 2.1.

Consider either the Stokes case $\operatorname{Re}=0$ with $d\leq 3$, or the
Navier–Stokes case $\operatorname{Re}=1$ with $d=2$. Given $\varepsilon>0$,
given
$h\in\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};\operatorname{L}^{2}(\mathbb{T}^{d}))$,
given an initial condition
$f^{\circ}_{\varepsilon}\in\operatorname{L}^{2}\cap\mathcal{P}(\mathbb{T}^{d}\times\mathbb{S}^{d-1})$,
and given also $u^{\circ}\in\operatorname{L}^{2}(\mathbb{T}^{d})^{d}$ with
${\operatorname{div}}(u^{\circ})=0$ in the Navier–Stokes case, the Cauchy
problem (2.5)–(2.6) admits a global weak solution
$(u_{\varepsilon},f_{\varepsilon})$ with:

  1. (i)

in the Stokes case $\operatorname{Re}=0$, $d\leq 3$,

$\displaystyle u_{\varepsilon}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{1}(\mathbb{T}^{d})^{d}),$
$\displaystyle f_{\varepsilon}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};\operatorname{L}^{2}\cap\mathcal{P}(\mathbb{T}^{d}\times\mathbb{S}^{d-1}))\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{1}(\mathbb{T}^{d}\times\mathbb{S}^{d-1}));$

  2. (ii)

in the Navier–Stokes case $\operatorname{Re}=1$, $d=2$,

$\displaystyle u_{\varepsilon}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};\operatorname{L}^{2}(\mathbb{T}^{2})^{2})\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{1}(\mathbb{T}^{2})^{2}),$
$\displaystyle f_{\varepsilon}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};\operatorname{L}^{2}\cap\mathcal{P}(\mathbb{T}^{2}\times\mathbb{S}^{1}))\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{1}(\mathbb{T}^{2}\times\mathbb{S}^{1})).$

In both cases, a weak-strong uniqueness principle further holds: if
$(u_{\varepsilon},f_{\varepsilon})$ and
$(u^{\prime}_{\varepsilon},f^{\prime}_{\varepsilon})$ are two such global weak
solutions with identical initial conditions, and if
$(u^{\prime}_{\varepsilon},f^{\prime}_{\varepsilon})$ has the following
additional regularity,

$\displaystyle u^{\prime}_{\varepsilon}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{2}_{\mathrm{loc}}(\mathbb{R}^{+};W^{1,\infty}(\mathbb{T}^{d})^{d}),$
$\displaystyle f^{\prime}_{\varepsilon}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\mathrm{loc}}(\mathbb{R}^{+};\operatorname{L}^{\infty}(\mathbb{T}^{d}\times\mathbb{S}^{d-1})),$

then we have
$(u_{\varepsilon},f_{\varepsilon})=(u^{\prime}_{\varepsilon},f^{\prime}_{\varepsilon})$.
∎

##  3\. Ordered fluid models

Various non-Newtonian fluid models have been considered in the literature,
taking into account nonlinear and memory effects in different ways on the
macroscopic scale. We focus here on so-called ordered fluid models, which will
indeed be shown to naturally appear as hydrodynamic approximations of the
Doi–Saintillan–Shelley theory. We start by defining such models and reviewing
their non-Newtonian properties; for more details on modelling aspects and
applications, we refer to [DR95, Böh87, BAH87, Jos13, PTMD13]. Next, we
develop a perturbative well-posedness theory and we introduce (non-standard)
inhomogeneous, diffusive versions of these models.

###  3.1. Standard ordered fluid models

Ordered fluid equations are of the form

$\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+u\cdot\nabla)u-{\operatorname{div}}(\sigma)+\nabla
p=h,\\\\[2.84526pt] {\operatorname{div}}(u)=0,\end{array}\right.$ (3.1)

where $\sigma$ is a function of the Rivlin–Ericksen tensors
$\\{A_{k}(u)\\}_{k}$ defined iteratively through

$\displaystyle A_{1}(u)\,:=\,2\\!\operatorname{D}(u)\,=\,\nabla u+(\nabla
u)^{T},$ $\displaystyle
A_{n+1}(u)\,:=\,(\partial_{t}+u\cdot\nabla)A_{n}(u)+(\nabla
u)^{T}A_{n}(u)+A_{n}(u)\nabla u,\qquad n\geq 1.$ (3.2)

In other words, $A_{n+1}(u)$ is the so-called lower-convected derivative of
$A_{n}(u)$, which ensures frame indifference of the equations.222In the
literature, _upper_ -convected instead of lower-convected derivatives are
sometimes used to define ordered fluids (see e.g. [BCAH87, Section 6]). This
is merely a choice of convention, as both lead to equivalent fluid equations
(although the value and interpretation of parameters of course depends on the
chosen convention). Indeed, the upper-convected derivative of $A_{n}(u)$ can
be rewritten in terms of the lower-convected derivative as
$(\partial_{t}+u\cdot\nabla)A_{n}(u)-A_{n}(u)(\nabla u)^{T}-(\nabla
u)A_{n}(u)\,=\,A_{n+1}(u)-\big{(}A_{1}(u)A_{n}(u)+A_{n}(u)A_{1}(u)\big{)}.$
While first-order fluids are of the form (3.1) with constitutive law
$\sigma=\eta_{0}A_{1}(u)$, thus coinciding with standard Newtonian fluids,
second-order fluids are of the form (3.1) with

$\displaystyle\sigma\,:=\,\eta_{0}A_{1}(u)+\alpha_{1}A_{2}(u)+\alpha_{2}A_{1}(u)^{2},$
(3.3)

see e.g. [Böh87, Eq. (8.48)], and third-order fluids amount to

$\sigma\,:=\,\eta_{0}A_{1}(u)+\alpha_{1}A_{2}(u)+\alpha_{2}A_{1}(u)^{2}\\\
+\beta_{1}A_{3}(u)+\beta_{2}\big{(}A_{1}(u)A_{2}(u)+A_{2}(u)A_{1}(u)\big{)}+\beta_{3}A_{1}(u){\operatorname{tr}}\big{(}A_{1}(u)^{2}\big{)},$
(3.4)

for some coefficients
$\eta_{0},\alpha_{1},\alpha_{2},\beta_{1},\beta_{2},\beta_{3}\in\mathbb{R}$.333In
arbitrary dimension $d$, a further term $\beta_{4}A_{1}(u)^{3}$ should be
included in general in the stress $\sigma$ of third-order fluids, but it is
redundant in dimensions $d\leq 3$ as we then have
$B^{3}=\frac{1}{2}B\,{\operatorname{tr}}(B^{2})+\frac{1}{3}{\operatorname{tr}}(B^{3})$
by the Cayley–Hamilton theorem for any symmetric trace-free matrix $B$. Since
third- and higher-order fluid models are fairly complicated and involve a lot
of terms, the most widely used ordered fluid model for viscoelastic fluids is
the second-order model (3.3).

We briefly recall the basic rheological properties of these ordered fluid
models depending on the different parameter values; we focus here on the 3D
setting. The coefficient $\eta_{0}>0$ in (3.3) and (3.4) is the zero-shear
viscosity, as in the first-order (Newtonian) model, while other parameters
account for various non-Newtonian behaviors.

  1. —

Shear-dependent viscosity: In a simple shear flow $u_{0}(t,x)=\kappa
x_{2}e_{1}$ with shear rate $\kappa>0$, the shear-dependent viscosity is
defined as

$\eta(\kappa)\,:=\,\frac{\sigma_{12}}{\kappa}.$ (3.5)

Noting that

$\displaystyle\qquad A_{1}(u_{0})=\kappa\begin{pmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix},$ $\displaystyle
A_{1}(u_{0})^{2}=\kappa^{2}\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&0\end{pmatrix},$ $\displaystyle
A_{2}(u_{0})=\kappa^{2}\begin{pmatrix}0&0&0\\\ 0&2&0\\\ 0&0&0\end{pmatrix},$

we compute for the second-order fluid that the shear-dependent viscosity
simply coincides with the zero-shear viscosity, $\eta(\kappa)=\eta_{0}$. For
the third-order model, in contrast, we find a nontrivial shear-dependent
relation,

$\eta(\kappa)=\eta_{0}+2\kappa^{2}(\beta_{2}+\beta_{3}).$

Most real-life viscoelastic fluids, and in particular passive dilute
suspensions, happen to be shear-thinning, meaning that the map
$\kappa\mapsto\eta(\kappa)$ is decreasing: this holds for the third-order
fluid model provided that the coefficients satisfy $\beta_{2}+\beta_{3}<0$.

  2. —

Normal-stress differences: In a simple shear flow $u_{0}(t,x)=\kappa
x_{2}e_{1}$ with shear rate $\kappa\in\mathbb{R}$, non-Newtonian fluids
typically display non-zero normal stresses. This is responsible for a number
of phenomena, of which the rod-climbing effect is the best known; see e.g.
[Jos13, Chapter 17] and [BAH87, Chapter 6]. Normal stress coefficients are
defined as444Beware of different sign conventions for the normal stress
coefficients. We follow here the definition of [Böh87, Chapter 2.2].

$\displaystyle\nu_{10}\,:=\,\frac{\sigma_{11}-\sigma_{22}}{\kappa^{2}},$
$\displaystyle\nu_{20}\,:=\,\frac{\sigma_{22}-\sigma_{33}}{\kappa^{2}},$ (3.6)

and thus, for second- and third-order models,

$\displaystyle\nu_{10}\,=\,-2\alpha_{1},$
$\displaystyle\nu_{20}\,=\,2\alpha_{1}+\alpha_{2}.$ (3.7)

In other words, $\alpha_{1},\alpha_{2}$ are related to normal stress
coefficients via

$\displaystyle\alpha_{1}\,=\,-\tfrac{1}{2}\nu_{10},$
$\displaystyle\alpha_{2}\,=\,\nu_{10}+\nu_{20}.$ (3.8)

For most real-life viscoelastic fluids, and in particular for polymer
solutions, it is found experimentally that $\nu_{10}>0$, $\nu_{20}<0$, and
that $|\nu_{20}|$ is considerably smaller than $|\nu_{10}|$ (by a factor of
around $10$, see e.g. [Böh87, Section 2.2] and [PTMD13, Section2.2]), which
means in particular, in terms of second-order fluid coefficients,

$\displaystyle\alpha_{1}<0,$ $\displaystyle\alpha_{2}>0.$ (3.9)

  3. —

Elongational viscosity: The apparent viscosity of a non-Newtonian fluid may be
completely different in an elongational flow. Given a uniaxial elongational
flow $u_{0}(t,x)=\kappa(x_{1}e_{1}-\frac{1}{2}(x_{2}e_{2}+x_{3}e_{3}))$ in the
direction $e_{1}$, the elongational viscosity is defined as

$\displaystyle\eta_{E}\,:=\,\frac{\sigma_{11}-\frac{1}{2}(\sigma_{22}+\sigma_{33})}{\kappa}.$
(3.10)

Noting that

$\displaystyle A_{1}(u_{0})=\kappa\begin{pmatrix}2&0&0\\\ 0&-1&0\\\
0&0&-1\end{pmatrix},$ $\displaystyle
A_{2}(u_{0})=A_{1}(u_{0})^{2}=\kappa^{2}\begin{pmatrix}4&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix},$ (3.11)

we compute for the second-order model,

$\displaystyle\eta_{E}\,=\,3\eta_{0}+3\kappa(\alpha_{1}+\alpha_{2}).$

For real-life viscoelastic fluids, it is typically observed that the
elongational viscosity increases with the strain rate (so-called strain-
thickening behavior), which holds provided that coefficients satisfy
$\alpha_{1}+\alpha_{2}>0$.

  4. —

Retardation phase shift in oscillatory flow: In a simple shear flow
$u_{0}(t,x)=\kappa(t)x_{2}e_{1}$ with oscillatory shear rate $\kappa(t)=\sin
t$, we compute

$A_{2}(u_{0})=\dot{\kappa}\begin{pmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix}+\kappa^{2}\begin{pmatrix}0&0&0\\\ 0&2&0\\\
0&0&0\end{pmatrix},$

which gives rise to a phase shift for $\sigma_{12}$ in the second-order model,
in form of

$\sigma_{12}\,=\,\eta_{0}\kappa+\alpha_{1}\dot{\kappa}\,=\,(\eta_{0}^{2}+\alpha_{1}^{2})^{\frac{1}{2}}\sin\big{(}t+\arctan(\tfrac{\alpha_{1}}{\eta_{0}})\big{)}.$

This constitutes another typical (time-dependent) non-Newtonian feature, in
link with the dependence of the stress on the flow history.

######  Remark 3.1 (Connection to Oldroyd–B model).

The Oldroyd–B model is a special case of a simple fluid model, which is
particularly popular as a formal exact closure of the kinetic Hookean dumbbell
model. It is characterized by the following constitutive equation for the
stress tensor,

$\displaystyle\sigma=2\eta_{s}\operatorname{D}(u)+\eta_{p}\tau,$
$\displaystyle\tau+\operatorname{Wi}\left((\partial_{t}+u\cdot\nabla)\tau-\tau(\nabla
u)^{T}-(\nabla u)\tau\right)\,=\,2\operatorname{D}(u).$ (3.12)

By a formal expansion with respect to $\operatorname{Wi}\ll 1$, this model is
found to agree to order $\operatorname{Wi}^{k}$ with a $k$th-order fluid model
with some specific choice of parameters. In particular, to order
$\operatorname{Wi}^{2}$, we recover the second-order fluid model (3.3) with
$\alpha_{2}=-2\alpha_{1}$, which means in particular that the second normal
stress coefficient vanishes, $\nu_{20}=0$. We emphasize that this is not the
case for the second-order fluid model that we shall derive from the
Doi–Saintillan–Shelley theory. Hence, at small $\operatorname{Wi}$, the
viscoelastic effects of suspensions of rigid Brownian particles are not well
described by an Oldroyd–B model. ∎

###  3.2. Non-standard ordered fluid models at $\operatorname{Pe}<\infty$

In the case of a particle suspension with finite Péclet number
$\operatorname{Pe}<\infty$, as considered in this work, cf. (2.5), the
Rivlin–Ericksen tensors (3.2) naturally need to be modified as follows,

$\displaystyle A_{1}^{\prime}(u)\,:=\,A_{1}(u)\,=\,2\operatorname{D}(u),$
(3.13) $\displaystyle
A_{n+1}^{\prime}(u)\,:=\,(\partial_{t}-{\tfrac{1}{\operatorname{Pe}}}\Delta+u\cdot\nabla)A_{n}^{\prime}(u)+(\nabla
u)^{T}A_{n}^{\prime}(u)+A_{n}^{\prime}(u)(\nabla u),\qquad n\geq 1,$

hence in particular
$A_{2}^{\prime}(u):=A_{2}(u)-\frac{1}{\operatorname{Pe}}\Delta A_{1}(u)$. The
second-order fluid equations (3.1)–(3.3) are then replaced by

$\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+u\cdot\nabla)u-{\operatorname{div}}(\sigma)+\nabla
p=h,\\\\[2.84526pt]
\sigma=\eta_{0}A_{1}(u)+\alpha_{1}A_{2}^{\prime}(u)+\alpha_{2}A_{1}(u)^{2},\\\\[2.84526pt]
{\operatorname{div}}(u)=0.\end{array}\right.$ (3.14)

For higher-order fluid models, in this diffusive setting
$\operatorname{Pe}<\infty$, several additional tensors actually need to be
included in the stress: for the third-order model, instead of (3.4), the
stress rather needs to be chosen in general as

$\sigma\,=\,\eta_{0}A_{1}(u)+\alpha_{1}A_{2}^{\prime}(u)+\alpha_{2}A_{1}(u)^{2}+\beta_{1}A_{3}^{\prime}(u)+\beta_{1}^{\prime}B_{3}^{\prime}(u)\\\
+\beta_{2}\big{(}A_{1}(u)A_{2}^{\prime}(u)+A_{2}^{\prime}(u)A_{1}(u)\big{)}+\beta_{3}A_{1}(u){\operatorname{tr}}\big{(}A_{1}(u)^{2}\big{)},$
(3.15)

in terms of the following additional tensor,

$B_{3}^{\prime}(u)\,:=\,(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta+u\cdot\nabla)A_{1}(u)^{2}+\big{(}(\nabla
u)^{T}A_{1}(u)^{2}+A_{1}(u)^{2}(\nabla u)\big{)}.$

At infinite Péclet number $\operatorname{Pe}=\infty$, this additional tensor
is indeed redundant as it reduces to
$B_{3}^{\prime}(u)=A_{1}(u)A_{2}(u)+A_{2}(u)A_{1}(u)-A_{1}(u)^{3}$, so we
recover (3.4).555Recall that, as in (3.4), a further term
$\beta_{4}A_{1}(u)^{3}$ should always be included in the stress in arbitrary
dimension $d$, but it reduces to
$\frac{1}{2}\beta_{4}A_{1}(u){\operatorname{tr}}(A_{1}(u)^{2})$ in dimensions
$d\leq 3$.

###  3.3. Well-posedness of ordered fluid models

We focus for shortness on the second-order fluid model. There has actually
been a fair amount of confusion on the relevant range of parameters
$\alpha_{1},\alpha_{2}$: the sign condition $\alpha_{1}<0$ in (3.9) is
motivated by experimental normal stress measurements, but it actually appears
inconsistent with thermodynamics; see e.g. [DR95] for a detailed discussion
from the physics perspective. From the mathematical point of view, this
inconsistency leads to ill-posedness issues. The matter was investigated by
Galdi [Gal95], who showed the following for the second-order fluid equations
(3.1)–(3.3) at infinite Péclet number $\operatorname{Pe}=\infty$:

  1. —

the local-in-time well-posedness holds whenever $1/\alpha_{1}>-\lambda_{1}$,
where $\lambda_{1}$ stands for the Poincaré constant in $\mathbb{T}^{d}$,
which is quite consistent with the choice (3.9) (although the case of a
negative $\alpha_{1}$ with a small absolute value is prohibited);

  2. —

the long-time well-posedness, as well as the stability of steady solutions,
can only hold provided that $\alpha_{1}>0$.

For the corresponding system (3.14) with finite Péclet number
$\operatorname{Pe}<\infty$, the situation is even worse: even local-in-time
well-posedness actually fails whenever $\alpha_{1}<0$ because the equation
then behaves like a backwards heat equation. Yet, even though the kinetic Doi
model itself is known to be thermodynamically consistent (see [DE88, Chapter
8]), our analysis will precisely lead us to a second-order fluid with
$\alpha_{1}<0$, and it is thus crucial to determine what meaning should be
given to the model in that case. In fact, we shall derive (3.14) in the
small-$\operatorname{Wi}$ limit, $\varepsilon:=\operatorname{Wi}\ll 1$, with

$(\alpha_{1},\alpha_{2})=(\varepsilon\gamma_{1},\varepsilon\gamma_{2}),\qquad\text{for
some $\gamma_{1}<0,\leavevmode\nobreak\ \gamma_{2}\in\mathbb{R}$.}$ (3.16)

In this perturbative setting $\varepsilon\ll 1$, although the equation is ill-
posed for fixed $\varepsilon$, there are several ways to rearrange the
nonlinearity and define well-posed notions of approximate solutions that only
satisfy the system (3.14) up to a higher-order $O(\varepsilon^{2})$ remainder:

$\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+\bar{u}_{\varepsilon}\cdot\nabla)\bar{u}_{\varepsilon}-{\operatorname{div}}(\bar{\sigma}_{\varepsilon})+\nabla\bar{p}_{\varepsilon}=h+O(\varepsilon^{2}),\\\\[2.84526pt]
\bar{\sigma}_{\varepsilon}\,:=\,\eta_{0}A_{1}(\bar{u}_{\varepsilon})+\varepsilon\gamma_{1}A_{2}^{\prime}(\bar{u}_{\varepsilon})+\varepsilon\gamma_{2}A_{1}(\bar{u}_{\varepsilon})^{2},\\\\[2.84526pt]
{\operatorname{div}}(\bar{u}_{\varepsilon})=0.\end{array}\right.$ (3.17)

These notions of approximate solutions will be viewed as proper perturbative
ways to interpret second-order fluid equations and settle ill-posedness and
instability issues.

The simplest way is to define a notion of approximate _hierarchical
solutions_. It is particularly convenient for us in this work due to its close
relation to the Hilbert expansion method that we use for the hydrodynamic
approximation. The proof is straightforward and is postponed to Appendix B.

######  Proposition 3.2 (Hierarchical solutions).

Consider the system (3.17) in the regime $\varepsilon\ll 1$ with parameters
$\eta_{0},\operatorname{Pe}>0$ and $\gamma_{1},\gamma_{2}\in\mathbb{R}$. If
$v_{0},v_{1}$ are smooth solutions of the following two auxiliary systems,

$\left\\{\begin{array}[]{l}\operatorname{Re}\big{(}\partial_{t}v_{0}+(v_{0}\cdot\nabla)v_{0}\big{)}-\eta_{0}\Delta
v_{0}+\nabla p_{0}=h,\\\\[1.42262pt] {\operatorname{div}}(v_{0})=0,\\\
v_{0}|_{t=0}=u^{\circ}\quad\text{if $\operatorname{Re}\neq
0$},\quad\int_{\mathbb{T}^{d}}v_{0}=0\quad\text{if
$\operatorname{Re}=0$},\end{array}\right.$ (3.18)
$\left\\{\begin{array}[]{l}\operatorname{Re}\big{(}\partial_{t}v_{1}+(v_{0}\cdot\nabla)v_{1}+(v_{1}\cdot\nabla)v_{0}\big{)}-\eta_{0}\Delta
v_{1}+\nabla p_{1}\\\\[2.84526pt] \hskip
113.81102pt={\operatorname{div}}\big{(}\gamma_{1}A_{2}^{\prime}(v_{0})+\gamma_{2}A_{1}(v_{0})^{2}\big{)},\\\\[1.42262pt]
{\operatorname{div}}(v_{1})=0,\\\\[1.42262pt] v_{1}|_{t=0}=0\quad\text{if
$\operatorname{Re}\neq 0$},\quad\int_{\mathbb{T}^{d}}v_{1}=0\quad\text{if
$\operatorname{Re}=0$},\end{array}\right.$ (3.19)

then the superposition $\bar{u}_{\varepsilon}:=v_{0}+\varepsilon v_{1}$ indeed
satisfies the system (3.17) with some controlled error term
$O(\varepsilon^{2})$ (and with initial condition
$\bar{u}_{\varepsilon}|_{t=0}=u^{\circ}$ if $\operatorname{Re}\neq 0$). For
the well-posedness of (3.18) and (3.19), we separately consider the Stokes and
Navier–Stokes cases:

  1. (i)

_Stokes case $\operatorname{Re}=0$, $d\leq 3$:_
Given $s>\frac{d}{2}-1$ and
$h\in\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{d})^{d})\cap
W^{1,\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s-2}(\mathbb{T}^{d})^{d})$,
there is a unique global solution
$v_{0}\in\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+2}(\mathbb{T}^{d})^{d})$
of (3.18), and a unique global solution
$v_{1}\in\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})^{d})$
of (3.19), leading to $\bar{u}_{\varepsilon}=v_{0}+\varepsilon
v_{1}\in\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})^{d})$.

  2. (ii)

_Navier–Stokes case $\operatorname{Re}=1$, $d=2$:_
Given $s>0$,
$h\in\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{2})^{2})$,
and $u^{\circ}\in H^{s+1}(\mathbb{T}^{2})^{2}$ with
${\operatorname{div}}(u^{\circ})=0$, there is a unique global solution
$v_{0}\in\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{2})^{2})\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+2}(\mathbb{T}^{2})^{2})$
of (3.18), and a unique global solution
$v_{1}\in\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{2})^{2})\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{2})^{2})$
of (3.19), leading to $\bar{u}_{\varepsilon}=v_{0}+\varepsilon
v_{1}\in\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{2})^{2})\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{2})^{2})$.
∎

There are also non-hierarchical ways to perturbatively make sense of the ill-
posed second-order fluid model (3.17), which may be more desirable in
particular for stability issues. Comparing to corresponding ill-posedness
issues in the Boussinesq theory for water waves, we recall that there is a
standard way to rearrange the ill-posed Boussinesq equation perturbatively and
make it well-posed, see [CMV96]: in a nutshell, the idea is to replace
indefinite operators like $1+\varepsilon\Delta$ by corresponding positive
operators like $(1-\varepsilon\Delta)^{-1}$ up to $O(\varepsilon^{2})$ errors.
We show that a similar so-called Boussinesq trick can be used in the present
situation as well: for any value of $\gamma_{1},\gamma_{2}$, both at finite
and infinite Péclet number, it leads us to a perturbative rearrangement of the
second-order fluid equation that is well-posed and is indeed equivalent to
(3.14) up to $O(\varepsilon^{2})$ terms. The so-defined solution is easily
checked to differ from the corresponding hierarchical solution only by
$O(\varepsilon^{2})$. We focus here on the relevant ill-posed case
$\gamma_{1}\leq 0$. Note that the procedure is easily extended to higher-order
fluid equations (see also [ABV16, DGR23] in a different context). The proof is
postponed to Appendix B.

######  Proposition 3.3 (Boussinesq-like solutions).

Consider the system (3.17) in the regime $\varepsilon\ll 1$ with parameters
$\eta_{0},\operatorname{Pe}>0$, $\gamma_{1}\leq 0$, and
$\gamma_{2}\in\mathbb{R}$.

  1. (i)

_Stokes case $\operatorname{Re}=0$, $d\leq 3$:_
Given $s>\frac{d}{2}$, $T_{0}>0$, and
$h\in\operatorname{L}^{\infty}(0,T_{0};H^{s+1}(\mathbb{T}^{d})^{d})\cap
W^{1,\infty}(0,T_{0};H^{s-1}(\mathbb{T}^{d})^{d})$, provided that
$\varepsilon\ll 1$ is small enough (depending on $s,T_{0},h$ and on all
parameters), the following nonlinear problem admits a unique solution
$\bar{u}_{\varepsilon}\in\operatorname{L}^{\infty}(0,T_{0};H^{s+1}(\mathbb{T}^{d})^{d})$,

$\qquad\left\\{\begin{array}[]{l}-\eta_{0}\Delta\bar{u}_{\varepsilon}+\nabla\bar{p}_{\varepsilon}=\big{(}1-\varepsilon\tfrac{\gamma_{1}}{\eta_{0}}(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta)\big{)}h+\varepsilon{\operatorname{div}}(F_{0}(\bar{u}_{\varepsilon})),\\\\[2.84526pt]
{\operatorname{div}}(\bar{u}_{\varepsilon})=0,\quad\int_{\mathbb{T}^{d}}\bar{u}_{\varepsilon}=0,\end{array}\right.$
(3.20)

where we have set for abbreviation

$\qquad
F_{0}(u)\,:=\,\gamma_{1}(u\cdot\nabla)2\\!\operatorname{D}(u)+\gamma_{1}\big{(}(\nabla
u)^{T}2\\!\operatorname{D}(u)+2\\!\operatorname{D}(u)(\nabla
u)\big{)}+\gamma_{2}(2\\!\operatorname{D}(u))^{2}.$ (3.21)

Moreover, the so-defined solution $\bar{u}_{\varepsilon}$ satisfies (3.17)
with $\operatorname{Re}=0$ for some controlled error term
$O(\varepsilon^{2})$.

  2. (ii)

_Navier–Stokes case $\operatorname{Re}=1$, $d=2$:_
Further assume $\eta_{0}\geq\frac{1}{\operatorname{Pe}}$.666For
$\eta_{0}<\frac{1}{\operatorname{Pe}}$, the equations would need to be
rearranged differently; we skip it here for shortness. Given $s>1$, $T_{0}>0$,
$u^{\circ}\in H^{s}(\mathbb{T}^{2})^{2}$ with
${\operatorname{div}}(u^{\circ})=0$, and
$h\in\operatorname{L}^{2}(0,T_{0};H^{s+1}(\mathbb{T}^{2})^{2})$, provided that
$\varepsilon\ll 1$ is small enough (depending on $s,T_{0},u^{\circ},h$ and on
all parameters), the following nonlinear problem admits a unique solution
$\bar{u}_{\varepsilon}\in\operatorname{L}^{\infty}(0,T_{0};H^{s}(\mathbb{T}^{2})^{2})\cap\operatorname{L}^{2}(0,T_{0};H^{s+1}(\mathbb{T}^{2})^{2})$,

$\qquad\left\\{\begin{array}[]{l}(\partial_{t}+\bar{u}_{\varepsilon}\cdot\nabla)\bar{u}_{\varepsilon}-\eta_{0}\Delta\bar{u}_{\varepsilon}-\varepsilon\gamma_{1}(\eta_{0}-\frac{1}{\operatorname{Pe}})\Delta^{2}\bar{u}_{\varepsilon}+\nabla\bar{p}_{\varepsilon}\\\
\hskip
142.26378pt\,=\,(1+\varepsilon\gamma_{1}\Delta)h+\varepsilon{\operatorname{div}}(F_{1}(\bar{u}_{\varepsilon})),\\\
{\operatorname{div}}(\bar{u}_{\varepsilon})=0,\\\
\bar{u}_{\varepsilon}|_{t=0}=u^{\circ},\end{array}\right.$ (3.22)

where we have set for abbreviation

$\qquad F_{1}(u)\,:=\,2\gamma_{1}(\nabla u)^{T}(\nabla
u)+\gamma_{2}(2\\!\operatorname{D}(u))^{2}.$

Moreover, the so-defined solution $\bar{u}_{\varepsilon}$ satisfies (3.17)
with $\operatorname{Re}=1$ with some controlled error term
$O(\varepsilon^{2})$ and with initial condition
$\bar{u}_{\varepsilon}|_{t=0}=u^{\circ}$.∎

###  3.4. Models for inhomogeneous suspensions

The above formulation of ordered fluid models describes the behavior of
spatially homogeneous suspensions, that is, the behavior of fluids with a
constant density of suspended particles. This can be naturally generalized to
an inhomogeneous setting to describe non-uniform particle suspensions.
Surprisingly, we were not able to find any account of this generalization in
the literature. We emphasize however that such an inhomogeneous setting should
also arise naturally from the small-$\operatorname{Wi}$ expansion of the
inhomogeneous Oldroyd–B model, which is the formal exact closure of the
kinetic Hookean dumbbell model. Note that we are still considering an
homogeneous solvent fluid. At infinite Péclet number
$\operatorname{Pe}=\infty$, the adaptation is straightforward: the fluid
equations (3.1) are simply coupled to a conservation equation for the particle
density $\rho:=\fint_{\mathbb{S}^{d-1}}f(\cdot,n)\,\mathrm{d}n$,777The
physical particle density is rather given by
$x\mapsto\int_{\mathbb{S}^{d-1}}f(x,n)\,\mathrm{d}n$, but for notational
convenience we choose to normalize it by the area of $\mathbb{S}^{d-1}$. In
particular, we have $\rho\in\frac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d})$.

$(\partial_{t}+u\cdot\nabla)\rho=0,$

while the stress $\sigma$ is now a function both of $\rho$ and of the
Rivlin–Ericksen tensors $\\{A_{k}(u)\\}_{k}$ where non-Newtonian corrections
to the pure solvent viscosity $\eta_{0}$ are taken proportional to the
suspended particle density $\rho$. More precisely, the inhomogeneous second-
order fluid model takes the form

$\displaystyle\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+u\cdot\nabla)u-{\operatorname{div}}(\sigma)+\nabla
p=h,\\\\[2.84526pt] (\partial_{t}+u\cdot\nabla)\rho=0,\\\\[2.84526pt]
\sigma\,=\,(\eta_{0}+\eta_{1}\rho)A_{1}(u)+\alpha_{1}\rho
A_{2}(u)+\alpha_{2}\rho A_{1}(u)^{2},\\\\[2.84526pt]
{\operatorname{div}}(u)=0,\end{array}\right.$ (3.27)

for some coefficients $\eta_{1},\alpha_{1},\alpha_{2}\in\mathbb{R}$.
Inhomogeneous versions of higher-order fluid models are formulated similarly.

At finite Péclet number $\operatorname{Pe}<\infty$, on the other hand, the
particle density $\rho$ is no longer simply transported by the fluid, but
rather solves a transport-diffusion equation,

$(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta+u\cdot\nabla)\rho=0.$

In this diffusive setting, the structure of ordered fluid models becomes
slightly more complicated: due to diffusion, the transport-diffusion operator
that appears in the Rivlin–Ericksen tensors (3.13) at finite Péclet number
does not commute with multiplication with the particle density $\rho$. The
second-order fluid model then rather takes the form

$\displaystyle\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+u\cdot\nabla)u-{\operatorname{div}}(\sigma)+\nabla
p=h,\\\\[2.84526pt]
(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta+u\cdot\nabla)\rho=0,\\\\[2.84526pt]
\sigma\,=\,(\eta_{0}+\eta_{1}\rho)A_{1}(u)+\alpha_{1}A_{2}^{\prime}(\rho,u)+\alpha_{2}\rho
A_{1}(u)^{2},\\\\[2.84526pt] {\operatorname{div}}(u)=0,\end{array}\right.$
(3.32)

for some coefficients $\eta_{1},\alpha_{1},\alpha_{2}\in\mathbb{R}$, in terms
of the modified inhomogeneous second-order Rivlin–Ericksen tensor

$A_{2}^{\prime}(\rho,u)\,:=\,(\partial_{t}-{\tfrac{1}{\operatorname{Pe}}}\Delta+u\cdot\nabla)(\rho
A_{1}(u))+\rho\big{(}(\nabla u)^{T}A_{1}(u)+A_{1}(u)(\nabla u)\big{)}.$ (3.33)

Indeed, due to the diffusion, the latter quantity does not reduce to the
Rivlin–Ericksen tensor defined in (3.13): we have
$A_{2}^{\prime}(\rho,u)\neq\rho A_{2}^{\prime}(u)$ in general along solutions
— in contrast with the case of infinite Péclet number.888At infinite Péclet
number $\operatorname{Pe}=\infty$, as the particle density $\rho$ satisfies
$(\partial_{t}+u\cdot\nabla)\rho=0$, we indeed obtain
$(\partial_{t}+u\cdot\nabla)(\rho
A_{1}(u))=\rho(\partial_{t}+u\cdot\nabla)A_{1}(u)$, so that the system (3.32)
reduces to (3.27). Similarly, the inhomogeneous third-order fluid model
amounts to (3.32) with stress

$\sigma\,=\,(\eta_{0}+\eta_{1}\rho)A_{1}(u)+\alpha_{1}A_{2}^{\prime}(\rho,u)+\alpha_{2}A_{1}(u)^{2}+\beta_{1}A_{3}^{\prime}(\rho,u)+\beta_{1}^{\prime}B_{3}^{\prime}(\rho,u)\\\
+\beta_{2}\big{(}A_{1}(u)A_{2}^{\prime}(\rho,u)+A_{2}^{\prime}(\rho,u)A_{1}(u)\big{)}+\beta_{3}\rho
A_{1}(u){\operatorname{tr}}(A_{1}(u)^{2}),$

in terms of the modified inhomogeneous third-order Rivlin–Ericksen tensor

$A_{3}^{\prime}(\rho,u)\,:=\,(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta+u\cdot\nabla)A_{2}^{\prime}(\rho,u)+(\nabla
u)^{T}A_{2}^{\prime}(\rho,u)+A_{2}^{\prime}(\rho,u)(\nabla u),$

and in terms of the following additional quantity, which needs to be included
similarly as in (3.15) at finite Péclet number,

$B_{3}^{\prime}(\rho,u)\,:=\,(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta+u\cdot\nabla)(\rho
A_{1}(u)^{2})+\rho\big{(}(\nabla u)^{T}A_{1}(u)^{2}+A_{1}(u)^{2}(\nabla
u)\big{)}.$

We turn to the corresponding well-posedness question for the above
inhomogeneous models. As in the homogeneous setting, we focus for shortness on
the second-order model (3.32), which is again ill-posed whenever
$\alpha_{1}<0$. We consider the perturbative case of a weak nonlinearity,

$\alpha_{1}=\varepsilon\gamma_{1},\qquad\alpha_{2}=\varepsilon\gamma_{2},\qquad\varepsilon\ll
1,$

and we shall define a well-posed notion of approximate solutions that only
satisfy equation (3.32) up to higher-order $O(\varepsilon^{2})$ remainder:

$\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+\bar{u}_{\varepsilon}\cdot\nabla)\bar{u}_{\varepsilon}-{\operatorname{div}}(\bar{\sigma}_{\varepsilon})+\nabla\bar{p}_{\varepsilon}=h+O(\varepsilon^{2}),\\\\[2.84526pt]
(\partial_{t}+\bar{u}_{\varepsilon}\cdot\nabla)\bar{\rho}_{\varepsilon}-(\tfrac{1}{\operatorname{Pe}}+\varepsilon\mu_{0})\Delta\bar{\rho}_{\varepsilon}=O(\varepsilon^{2}),\\\\[2.84526pt]
\bar{\sigma}_{\varepsilon}\,=\,(\eta_{0}+\eta_{1}\bar{\rho}_{\varepsilon})A_{1}(\bar{u}_{\varepsilon})+\varepsilon\gamma_{1}A_{2}^{\prime}(\bar{\rho}_{\varepsilon},\bar{u}_{\varepsilon})+\varepsilon\gamma_{2}A_{1}(\bar{u}_{\varepsilon})^{2},\\\\[2.84526pt]
{\operatorname{div}}(\bar{u}_{\varepsilon})=0.\end{array}\right.$ (3.34)

For later purposes, note that we henceforth increase the diffusion of the
particle density by an additional constant $\mu_{0}\geq 0$, which will appear
in our setting as a possible effect of particle swimming velocities. Similarly
as in Proposition 3.2 for the homogeneous case, the simplest notion of well-
posed solutions takes the form of hierarchical solutions as described in the
following statement. The proof is postponed to Appendix B. Note that the
regularity theory for (3.35) below is quite delicate in the 3D Stokes case, as
a particularly careful stepwise argument is needed to first cover low-
regularity situations. The notion of Boussinesq-type solutions of Proposition
3.3 could also be easily extended to the present inhomogeneous setting, but we
skip the detail for conciseness.

######  Proposition 3.4 (Hierarchical solutions).

Consider equation (3.34) with parameters $\eta_{0},\operatorname{Pe}>0$,
$\eta_{1},\mu_{0}\geq 0$, and $\gamma_{1},\gamma_{2}\in\mathbb{R}$. If
$(u_{0},\rho_{0}),(u_{1},\rho_{1})$ are smooth solutions of the following two
auxiliary systems,

$\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+u_{0}\cdot\nabla)u_{0}-{\operatorname{div}}\big{(}2(\eta_{0}+\eta_{1}\rho_{0})\operatorname{D}(u_{0})\big{)}+\nabla
p_{0}=h,\\\\[2.84526pt]
(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta+u_{0}\cdot\nabla)\rho_{0}=0,\\\\[2.84526pt]
{\operatorname{div}}(u_{0})=0,\\\\[2.84526pt]
u_{0}|_{t=0}=u^{\circ}\quad\text{if $\operatorname{Re}\neq
0$},\quad\int_{\mathbb{T}^{d}}u_{0}=0\quad\text{if
$\operatorname{Re}=0$},\\\\[2.84526pt]
\rho_{0}|_{t=0}=\rho^{\circ},\end{array}\right.$ (3.35)
$\left\\{\begin{array}[]{l}\operatorname{Re}\big{(}(\partial_{t}+u_{0}\cdot\nabla)u_{1}+(u_{1}\cdot\nabla)u_{0}\big{)}-{\operatorname{div}}\big{(}2(\eta_{0}+\eta_{1}\rho_{0})\operatorname{D}(u_{1})\big{)}+\nabla
p_{1}\\\\[1.42262pt] \hskip
99.58464pt={\operatorname{div}}\big{(}2\eta_{1}\rho_{1}\operatorname{D}(u_{0})+\gamma_{1}A_{2}^{\prime}(\rho_{0},u_{0})+\gamma_{2}(2\operatorname{D}(u_{0}))^{2}\big{)},\\\\[2.84526pt]
(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta+u_{0}\cdot\nabla)\rho_{1}=\mu_{0}\Delta\rho_{0}-u_{1}\cdot\nabla\rho_{0},\\\\[2.84526pt]
{\operatorname{div}}(u_{1})=0,\\\\[2.84526pt] u_{1}|_{t=0}=0\quad\text{if
$\operatorname{Re}\neq 0$},\quad\int_{\mathbb{T}^{d}}u_{1}=0\quad\text{if
$\operatorname{Re}=0$},\\\\[2.84526pt] \rho_{1}|_{t=0}=0,\end{array}\right.$
(3.36)

then the superposition
$(\bar{u}_{\varepsilon},\bar{\rho}_{\varepsilon})=(u_{0}+\varepsilon
u_{1},\rho_{0}+\varepsilon\rho_{1})$ satisfies equation (3.34) with some
controlled remainder $O(\varepsilon^{2})$ and with initial condition
$\bar{\rho}_{\varepsilon}|_{t=0}=\rho^{\circ}$ (and
$\bar{u}_{\varepsilon}|_{t=0}=u^{\circ}$ if $\operatorname{Re}\neq 0$). For
the well-posedness of (3.35) and (3.36), we separately consider the Stokes and
Navier–Stokes cases:

  1. (i)

_Stokes case $\operatorname{Re}=0$, $d\leq 3$:_
Given
$\rho^{\circ}\in\operatorname{L}^{2}(\mathbb{T}^{d})\cap\frac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d})$
and
$h\in\operatorname{L}^{\infty}_{{\operatorname{loc}}}(\mathbb{R}^{+};H^{-1}(\mathbb{T}^{d})^{d})$,
there exists a unique global solution $(u_{0},\rho_{0})$ of (3.35) with

$\displaystyle u_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{1}(\mathbb{T}^{d})^{d}),$
$\displaystyle\rho_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}\big{(}\mathbb{R}^{+};\operatorname{L}^{2}(\mathbb{T}^{d})\cap\tfrac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d})\big{)}\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{1}(\mathbb{T}^{d})).$

Moreover, for all integers $s\geq\frac{d}{2}+1$, provided that
$\rho^{\circ}\in H^{s}(\mathbb{T}^{d})$ and that $h$ belongs to
$\operatorname{L}^{\infty}_{{\operatorname{loc}}}(\mathbb{R}^{+};H^{s-1}(\mathbb{R}^{+}\times\mathbb{T}^{d})^{d})$,
this solution further satisfies

$\displaystyle u_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})^{d}),$
$\displaystyle\rho_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{d}))\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})).$

In that case, if furthermore $h\in
W^{1,\infty}_{{\operatorname{loc}}}(\mathbb{R}^{+};H^{s-3}(\mathbb{T}^{d})^{d})$,
there exists a unique global solution $(u_{1},\rho_{1})$ of (3.36) with

$\displaystyle u_{1}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s-1}(\mathbb{T}^{d})^{d}),$
$\displaystyle\rho_{1}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s-2}(\mathbb{T}^{d}))\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s-1}(\mathbb{T}^{d})).$

  2. (ii)

_Navier–Stokes case $\operatorname{Re}=1$, $d=2$:_
Given
$\rho^{\circ}\in\operatorname{L}^{2}(\mathbb{T}^{2})\cap\frac{1}{2\pi}\mathcal{P}(\mathbb{T}^{2})$,
$u^{\circ}\in\operatorname{L}^{2}(\mathbb{T}^{2})^{2}$ which satisfies
${\operatorname{div}}(u^{\circ})=0$, and
$h\in\operatorname{L}^{2}_{{\operatorname{loc}}}(\mathbb{R}^{+};H^{-1}(\mathbb{T}^{2})^{2})$,
there exists a unique global solution $(u_{0},\rho_{0})$ of (3.35) with

$\displaystyle u_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};\operatorname{L}^{2}(\mathbb{T}^{2})^{2})\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{1}(\mathbb{T}^{2})^{2}),$
$\displaystyle\rho_{0}$ $\displaystyle\in$
$\displaystyle{\operatorname{L}^{\infty}_{\operatorname{loc}}\big{(}\mathbb{R}^{+};\operatorname{L}^{2}(\mathbb{T}^{2})\cap\tfrac{1}{2\pi}\mathcal{P}(\mathbb{T}^{2})\big{)}\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{1}(\mathbb{T}^{2})).}$

Moreover, for all integers $s\geq 2$, provided that $\rho^{\circ}\in
H^{s}(\mathbb{T}^{2})$, $u^{\circ}\in H^{s}(\mathbb{T}^{2})^{2}$, and
$h\in\operatorname{L}^{2}_{{\operatorname{loc}}}(\mathbb{R}^{+};H^{s-1}(\mathbb{T}^{2})^{2})$,
this solution further satisfies

$\displaystyle u_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{2})^{2})\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{2})^{2}),$
$\displaystyle\rho_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{2}))\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{2})).$

In that case, there exists a unique global solution $(u_{1},\rho_{1})$ of
(3.36) with

$\displaystyle u_{1}$ $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\in\leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s-2}(\mathbb{T}^{2})^{2})\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s-1}(\mathbb{T}^{2})^{2}),$
$\displaystyle\rho_{1}$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \in\leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s-2}(\mathbb{T}^{2}))\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s-1}(\mathbb{T}^{2})).\qed$

##  4\. Statement of main results

We turn to the precise statement of our main result, that is, the rigorous
hydrodynamic approximation of the Doi–Saintillan–Shelley theory in the
small-$\operatorname{Wi}$ regime. For simplicity, we focus on the first-order
approximation and the emergence of the second-order fluid model, but the same
analysis can be pursued to arbitrary order (see Section 4.2). More precisely,
we derive the second-order fluid model (3.34) with explicit coefficients

$\displaystyle\begin{split}\eta_{0}=1,\qquad\eta_{1}=\lambda\tfrac{(\theta+2)\omega_{d}}{2d(d+2)},\qquad\mu_{0}=\tfrac{1}{d(d-1)}U_{0}^{2},\\\\[2.84526pt]
\gamma_{1}=-\lambda\theta\tfrac{\omega_{d}}{4d^{2}(d+2)},\qquad\gamma_{2}=\lambda\tfrac{\omega_{d}}{2d^{2}(d+4)}(\theta+\tfrac{2d}{d+2}),\end{split}$
(4.1)

where we recall that $\lambda,\theta,U_{0}$ are parameters from the
Doi–Saintillan–Shelley system (2.5). Before formulating a precise result, we
introduce a suitable well-preparedness assumption for initial data. More
precisely, in order to avoid initial boundary layers due to the
$O(\frac{1}{\varepsilon})$ rotational diffusion, we first need to assume that
to leading order the initial density is invariant under this rotational
diffusion, which means that it is isotropic to leading order,

$f_{\varepsilon}(x,n)|_{t=0}\,=\,f_{\varepsilon}^{\circ}(x,n)\,=\,\rho^{\circ}(x)+O(\varepsilon),\qquad{\rho^{\circ}\in\tfrac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d})}.$

Yet, this is not sufficient: in order to avoid initial boundary layers in the
higher-order $\varepsilon$-expansion, we further need to assume that initial
data are compatible with the formal limiting hierarchy, which is precisely the
content of the assumption below. Similar issues are well known for higher-
order hydrodynamic expansions in the Boltzmann theory, see e.g. [Caf80, SK83,
Lac87].

######  Assumption 4.1 (Well-preparedness).

Let $h\in C^{\infty}(\mathbb{R}^{+}\times\mathbb{T}^{d})^{d}$, let
$\rho^{\circ}\in
H^{s}(\mathbb{T}^{d})\cap\frac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d})$ for
some $s\gg 1$, and let also $u^{\circ}\in H^{s}(\mathbb{T}^{d})^{d}$ with
${\operatorname{div}}(u^{\circ})=0$ in the Navier–Stokes case. We assume that
the initial condition $f_{\varepsilon}|_{t=0}=f_{\varepsilon}^{\circ}$ for the
Doi–Saintillan–Shelley system (2.5) is well-prepared in the following sense:
decomposing

$f_{\varepsilon}^{\circ}(x,n)\,=\,\rho_{\varepsilon}^{\circ}(x)+g_{\varepsilon}^{\circ}(x,n),\qquad\rho_{\varepsilon}^{\circ}\,:=\,\langle
f_{\varepsilon}^{\circ}\rangle\,=\,\fint_{\mathbb{S}^{d-1}}f_{\varepsilon}^{\circ}(\cdot,n)\,\mathrm{d}n\,\in\,\tfrac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d}),$

we have

$\rho_{\varepsilon}^{\circ}=\rho^{\circ}\qquad\text{and}\qquad\varepsilon^{\frac{1}{2}}\|g_{\varepsilon}^{\circ}-(\varepsilon
g_{1}+\varepsilon^{2}g_{2})|_{t=0}\|_{\operatorname{L}^{2}_{x,n}}\leq\,C_{0}\varepsilon^{3},$

for some constant $C_{0}<\infty$, where $g_{1}$ and $g_{2}$ are the solutions
of the hierarchical equations (5.3) and (5.4) below with initial data
$\rho^{\circ}$ in the Stokes case and $(u^{\circ},\rho^{\circ})$ in the
Navier–Stokes case. ∎

Note that this well-preparedness assumption is compatible with the positivity
$f_{\varepsilon}|_{t=0}\geq 0$ for $\varepsilon$ small enough, which is
necessary to ensure well-posedness of the Doi–Saintillan–Shelley system (2.5),
cf. Proposition 2.1. In these terms, we are now in position to state our main
result, thus finally providing a more detailed statement of Theorem 1.1. The
proof is given in Section 5.

######  Theorem 4.2 (Small-$\operatorname{Wi}$ expansion).

Let $h\in C^{\infty}(\mathbb{R}^{+}\times\mathbb{T}^{d})^{d}$, let
$\rho^{\circ}\in
H^{s}(\mathbb{T}^{d})\cap\frac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d})$ for
some $s\gg 1$, and let also $u^{\circ}\in H^{s}(\mathbb{T}^{d})^{d}$ with
${\operatorname{div}}(u^{\circ})=0$ in the Navier–Stokes case. Denote by
$(u_{\varepsilon},f_{\varepsilon})$ the global solution of the
Doi–Saintillan–Shelley model (2.5) as given by Proposition 2.1 with initial
condition
$f_{\varepsilon}|_{t=0}=f_{\varepsilon}^{\circ}\in\operatorname{L}^{2}\cap\mathcal{P}(\mathbb{T}^{d}\times\mathbb{S}^{d-1})$
satisfying the well-preparedness of Assumption 4.1. Further assume that

$\displaystyle\lambda\theta(1+\operatorname{Pe})\|\rho^{\circ}\|_{\operatorname{L}^{\infty}_{x}}\,\ll\,1$

is smaller than some universal constant.

  1. (i)

_Stokes case $\operatorname{Re}=0$, $d\leq 3$:_
Let $(\bar{u}_{\varepsilon},\bar{\rho}_{\varepsilon})$ be the unique global
hierarchical solution of (3.34) as given by Proposition 3.4(i) with explicit
coefficients (4.1). Then we have for all $t\geq 0$,

$\displaystyle\|\nabla(u_{\varepsilon}-\bar{u}_{\varepsilon})\|_{\operatorname{L}^{2}_{t}\operatorname{L}^{2}_{x}}$
$\displaystyle\lesssim$ $\displaystyle\varepsilon^{2},$
$\displaystyle\|\rho_{\varepsilon}-\bar{\rho}_{\varepsilon}\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x}}+\|\nabla(\rho_{\varepsilon}-\bar{\rho}_{\varepsilon})\|_{\operatorname{L}^{2}_{t}\operatorname{L}^{2}_{x}}$
$\displaystyle\lesssim$ $\displaystyle\varepsilon^{2},$

where multiplicative constants depend on $\operatorname{Pe}$ and on an upper
bound on $t,\lambda,U_{0},$ and on controlled norms of $h$ and $\rho^{\circ}$.

  2. (ii)

_Navier–Stokes case $\operatorname{Re}=0$, $d=2$:_
Let $(\bar{u}_{\varepsilon},\bar{\rho}_{\varepsilon})$ be the unique global
hierarchical solution of (3.34) as given by Proposition 3.4(ii) with explicit
coefficients (4.1). Then we have for all $t\geq 0$,

$\displaystyle\|u_{\varepsilon}-\bar{u}_{\varepsilon}\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x}}+\|\nabla(u_{\varepsilon}-\bar{u}_{\varepsilon})\|_{\operatorname{L}^{2}_{t}\operatorname{L}^{2}_{x}}$
$\displaystyle\leavevmode\nobreak\ \lesssim\leavevmode\nobreak\
\varepsilon^{2},$
$\displaystyle\|\rho_{\varepsilon}-\bar{\rho}_{\varepsilon}\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x}}+\|\nabla(\rho_{\varepsilon}-\bar{\rho}_{\varepsilon})\|_{\operatorname{L}^{2}_{t}\operatorname{L}^{2}_{x}}$
$\displaystyle\leavevmode\nobreak\ \lesssim\leavevmode\nobreak\
\varepsilon^{2},$

where multiplicative constants depend on $\operatorname{Pe}$ and on an upper
bound on $t,\lambda,U_{0},$ and on controlled norms of $h,\rho^{\circ},$ and
$u^{\circ}$. ∎

###  4.1. Non-Newtonian properties of hydrodynamic approximation

We comment on the rheological features of the obtained second-order fluid
system (3.34) with the explicit parameters
$\eta_{0},\eta_{1},\mu_{0},\gamma_{1},\gamma_{2}$ as defined in (4.1), briefly
describing the resulting non-Newtonian properties and how they depend on
microscopic features. We focus here on the physically relevant 3D setting, and
we point out that in the passive case the parameter values agree with [HL72,
Eqn (41) and Table 2] and [Bre74, Eqn (7.4)].

  1. $\bullet$

Effective spatial diffusion:
The spatial diffusion $\frac{1}{\operatorname{Pe}}$ of the suspended particle
density is enhanced by particles’ activity even at infinite Péclet number: it
is replaced by $\frac{1}{\operatorname{Pe}}+\varepsilon\mu_{0}$ with
$\mu_{0}=\frac{1}{6}U_{0}^{2}$. This naturally follows from the coupling of
particles’ swimming velocity with their rotational diffusion. This phenomenon
of increased mixing has been observed in studies such as [SS07, Sai18].

  2. $\bullet$

Modified zero-shear viscosity:
The presence of suspended particles leads to a non-trivial contribution to the
zero-shear viscosity: in the homogeneous setting
$\bar{\rho}_{\varepsilon}\equiv\frac{1}{\omega_{d}}$, we obtain a zero-shear
viscosity

$\displaystyle\tilde{\eta}_{0}\,:=\,\eta_{0}+\tfrac{1}{\omega_{d}}\eta_{1}=1+\tfrac{1}{30}\lambda(2+\theta).$

In particular, in case of passive particles, this zero-shear viscosity is
always larger than the plain fluid viscosity, $\tilde{\eta}_{0}>1$, while
particles’ activity can reverse this effect. For a precise description, first
recall that
$\theta=6+2\alpha\tfrac{|V_{0}|\ell\mu_{\operatorname{fl}}}{k_{B}\Theta}$, cf.
(2.4), where $\alpha$ characterizes the swimming mechanism:

    1. —

for passive particles $\alpha=0$, the zero-shear viscosity is
$\tilde{\eta}_{0}=1+\frac{4}{15}\lambda>1$;

    2. —

for so-called puller particles $\alpha>0$, the zero-shear viscosity is even
larger than for passive particles;

    3. —

for so-called pusher particles $\alpha<0$, the zero-shear viscosity is smaller
than for passive particles, and it can even be smaller than the plain fluid
viscosity provided that the activity of the particles is strong enough: we
find $\tilde{\eta}_{0}<1$ if
$\alpha<-\frac{4k_{B}\Theta}{|V_{0}|\ell\mu_{\operatorname{fl}}}$.

This prediction is consistent with well-known experimental results, see e.g.
[SA09, RJP10, Lóp+15], and it has been largely confirmed in the literature
[HABK08, Hai+09, Sai10a, Sai10, GLAB11, AMES16]. In particular, for E. coli
bacteria (a typical pusher particle), we can assess the value of the parameter
$\alpha$ using the experimental measurements performed in [Dre+11]: this
yields $\alpha<-\frac{4k_{B}\Theta}{|V_{0}|\ell\mu_{\operatorname{fl}}}$ and
the experimental findings of [Gac+13] then confirm our prediction that the
effective zero-shear viscosity is smaller than the plain fluid viscosity.

  3. $\bullet$

Normal-stress differences:
In the homogeneous setting
$\bar{\rho}_{\varepsilon}\equiv\frac{1}{\omega_{d}}$, we obtain the following
values for first and second-normal stress coefficients, as defined in (3.6),

$\nu_{10}=\tfrac{\theta}{90}\varepsilon\lambda\qquad\text{and}\qquad\nu_{20}=\tfrac{3-\theta}{315}\varepsilon\lambda.$

For passive particles ($\theta=6$), we thus obtain $\nu_{10}>0$, $\nu_{20}<0$,
and $\nu_{10}/|\nu_{20}|=7$, which agrees with experiments as discussed in
Section 3. The amplitude of these normal stress coefficients is even increased
in case of puller particles. In contrast, for pusher particles, normal stress
coefficients are reduced, and a very large activity could even result into
opposite effects: we find $\nu_{10}<0$ if $\theta<0$, and $\nu_{20}>0$ if
$\theta<3$. This behavior was also predicted in [Sai10, PRB16], but has yet to
be experimentally verified.

  4. $\bullet$

Elongational viscosity:
In the homogeneous setting
$\bar{\rho}_{\varepsilon}\equiv\frac{1}{\omega_{d}}$, in a uniaxial
elongational flow in the direction $e_{1}$, that is,
$\bar{u}_{\varepsilon}=\kappa(x_{1}e_{1}-\tfrac{1}{2}(x_{2}e_{2}+x_{3}e_{3}))$,
we obtain the following value for the elongational viscosity, as defined in
(3.10),

$\eta_{E}=\big{(}3+\tfrac{1}{10}\lambda(2+\theta)\big{)}+\kappa\tfrac{1}{140}\varepsilon\lambda(4+\theta).$

This shows that passive suspensions ($\theta=6$) lead to a strain-thickening
behavior, which is even increased in case of puller particles. In contrast,
for pusher particles, the strain-thickening behavior is reduced, and a very
large activity could even result in the opposite effect: the system becomes
strain-thinning if $\theta<-4$ (that is,
$\alpha<-\tfrac{5k_{B}\Theta}{|V_{0}|\ell\mu_{\operatorname{fl}}}$).

###  4.2. Next-order description

The above result is easily pursued to higher orders in the
small-$\operatorname{Wi}$ expansion. For shortness, we stick here to a formal
discussion. First, we show in Appendix C that the next-order description of
the suspended particle density involves additional nontrivial transport and
anisotropic diffusion terms depending on the surrounding fluid flow: we find

$(\partial_{t}+\bar{u}_{\varepsilon}\cdot\nabla)\bar{\rho}_{\varepsilon}-{\operatorname{div}}\Big{(}\big{(}{\tfrac{1}{\operatorname{Pe}}}+\varepsilon\mu_{0}+\varepsilon^{2}\mu_{1}\operatorname{D}(\bar{u}_{\varepsilon})\big{)}\nabla\bar{\rho}_{\varepsilon}\Big{)}=\varepsilon^{2}\mu_{2}\,{\operatorname{div}}(\bar{\rho}_{\varepsilon}\Delta\bar{u}_{\varepsilon})+O(\varepsilon^{3}),$
(4.2)

with explicit coefficients

$\displaystyle\mu_{0}$ $\displaystyle:=$
$\displaystyle\tfrac{U_{0}^{2}}{d(d-1)},$ $\displaystyle\mu_{1}$
$\displaystyle:=$ $\displaystyle\tfrac{(3d+1)U_{0}^{2}}{d(d-1)^{2}(d+2)},$
$\displaystyle\mu_{2}$ $\displaystyle:=$
$\displaystyle\tfrac{U_{0}^{2}}{2d(d-1)(d+2)}.$

In particular, this shows that homogeneous spatial densities are still stable
to order $O(\varepsilon^{3})$, and we shall henceforth restrict for simplicity
to the homogeneous setting,

$\bar{\rho}_{\varepsilon}=\tfrac{1}{\omega_{d}}+O(\varepsilon^{3}).$

In addition, we shall focus on the case of infinite Péclet number and of
vanishing particle swimming velocity,

$\operatorname{Pe}=\infty,\qquad U_{0}=0,$

as this choice strongly simplifies the macroscopic equations and as it seems
anyhow to be the most relevant setting physically, cf. Section 2 (recall
however that our rigorous results do not hold for $\operatorname{Pe}=\infty$).
In this setting, we formally derive in Appendix C the following third-order
fluid equations,

$\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+\bar{u}_{\varepsilon}\cdot\nabla)\bar{u}_{\varepsilon}-{\operatorname{div}}(\bar{\sigma}_{\varepsilon})+\nabla\bar{p}_{\varepsilon}=h+O(\varepsilon^{3}),\\\\[2.84526pt]
{\operatorname{div}}(\bar{u}_{\varepsilon})=0,\end{array}\right.$ (4.3)

where the stress is given by

$\bar{\sigma}_{\varepsilon}=(1+\eta_{1})A_{1}(\bar{u}_{\varepsilon})+\varepsilon\gamma_{1}A_{2}(\bar{u}_{\varepsilon})+\varepsilon\gamma_{2}A_{1}(\bar{u}_{\varepsilon})^{2}\\\
\qquad\,+\varepsilon^{2}\kappa_{1}A_{3}(\bar{u}_{\varepsilon})+\varepsilon^{2}\kappa_{2}\big{(}A_{1}(\bar{u}_{\varepsilon})A_{2}(\bar{u}_{\varepsilon})+A_{2}(\bar{u}_{\varepsilon})A_{1}(\bar{u}_{\varepsilon})\big{)}+\varepsilon^{2}\kappa_{3}A_{1}(\bar{u}_{\varepsilon}){\operatorname{tr}}(A_{1}(\bar{u}_{\varepsilon})^{2})$

with explicit coefficients

$\displaystyle\eta_{1}$ $\displaystyle:=$
$\displaystyle\lambda\tfrac{1}{2d(d+2)}(\theta+2),$ $\displaystyle\gamma_{1}$
$\displaystyle:=$ $\displaystyle-\lambda\theta\tfrac{1}{4d^{2}(d+2)},$
$\displaystyle\gamma_{2}$ $\displaystyle:=$
$\displaystyle\lambda\tfrac{1}{2d^{2}(d+4)}(\theta+\tfrac{2d}{d+2}),$
$\displaystyle\kappa_{1}$ $\displaystyle:=$
$\displaystyle\lambda\theta\tfrac{1}{8d^{3}(d+2)},$ $\displaystyle\kappa_{2}$
$\displaystyle:=$
$\displaystyle-\lambda\tfrac{1}{8d^{3}(d+4)}(3\theta+\tfrac{2d}{d+2}),$
$\displaystyle\kappa_{3}$ $\displaystyle:=$
$\displaystyle\lambda\tfrac{1}{8d^{3}(d+2)^{2}(d+4)(d+6)}\Big{(}2d(3d^{2}+10d+6)+\theta(d+4)(3d^{2}+11d+12)\Big{)}.$

These third-order fluid coefficients coincide in the passive case ($\theta=6$)
with those computed by Brenner [Bre74, Eq. (7.4)] (once the notation is
properly identified). Regarding the non-Newtonian phenomena discussed in
Section 3, the main observation is that this third-order fluid model describes
the expected shear-thinning behavior of the suspension. Indeed, the shear-
dependent viscosity is given in 3D as follows, cf. Section 3.1,

$\displaystyle\kappa\leavevmode\nobreak\ \leavevmode\nobreak\
\mapsto\leavevmode\nobreak\ \leavevmode\nobreak\ $ $\displaystyle
1+\eta_{1}+2\varepsilon^{2}(\kappa_{2}+\kappa_{3})\kappa^{2}$
$\displaystyle\,\qquad\qquad=\leavevmode\nobreak\ $ $\displaystyle
1+\lambda\tfrac{\theta+2}{30}-\varepsilon^{2}\lambda\tfrac{19\theta-12}{18900}\kappa^{2},$

which is decreasing in $\kappa$ if and only if $\theta>\frac{12}{19}$. As
expected, this shows that passive suspensions ($\theta=6$) lead to a shear-
thinning behavior, which is even increased in case of puller particles. In
contrast, for pusher particles, the shear-thinning behavior is reduced, and a
very large activity could even result in the opposite effect: the system
becomes shear-thickening if $\theta<\frac{12}{19}$ (that is,
$\alpha<-\frac{51k_{B}\Theta}{19|V_{0}|\ell\mu_{\operatorname{fl}}}$). This
possible shear-thickening effect was indeed measured experimentally in
[Gac+13, Lóp+15] for suspensions of E. coli bacteria (pusher-type
microswimmers) with strong enough activity. We also refer to [Hai+09, GLAB11,
PRB16] for analytical and numerical results showing the same effect.

We note that for $U_{0}\neq 0$ the corresponding fluid equation for
$\bar{u}_{\varepsilon}$ would differ from the 3rd-order fluid model even at
infinite Péclet number and at uniform particle density. In particular, an
additional dispersive correction
$-\varepsilon^{2}\kappa_{4}\Delta^{2}\bar{u}_{\varepsilon}$ needs to be
included in the fluid equation. We skip the detail as the case $|U_{0}|\ll 1$
seems to be the most relevant physically.

##  5\. Small-$\operatorname{Wi}$ expansion of Doi–Saintillan–Shelley theory

This section is devoted to the small-$\varepsilon$ expansion of the solution
$(u_{\varepsilon},f_{\varepsilon})$ of the Doi–Saintillan–Shelley model (2.5),
as well as to the identification of second-order fluid equations satisfied by
the truncated expansion. We naturally split the particle density as

$f_{\varepsilon}(x,n)\,=\,\rho_{\varepsilon}(x)+g_{\varepsilon}(x,n),\qquad\rho_{\varepsilon}\,:=\,\langle
f_{\varepsilon}\rangle\,:=\,\fint_{\mathbb{S}^{d-1}}f_{\varepsilon}(\cdot,n)\,\mathrm{d}n\,\in\,\tfrac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d}),$

where $\rho_{\varepsilon}$ stands for the spatial density and where $\langle
g_{\varepsilon}\rangle=0$. Recall the well-preparedness condition of
Assumption 4.1: we assume in particular

$\rho_{\varepsilon}|_{t=0}=\rho^{\circ},\qquad
g_{\varepsilon}|_{t=0}=O(\varepsilon),\qquad\rho^{\circ}\in\tfrac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d}).$
(5.1)

meaning that we start from an initial density that is to leading order at
equilibrium with respect to the strong rotational diffusion in (2.5). In this
setting, we shall analyze the asymptotic behavior of the solution
$(u_{\varepsilon},f_{\varepsilon})$ and derive a hydrodynamic approximation in
the spirit of Hilbert’s expansion method in the Boltzmann theory [Hil12,
Caf80, Gol05, SR09]. We start from the ansatz

$\displaystyle\begin{split}u_{\varepsilon}&=u_{0}+\varepsilon
u_{1}+\varepsilon^{2}u_{2}+\ldots,\\\
\rho_{\varepsilon}&=\rho_{0}+\varepsilon\rho_{1}+\varepsilon^{2}\rho_{2}+\ldots,\\\
g_{\varepsilon}&=g_{0}+\varepsilon
g_{1}+\varepsilon^{2}g_{2}+\varepsilon^{3}g_{3}+\ldots,\end{split}$ (5.2)

with $\rho_{0}\in\frac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d})$, with
$\int_{\mathbb{T}^{d}}\rho_{n}=0$ for $n\geq 1$, and with $\langle
g_{n}\rangle=0$ for all $n\geq 0$. Inserting it into the system (2.5), and
identifying powers of $\varepsilon$, we are led formally to the following
hierarchy of coupled equations:

  1. $\bullet$

_Order $O(\varepsilon^{-1})$:_ The equation for the particle density yields
$\Delta_{n}g_{0}=0$, and therefore, as by definition $\langle g_{0}\rangle=0$,
we must have

$g_{0}=0.$

On the other hand, the fluid equation yields
${\operatorname{div}}_{x}(\sigma_{1}[\rho_{0}+g_{0}])=0$, which is then
automatically satisfied as $g_{0}=0$ and as $\rho_{0}$ does not depend on $n$.

  2. $\bullet$

_Order $O(\varepsilon^{0})$:_ The triplet $(u_{0},\rho_{0},g_{1})$ satisfies

$\left\\{\begin{array}[]{l}\operatorname{Re}(\partial_{t}+u_{0}\cdot\nabla)u_{0}-\Delta
u_{0}+\nabla
p_{0}=h+{\operatorname{div}}(\sigma_{1}[g_{1}])+{\operatorname{div}}(\sigma_{2}[\rho_{0},\nabla
u_{0}]),\\\\[2.84526pt]
\Delta_{n}g_{1}=U_{0}n\cdot\nabla_{x}\rho_{0}+{\operatorname{div}}_{n}\big{(}\pi_{n}^{\perp}(\nabla
u_{0})n\rho_{0}\big{)},\\\\[2.84526pt]
(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta+u_{0}\cdot\nabla)\rho_{0}=0,\\\\[2.84526pt]
{\operatorname{div}}(u_{0})=0,\leavevmode\nobreak\ \leavevmode\nobreak\
\langle g_{1}\rangle=0,\\\\[2.84526pt]
u_{0}|_{t=0}=u^{\circ}\leavevmode\nobreak\ \leavevmode\nobreak\ \text{if
$\operatorname{Re}\neq 0$},\ \
\int_{\mathbb{T}^{d}}u_{0}=0\leavevmode\nobreak\ \leavevmode\nobreak\ \text{if
$\operatorname{Re}=0$},\\\\[2.84526pt]
\rho_{0}|_{t=0}=\rho^{\circ}.\end{array}\right.$ (5.3)

  3. $\bullet$

_Order $O(\varepsilon^{1})$:_ The triplet $(u_{1},\rho_{1},g_{2})$ satisfies

$\left\\{\begin{array}[]{l}\operatorname{Re}\big{(}(\partial_{t}+u_{0}\cdot\nabla)u_{1}+(u_{1}\cdot\nabla)u_{0}\big{)}-\Delta
u_{1}+\nabla p_{1}\\\\[2.84526pt] \hskip
56.9055pt\,=\,{\operatorname{div}}(\sigma_{1}[g_{2}])+{\operatorname{div}}\big{(}\sigma_{2}[\rho_{0},\nabla
u_{1}]+\sigma_{2}[\rho_{1}+g_{1},\nabla u_{0}]\big{)},\\\\[2.84526pt]
\Delta_{n}g_{2}=(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta_{x}+u_{0}\cdot\nabla_{x})g_{1}+P_{1}^{\bot}\big{(}U_{0}n\cdot\nabla_{x}(\rho_{1}+g_{1})\big{)}\\\\[2.84526pt]
\hskip 113.81102pt+{\operatorname{div}}_{n}\big{(}\pi_{n}^{\bot}(\nabla
u_{1})n\rho_{0}+\pi_{n}^{\bot}(\nabla
u_{0})n(\rho_{1}+g_{1})\big{)},\\\\[2.84526pt]
(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta_{x}+u_{0}\cdot\nabla_{x})\rho_{1}+u_{1}\cdot\nabla_{x}\rho_{0}+\langle
U_{0}n\cdot\nabla_{x}g_{1}\rangle=0,\\\\[2.84526pt]
{\operatorname{div}}(u_{1})=0,\leavevmode\nobreak\ \leavevmode\nobreak\
\langle g_{2}\rangle=0,\\\\[2.84526pt] u_{1}|_{t=0}=0\leavevmode\nobreak\
\leavevmode\nobreak\ \text{if $\operatorname{Re}\neq 0$},\ \
\int_{\mathbb{T}^{d}}u_{1}=0\leavevmode\nobreak\ \leavevmode\nobreak\ \text{if
$\operatorname{Re}=0$},\\\\[2.84526pt] \rho_{1}|_{t=0}=0.\end{array}\right.$
(5.4)

  4. $\bullet$

_Order $O(\varepsilon^{2})$:_ The triplet $(u_{2},\rho_{2},f_{3})$ satisfies

$\left\\{\begin{array}[]{l}\operatorname{Re}\big{(}(\partial_{t}+u_{0}\cdot\nabla)u_{2}+(u_{1}\cdot\nabla)u_{1}+(u_{2}\cdot\nabla)u_{0}\big{)}-\Delta
u_{2}+\nabla p_{2}\\\\[2.84526pt] \hskip
28.45274pt=\,{\operatorname{div}}(\sigma_{1}[g_{3}])+{\operatorname{div}}\big{(}\sigma_{2}[\rho_{0},\nabla
u_{2}]+\sigma_{2}[\rho_{1}+g_{1},\nabla
u_{1}]+\sigma_{2}[\rho_{2}+g_{2},\nabla u_{0}]\big{)},\\\\[2.84526pt]
\Delta_{n}g_{3}=(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta_{x}+u_{0}\cdot\nabla_{x})g_{2}+u_{1}\cdot\nabla_{x}g_{1}+P_{1}^{\bot}\big{(}U_{0}n\cdot\nabla_{x}(\rho_{2}+g_{2})\big{)}\\\\[2.84526pt]
\hskip 42.67912pt\,+\,{\operatorname{div}}_{n}\big{(}\pi_{n}^{\bot}(\nabla
u_{2})n\rho_{0}+\pi_{n}^{\bot}(\nabla
u_{1})n(\rho_{1}+g_{1})+\pi_{n}^{\bot}(\nabla
u_{0})n(\rho_{2}+g_{2})\big{)},\\\\[2.84526pt]
(\partial_{t}-\tfrac{1}{\operatorname{Pe}}\Delta_{x}+u_{0}\cdot\nabla_{x})\rho_{2}+u_{1}\cdot\nabla_{x}\rho_{1}+u_{2}\cdot\nabla_{x}\rho_{0}+\langle
U_{0}n\cdot\nabla_{x}g_{2}\rangle=0,\\\\[2.84526pt]
{\operatorname{div}}(u_{2})=0,\leavevmode\nobreak\ \leavevmode\nobreak\
\langle g_{3}\rangle=0,\\\\[2.84526pt] u_{2}|_{t=0}=0\leavevmode\nobreak\
\leavevmode\nobreak\ \text{if $\operatorname{Re}\neq 0$},\ \
\int_{\mathbb{T}^{d}}u_{2}=0\leavevmode\nobreak\ \leavevmode\nobreak\ \text{if
$\operatorname{Re}=0$},\\\\[2.84526pt] \rho_{2}|_{t=0}=0.\end{array}\right.$
(5.5)

In view of this hierarchy, we understand that the condition (5.1) for initial
data needs to be further strengthened to avoid initial boundary layers: more
precisely, we need to assume that the initial condition
$f_{\varepsilon}|_{t=0}=f_{\varepsilon}^{\circ}$ is compatible with the above
hierarchy, meaning that $g_{\varepsilon}$ coincides initially with
$\varepsilon g_{1}+\varepsilon^{2}g_{2}+\varepsilon^{3}g_{3}$ up to higher
order errors. To accuracy $O(\varepsilon^{2})$, this well-preparedness is
precisely the content of Assumption 4.1.

The well-posedness and the propagation of regularity for the above hierarchy
are stated in the following result. We emphasize that we were not able to find
references for the types of systems that naturally appear in this
$\varepsilon$-expansion, and we believe that some of these new well-posedness
results may be of independent interest (see for instance the system (5.18)
below). Note that the hierarchy is triangular as we can eliminate the
densities $g_{1},g_{2},g_{3}$ in terms of the velocity fields
$u_{0},u_{1},u_{2}$ and of the spatial densities $\rho_{0},\rho_{1},\rho_{2}$.
We focus on the Stokes case $\operatorname{Re}=0$, while the 2D Navier–Stokes
case follows up to straightforward adaptations and is omitted for shortness.
The proof is displayed in Section 5.2.

######  Proposition 5.1 (Well-posedness of hierarchy).

Consider the Stokes case $\operatorname{Re}=0$, $d\leq 3$.

  1. (i)

_Well-posedness for $(u_{0},\rho_{0},g_{1})$:_
Given integer $s\geq\tfrac{d}{2}+1$, $\rho^{\circ}\in
H^{s}(\mathbb{T}^{d})\cap\frac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d})$, and
$h\in\operatorname{L}^{\infty}_{{\operatorname{loc}}}(\mathbb{R}^{+};H^{s-1}(\mathbb{T}^{d})^{d})$,
there exists a unique solution $(u_{0},\rho_{0},g_{1})$ of the Cauchy problem
(5.3) with

$\displaystyle u_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})^{d}),$
$\displaystyle\rho_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}\big{(}\mathbb{R}^{+};H^{s}(\mathbb{T}^{d})\cap\tfrac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d})\big{)}\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})),$
$\displaystyle g_{1}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s-1}(\mathbb{T}^{d}\times\mathbb{S}^{d-1})),$

and $g_{1}$ is given by the explicit formula

$\displaystyle
g_{1}(\cdot,n)=-\tfrac{1}{d-1}U_{0}n\cdot\nabla\rho_{0}+\tfrac{1}{2}(n\otimes
n):\rho_{0}\operatorname{D}(u_{0}).$ (5.6)

Moreover, for all $r\geq 0$, provided that $h\in
W^{2,\infty}(\mathbb{R}^{+};H^{s-1}(\mathbb{T}^{d})^{d})$ and that $s$ is
chosen large enough, we also have

$\displaystyle\partial_{t}u_{0},\,\partial_{t}^{2}u_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{{\operatorname{loc}}}(\mathbb{R}^{+};H^{r}(\mathbb{T}^{d})^{d}),$
$\displaystyle\partial_{t}\rho_{0},\,\partial_{t}^{2}\rho_{0}$
$\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{{\operatorname{loc}}}(\mathbb{R}^{+};H^{r}(\mathbb{T}^{d})).$

  2. (ii)

_Well-posedness for $(u_{1},\rho_{1},g_{2})$:_
Given $s\geq 0$, provided that the solution $(u_{0},\rho_{0})$ of item (i) is
such that

$\displaystyle u_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+3}(\mathbb{T}^{d})^{d})\cap
W^{1,\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})^{d}),$
$\displaystyle\rho_{0}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+2}(\mathbb{T}^{d}))\cap
W^{1,\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{d})),$

there exists a unique solution $(u_{1},\rho_{1},g_{2})$ of the Cauchy problem
(5.4) with

$\displaystyle u_{1}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})^{d}),$
$\displaystyle\rho_{1}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{d}))\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})),$
$\displaystyle g_{2}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{d}\times\mathbb{S}^{d-1})),$

and $g_{2}$ is given by the explicit formula

$\displaystyle\qquad g_{2}(\cdot,n)$ $\displaystyle=$
$\displaystyle-\tfrac{1}{d-1}U_{0}n\cdot\nabla\rho_{1}-\tfrac{1}{2d}(n\otimes
n-\tfrac{1}{d}\operatorname{Id}):\Big{(}\tfrac{1}{4}A_{2}^{\prime}(u_{0},\rho_{0})-\rho_{0}\operatorname{D}(u_{0})^{2}$
(5.7) $\displaystyle\hskip
113.81102pt-d\rho_{0}\operatorname{D}(u_{1})-d\rho_{1}\operatorname{D}(u_{0})-\tfrac{1}{d-1}U_{0}^{2}\nabla^{2}\rho_{0}\Big{)}$
$\displaystyle-\tfrac{1}{3(d-1)}U_{0}(\nabla\rho_{0})\cdot\Big{(}n\big{(}(n\otimes
n):\operatorname{D}(u_{0})\big{)}+\tfrac{4}{d-1}\operatorname{D}(u_{0})n\Big{)}$
$\displaystyle-\tfrac{1}{6(d+1)}U_{0}\rho_{0}\,{\operatorname{div}}_{x}\Big{(}n\big{(}(n\otimes
n):\operatorname{D}(u_{0})\big{)}+\tfrac{4}{d-1}\operatorname{D}(u_{0})n\Big{)}$
$\displaystyle+\tfrac{1}{8}\rho_{0}\Big{(}\big{(}(n\otimes
n):\operatorname{D}(u_{0})\big{)}^{2}-\tfrac{2}{d(d+2)}\operatorname{tr}(\operatorname{D}(u_{0})^{2})\Big{)},$

in terms of the (non-standard) Rivlin–Ericksen tensor $A_{2}^{\prime}$ defined
in (3.33). Moreover, for all $r\geq 0$, provided that $s$ is chosen large
enough, we also have

$\displaystyle\partial_{t}u_{1}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{{\operatorname{loc}}}(\mathbb{R}^{+};H^{r}(\mathbb{T}^{d})^{d}),$
$\displaystyle\partial_{t}\rho_{1}$ $\displaystyle\in$
$\displaystyle\operatorname{L}^{\infty}_{{\operatorname{loc}}}(\mathbb{R}^{+};H^{r}(\mathbb{T}^{d})).$

  3. (iii)

_Well-posedness for $(u_{2},\rho_{2},g_{3})$:_
Given $s\geq 0$, provided that the solutions $(u_{0},\rho_{0})$ and
$(u_{1},\rho_{1})$ of items (i) and (ii) are such that

$\displaystyle(u_{0},\rho_{0})$ $\displaystyle\in$ $\displaystyle
W^{2,\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{r}(\mathbb{T}^{d})^{d+1}),$
$\displaystyle(u_{1},\rho_{1})$ $\displaystyle\in$ $\displaystyle
W^{1,\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{r}(\mathbb{T}^{d})^{d+1}),$

for some $r$ large enough, then there exists a unique solution
$(u_{2},\rho_{2},g_{3})$ of the Cauchy problem (5.5) with

$\displaystyle u_{2}$ $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\in\leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})^{d}),$
$\displaystyle\rho_{2}$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \in\leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname{L}^{\infty}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s}(\mathbb{T}^{d}))\cap\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s+1}(\mathbb{T}^{d})),$
$\displaystyle g_{3}$ $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\in\leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname{L}^{2}_{\operatorname{loc}}(\mathbb{R}^{+};H^{s-1}(\mathbb{T}^{d}\times\mathbb{S}^{d-1})).$

Associated to all the above well-posedness results are estimates of the
corresponding norms of the solutions in terms of all the parameters and of the
controlled norms of the data. ∎

With the above construction of the hierarchy
$\\{u_{n},\rho_{n},g_{n+1}\\}_{n\geq 0}$, we can now turn to the justification
of the formal expansion (5.2). We stick to order $O(\varepsilon^{3})$ for
conciseness, but the proof could be pursued to arbitrary order. The proof is
displayed in Section 5.3. Note that the well-preparedness assumption (5.8)
below for initial data is one order stronger than in Assumption 4.1: indeed,
while our main result in Section 4 focusses on $O(\varepsilon)$ effects, only
deriving second-order fluid models, the present result further describes
$O(\varepsilon^{2})$ effects and therefore requires this strengthened well-
preparedness condition. Although not needed for the purposes of our main
result, the present next-order analysis is included to illustrate how the
$\varepsilon$-expansion can be pursued to arbitrary order without additional
mathematical difficulties, then leading to higher-order fluid models; we refer
to Appendix C for the corresponding derivation of third-order fluid models.

######  Proposition 5.2 (Error estimates for $\varepsilon$-expansion).

Let $h\in C^{\infty}(\mathbb{R}^{+}\times\mathbb{T}^{d})^{d}$, let
$\rho^{\circ}\in
H^{s}(\mathbb{T}^{d})\cap\frac{1}{\omega_{d}}\mathcal{P}(\mathbb{T}^{d})$ for
some $s\gg 1$, and let also $u^{\circ}\in H^{s}(\mathbb{T}^{d})^{d}$ with
${\operatorname{div}}(u^{\circ})=0$ in the Navier–Stokes case. Denote by
$(u_{\varepsilon},f_{\varepsilon})$ the solution of the Doi–Saintillan–Shelley
model (2.5) as given by Proposition 2.1, and assume that the initial condition
$f_{\varepsilon}|_{t=0}=f_{\varepsilon}^{\circ}$ is well-prepared in the
following sense: decomposing
$f_{\varepsilon}^{\circ}=\rho_{\varepsilon}^{\circ}+g_{\varepsilon}^{\circ}$
with $\rho_{\varepsilon}^{\circ}:=\langle f_{\varepsilon}^{\circ}\rangle$, we
have in terms of the functions $u_{0},u_{1},u_{2},g_{1},g_{2},g_{3}$ defined
in Proposition 5.1 with data $(h,u^{\circ},\rho^{\circ})$,

${\rho_{\varepsilon}^{\circ}\,=\,\rho^{\circ}}\qquad\text{and}\qquad\varepsilon^{\frac{1}{2}}\|g_{\varepsilon}^{\circ}-(\varepsilon
g_{1}+\varepsilon^{2}g_{2}+\varepsilon^{3}g_{3})|_{t=0}\|_{\operatorname{L}^{2}_{x,n}}\,\leq\,C_{0}\varepsilon^{4},$
(5.8)

for some constant $C_{0}<\infty$. Further assume that

$\displaystyle\lambda\theta(1+\operatorname{Pe})\|\rho^{\circ}\|_{\operatorname{L}^{\infty}_{x}}\,\ll\,1$

is smaller than some universal constant.

  1. (i)

_Stokes case $\operatorname{Re}=0$, $d\leq 3$:_
For all $t\geq 0$, we have

$\displaystyle\hskip
28.45274pt\varepsilon^{\frac{1}{2}}\|\nabla(u_{\varepsilon}-u_{0}-\varepsilon
u_{1}-\varepsilon^{2}u_{2})\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x}}+\|\nabla(u_{\varepsilon}-u_{0}-\varepsilon
u_{1}-\varepsilon^{2}u_{2})\|_{\operatorname{L}^{2}_{t}\operatorname{L}^{2}_{x}}$
$\displaystyle\\!\\!\leq\\!\\!$ $\displaystyle\mathcal{C}(t)\varepsilon^{3},$
$\displaystyle\hskip
28.45274pt\varepsilon^{\frac{1}{2}}\|g_{\varepsilon}-\varepsilon
g_{1}-\varepsilon^{2}g_{2}-\varepsilon^{3}g_{3}\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x,n}}\\!+\|\nabla_{n}(g_{\varepsilon}-\varepsilon
g_{1}-\varepsilon^{2}g_{2}-\varepsilon^{3}g_{3})\|_{\operatorname{L}^{2}_{t}\operatorname{L}^{2}_{x,n}}$
$\displaystyle\\!\\!\leq\\!\\!$ $\displaystyle\mathcal{C}(t)\varepsilon^{4},$
$\displaystyle\hskip
28.45274pt\|\rho_{\varepsilon}-\rho_{0}-\varepsilon\rho_{1}-\varepsilon^{2}\rho_{2}\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x}}+\|\nabla(\rho_{\varepsilon}-\rho_{0}-\varepsilon\rho_{1}-\varepsilon^{2}\rho_{2})\|_{\operatorname{L}^{2}_{t}\operatorname{L}^{2}_{x}}$
$\displaystyle\\!\\!\leq\\!\\!$ $\displaystyle\mathcal{C}(t)\varepsilon^{3},$

provided that $\varepsilon\mathcal{C}(t)\ll 1$ is small enough, where the
multiplicative constant $\mathcal{C}(t)$ depends on $\operatorname{Pe}$ and on
an upper bound on

$\displaystyle\qquad t,\leavevmode\nobreak\ {C_{0},\leavevmode\nobreak\
\lambda,\leavevmode\nobreak\ U_{0},}\quad\|(\nabla u_{0},\nabla u_{1},\nabla
u_{2})\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{\infty}_{x}},\quad\|(\rho_{0},\rho_{1},\rho_{2})\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{\infty}_{x}},$
$\displaystyle\qquad\|(g_{1},g_{2},g_{3})\|_{\operatorname{L}^{\infty}_{t}W^{1,\infty}_{x,n}},\quad\|(\partial_{t}-\Delta_{x})g_{3}\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x,n}}.$

  2. (ii)

_Navier–Stokes case $\operatorname{Re}=0$, $d=2$:_
For all $t\geq 0$, we have

$\displaystyle\hskip 28.45274pt\|u_{\varepsilon}-u_{0}-\varepsilon
u_{1}-\varepsilon^{2}u_{2}\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x}}+\|\nabla(u_{\varepsilon}-u_{0}-\varepsilon
u_{1}-\varepsilon^{2}u_{2})\|_{\operatorname{L}^{2}_{t}\operatorname{L}^{2}_{x}}$
$\displaystyle\\!\\!\leq\\!\\!$ $\displaystyle\mathcal{C}(t)\varepsilon^{3},$
$\displaystyle\hskip
28.45274pt\varepsilon^{\frac{1}{2}}\|g_{\varepsilon}-\varepsilon
g_{1}-\varepsilon^{2}g_{2}-\varepsilon^{3}g_{3}\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x,n}}\\!+\|\nabla_{n}(g_{\varepsilon}-\varepsilon
g_{1}-\varepsilon^{2}g_{2}-\varepsilon^{3}g_{3})\|_{\operatorname{L}^{2}_{t}\operatorname{L}^{2}_{x,n}}$
$\displaystyle\\!\\!\leq\\!\\!$ $\displaystyle\mathcal{C}(t)\varepsilon^{4},$
$\displaystyle\hskip
28.45274pt\|\rho_{\varepsilon}-\rho_{0}-\varepsilon\rho_{1}-\varepsilon^{2}\rho_{2}\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x}}+\|\nabla(\rho_{\varepsilon}-\rho_{0}-\varepsilon\rho_{1}-\varepsilon^{2}\rho_{2})\|_{\operatorname{L}^{2}_{t}\operatorname{L}^{2}_{x}}$
$\displaystyle\\!\\!\leq\\!\\!$ $\displaystyle\mathcal{C}(t)\varepsilon^{3},$

provided that $\varepsilon\mathcal{C}(t)\ll 1$ is small enough, where
$\mathcal{C}(t)$ now depends $\operatorname{Pe}$ and on an upper bound on

$\displaystyle\qquad t,\leavevmode\nobreak\ {C_{0},\leavevmode\nobreak\
\lambda,\leavevmode\nobreak\
U_{0},}\quad\|(u_{0},u_{1},u_{2})\|_{\operatorname{L}^{\infty}_{t}W^{1,\infty}_{x}},\quad\|(\rho_{0},\rho_{1},\rho_{2})\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{\infty}_{x}},$
$\displaystyle\qquad\|(g_{1},g_{2},g_{3})\|_{\operatorname{L}^{\infty}_{t}W^{1,\infty}_{x,n}},\quad\|(\partial_{t}-\Delta_{x})g_{3}\|_{\operatorname{L}^{\infty}_{t}\operatorname{L}^{2}_{x,n}}.\qed$

Finally, it remains to identify the equation satisfied by the truncated
$\varepsilon$-expansion
$(\bar{u}_{\varepsilon},\bar{\rho}_{\varepsilon}):=(u_{0}+\varepsilon
u_{1},\rho_{0}+\varepsilon\rho_{1})$: we show that this expansion precisely
coincides with the hierarchical solution of the second-order fluid model
(3.34) with some explicit choice of parameters. The proof is displayed in
Section 5.4.

######  Proposition 5.3 (From hierarchy to second-order fluids).

Given the solutions $(u_{0},\rho_{0},g_{1})$ and $(u_{1},\rho_{1},g_{2})$ of
the hierarchy (5.3)–(5.4), as constructed in Proposition 5.1 in the Stokes
case, the superposition
$(\bar{u}_{\varepsilon},\bar{\rho}_{\varepsilon}):=(u_{0}+\varepsilon
u_{1},\rho_{0}+\varepsilon\rho_{1})$ coincides with the unique hierarchical
solution of the second-order fluid model (3.34) in the sense of Proposition
3.2, with coefficients $\eta_{0},\eta_{1},\mu_{0},\gamma_{1},\gamma_{2}$
explicitly given by (4.1). ∎

The combination of Propositions 5.2 and 5.3 completes the proof of Theorem
4.2. Note however that we only appeal to Assumption 4.1 in the statement of
Theorem 4.2, which is one order weaker than the well-preparedness assumption
(5.8) required in Proposition 5.2. Indeed, as explained, we focus in Theorem
4.2 on $O(\varepsilon)$ effects, only deriving the second-order fluid model,
while in Proposition 5.2 we took care to further describe $O(\varepsilon^{2})$
effects. For the purposes of Theorem 4.2, the well-preparedness assumption
(5.8) can therefore simply be replaced by Assumption 4.1.

###  5.1. Computational tools for spherical calculus

In this section, we briefly recall several computational tools that will be
used throughout this work to compute derivatives and integrals on the sphere.
First, we recall that the Laplace–Beltrami operator $\Delta_{n}$ on the sphere
$\mathbb{S}^{d-1}$ ($d\geq 2$) can be computed as follows: given a smooth
function $g:\mathbb{S}^{d-1}\to\mathbb{R}$, we can extend it to
$\mathbb{R}^{d}\setminus\\{0\\}$ by setting
$G(x):=g\big{(}\tfrac{x}{|x|}\big{)}$, and we then have

$\Delta_{n}g\,=\,(\Delta_{x}G)|_{\mathbb{S}^{d-1}}.$ (5.9)

In particular, we can compute in this way

$\displaystyle\Delta_{n}(n_{i})$ $\displaystyle=$ $\displaystyle(1-d)n_{i},$
(5.10) $\displaystyle\Delta_{n}(n_{i}n_{j})$ $\displaystyle=$ $\displaystyle
2\delta_{ij}-2dn_{i}n_{j},$ $\displaystyle\Delta_{n}(n_{i}n_{j}n_{k})$
$\displaystyle=$ $\displaystyle
2(\delta_{ij}n_{k}+\delta_{ik}n_{j}+\delta_{jk}n_{i})-3(d+1)n_{i}n_{j}n_{k},$
$\displaystyle\Delta_{n}(n_{i}n_{j}n_{k}n_{l})$ $\displaystyle=$
$\displaystyle
2(\delta_{ij}n_{k}n_{l}+\delta_{kj}n_{i}n_{l}+\delta_{ki}n_{j}n_{l}+\delta_{il}n_{j}n_{k}+\delta_{lk}n_{i}n_{j}+\delta_{jl}n_{i}n_{k})$
$\displaystyle-4(d+2)n_{i}n_{j}n_{k}n_{l},$

and so on for higher-order polynomials. These formulas can be used to
explicitly invert $\Delta_{n}$ on mean-zero polynomial expressions: for any
trace-free symmetric matrix $A\in\mathbb{R}^{d\times d}$, we find for
instance,

$\displaystyle\Delta_{n}^{-1}(n)$ $\displaystyle=$
$\displaystyle-\tfrac{1}{d-1}n,$ (5.11)
$\displaystyle\Delta_{n}^{-1}\big{(}n\otimes n:A\big{)}$ $\displaystyle=$
$\displaystyle-\tfrac{1}{2d}\,n\otimes n:A,$ (5.12)
$\displaystyle\Delta_{n}^{-1}\Big{(}n(n\otimes n:A)\Big{)}$ $\displaystyle=$
$\displaystyle-\tfrac{1}{3(d+1)}\Big{(}n(n\otimes
n:A)+\tfrac{4}{d-1}An\Big{)},$ (5.13)
$\displaystyle\Delta_{n}^{-1}\left((n\otimes
n:A)^{2}-\tfrac{2}{d(d+2)}\operatorname{tr}(A^{2})\right)$ $\displaystyle=$
$\displaystyle-\tfrac{1}{4(d+2)}\Big{(}(n\otimes n:A)^{2}+\tfrac{4}{d}n\otimes
n:A^{2}\hskip 28.45274pt$ (5.14) $\displaystyle\hskip
71.13188pt-\tfrac{2}{d}(\tfrac{1}{d+2}+\tfrac{2}{d})\,\operatorname{tr}(A^{2})\Big{)}.$

Henceforth, the pseudo-inverse $\Delta_{n}^{-1}$ is chosen to be defined as an
operator from mean-zero fields to mean-zero fields.

We also recall that the divergence of functions on $\mathbb{S}^{d-1}$ can be
computed similarly as the Laplace–Beltrami operator (5.9) by an extension
procedure: for any trace-free matrix $A$ and any smooth function
$g:\mathbb{S}^{d-1}\to\mathbb{R}$, we find for instance

$\displaystyle{\operatorname{div}}_{n}(\pi_{n}^{\bot}Ang)\,=\,\nabla_{n}g\cdot
An-d\,(n\otimes n):Ag.$ (5.15)

Finally, we further note that the above differential formulas (5.10) imply by
direct integration the following elementary integral identities for polynomial
expressions on the sphere,

$\displaystyle\int_{\mathbb{S}^{d-1}}n_{i}n_{j}\,\mathrm{d}n$ $\displaystyle=$
$\displaystyle\tfrac{\omega_{d}}{d}\delta_{ij},$ (5.16)
$\displaystyle\int_{\mathbb{S}^{d-1}}n_{i}n_{j}n_{k}n_{l}\,\mathrm{d}n$
$\displaystyle=$
$\displaystyle\tfrac{\omega_{d}}{d(d+2)}\big{(}\delta_{ij}\delta_{kl}+\delta_{kj}\delta_{il}+\delta_{ki}\delta_{jl}\big{)},$
$\displaystyle\int_{\mathbb{S}^{d-1}}n_{i}n_{j}n_{k}n_{l}n_{m}n_{p}\,\mathrm{d}n$
$\displaystyle=$
$\displaystyle\tfrac{\omega_{d}}{d(d+2)(d+4)}\big{(}\delta_{ij}\delta_{kl}\delta_{mp}+\ldots\big{)},$

and so on, where we recall the notation $\omega_{d}=|\mathbb{S}^{d-1}|$. These
identities imply in particular, for any trace-free symmetric matrices
$A,B\in\mathbb{R}^{d\times d}$,

$\displaystyle\int_{\mathbb{S}^{d-1}}\big{(}n\otimes
n-\tfrac{1}{d}{\operatorname{Id}}\big{)}\,\mathrm{d}n$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle\int_{\mathbb{S}^{d-1}}\big{(}n\otimes
n-\tfrac{1}{d}\operatorname{Id}\big{)}\,(n\otimes n:A)\,\mathrm{d}n$