1. Introduction
The problem of understanding the rheological properties of complex fluids has a rich history. Complex fluids include colloidal suspensions, solutions of macromolecules, emulsions and suspensions of cells or microorganisms. In such systems, the presence of supramolecular particles may not only modify the shear viscosity, but also generate non-Newtonian effects, depending on the shape of the particles and their interactions with the dispersion fluid. In this context, a central role is played by the characteristic relationship between the stress tensor and the macroscopic fields describing the coarse-grained velocity and composition in the complex fluid of interest. Such relationships can be established by empirical approaches or predicted using theoretical methods as pioneered by Einstein (Reference Einstein1906) for colloidal suspensions composed of micrometric solid particles. On spatial scales larger than the size of atoms and molecules, we can use the methods of hydrodynamics combined with suitable boundary conditions between the micrometric particles and the surrounding fluid since the interactions between atoms and molecules have subnanometric or nanometric ranges, which are significantly smaller than the size of the colloids. These theoretical methods show that the shape of the colloidal particles, as well as their specific properties have an important influence on the rheological properties (Landau & Lifshitz Reference Landau and Lifshitz1959; Batchelor Reference Batchelor1970; Bedeaux, Kapral & Mazur Reference Bedeaux, Kapral and Mazur1977; Hinch Reference Hinch1977; Brady, Khair & Swaroop Reference Brady, Khair and Swaroop2006).
In addition, suspensions may be passive or active. In the latter case, the colloidal particles composing the suspension move in the fluid by harnessing energy from their local environment. Like in other kinds of active matter, the propulsion of these particles results from the local conversion of the energy supplied by the surroundings into mechanical motion (Marchetti et al. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Aditi Simha2013; Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016). Instead, in passive suspensions, the non-equilibrium driving forces are induced by the gradients of the macrofields, which are defined on spatial scales larger than the distance separating the colloidal particles.
Recently, the dynamics and the rheology of passive non-uniform-slip Janus spheres (Ramachandran & Khair Reference Ramachandran and Khair2009; Premlata & Wei Reference Premlata and Wei2021, Reference Premlata and Wei2022) and Janus drops (Díaz-Maldonado & Córdova-Figueroa Reference Díaz-Maldonado and Córdova-Figueroa2015, Reference Díaz-Maldonado and Córdova-Figueroa2016) have been studied with these methods. (These particles or drops are named Janus because they are made of two parts side by side.) Similar studies have also been carried out for active suspensions of stick-slip squirmers (Yang et al. Reference Yang, Chen, Premlata and Wei2023) and rodlike microswimmers modelling unicellular microorganisms (Saintillan Reference Saintillan2010; Saintillan & Shelley Reference Saintillan and Shelley2013, Reference Saintillan and Shelley2015; Saintillan Reference Saintillan2018). In bacterial suspensions, the viscosity has even been observed to drop to zero in a kind of superfluid state (López et al. Reference López, Gachelin, Douarche, Auradou and Clément2015).
Remarkable effects also exist in active colloidal suspensions composed of synthetic solid particles (Paxton et al. Reference Paxton, Kistler, Olmeda, Sen, St. Angelo, Cao, Mallouk, Lammert and Crespi2004; Fourniez-Bidoz et al. Reference Fournier-Bidoz, Arsenault, Manners and Ozin2005; Wang et al. Reference Wang, Duan, Ahmed, Mallouk and Sen2013; Sánchez et al. Reference Sánchez, Soler and Katuri2014). In an important family of such systems, the solid particles are partly coated with some catalytic material, where the fuel coming from the surrounding fluid can react by heterogeneous catalysis and generate local concentration gradients, driving the particle by self-diffusiophoresis or self-electrophoresis (Anderson Reference Anderson1989; Anderson & Prieve Reference Anderson and Prieve1991; Würger Reference Würger2015; Ibrahim, Golestanian & Liverpool Reference Ibrahim, Golestanian and Liverpool2017; Moran & Posner Reference Moran and Posner2017; Oshanin, Popescu & Dietrich Reference Oshanin, Popescu and Dietrich2017; Kuron, Kreissl & Holm Reference Kuron, Kreissl and Holm2018).
Examples of such active Janus particles are silica spheres with a platinum cap and moving with speeds of the order of
$10\,\unicode{x03BC}$
m s−1 in an aqueous solution of hydrogen peroxide (Ke et al. Reference Ke, Ye, Carroll and Showalter2010; Valadares et al. Reference Valadares, Tao, Zacharia, Kitaev, Galembeck, Kapral and Ozin2010). In another example, hydrazine is used as fuel and iridium as catalyst, leading to higher speeds of approximately
$20\,\unicode{x03BC}$
m s−1 (Gao et al. Reference Gao, Pei, Dong and Wang2014).
In these examples, the propulsion mechanism is based on diffusiophoresis, in which a velocity slippage is generated between the fluid and the solid particle, as a consequence of local molecular concentration gradients along the fluid–solid interface and, thus, gradients in the interaction forces between the molecules of the fluid and the atoms composing the different parts of the particle surface (Golestanian et al. Reference Golestanian, Liverpool and Ajdari2005, Reference Golestanian, Liverpool and Ajdari2007; Kapral Reference Kapral2013; Colberg et al. Reference Colberg, Reigh, Robertson and Kapral2014). The so-generated interfacial velocity slippage induces a local flow of the fluid around the particle, which leads to its propulsion (Reigh et al. Reference Reigh, Huang, Schofield and Kapral2017; Campbell et al. Reference Campbell, Ebbens, Illien and Golestanian2019).
These effects can be described in terms of the constitutive properties of the aqueous solution and the fluid–solid interfaces, combined with the geometry of the Janus particles. The aqueous solution is an incompressible fluid characterised by its shear viscosity and the diffusion coefficients of the fuel and product molecules. The constitutive properties of the fluid–solid interfaces include the rate constants of the catalytic surface reactions, the diffusiophoretic constants and the possible Navier slip lengths. Since Janus particles have sizes ranging from hundreds of nanometres to a few micrometres, there is a large separation of length scales with respect to the diameters of the atoms and molecules composing the system, so that the fluid can be described in terms of the Navier–Stokes equation of hydrodynamics coupled by interfacial reactions and diffusiophoresis to the advection–diffusion equations for the molecular concentrations, while the active colloids can be described as rigid Janus particles. In this way, the macroscopic coarse-grained equations ruling the active suspension as a continuous medium can be deduced from the mesoscopic description based on the chemohydrodynamics of the fluid containing the colloidal particles. Such an approach allows us to predict the macroscopic properties of active matter and to justify the coarse-grained theories from the underlying structure and dynamics of its components.
This approach has already been developed for dilute active suspensions, where the particles move in a fluid that is globally at rest (Gaspard & Kapral Reference Gaspard and Kapral2020). In dilute suspensions, the particles are largely separated from each other, so that their mutual mechanical interactions can be neglected and they can be considered as isolated in the fluid. As a consequence, the macroscopic diffusion–reaction equations ruling the colloidal and molecular species can be deduced together with the complete expression for the entropy production rate density of such dilute active suspensions, as obtained by Gaspard & Kapral (Reference Gaspard and Kapral2020).
However, the extension of these results to active suspensions with a fluid flowing on large scales requires the knowledge of the stress tensor, including its active part (Hatwalne et al. Reference Hatwalne, Ramaswamy, Rao and Simha2004; Kruse et al. Reference Kruse, Joanny, Jülicher, Prost and Sekimoto2005; Ramaswamy Reference Ramaswamy2010; Marchetti et al. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Aditi Simha2013; Cates & Tjhung Reference Cates and Tjhung2018; Jülicher et al. Reference Jülicher, Grill and Salbreux2018). In principle, the stress tensor can be deduced from the underlying fluid mechanics of the suspension, using the methods developed since 1906 in the cases of passive and active systems (Einstein Reference Einstein1906; Landau & Lifshitz Reference Landau and Lifshitz1959; Batchelor Reference Batchelor1970; Bedeaux et al. Reference Bedeaux, Kapral and Mazur1977; Hinch Reference Hinch1977; Brady et al. Reference Brady, Khair and Swaroop2006; Ramachandran & Khair Reference Ramachandran and Khair2009; Lauga & Michelin Reference Lauga and Michelin2016; Saintillan Reference Saintillan2018; Premlata & Wei Reference Premlata and Wei2022). In this regard, an open issue is to understand how the stress tensor may depend both on the concentrations of the molecular species supplying chemical energy and on the gradients of these concentrations.
Schematic representation of an active suspension of Janus particles immersed in a shear flow.

The purpose of this paper is to carry out this programme for a dilute suspension of axisymmetric Janus particles immersed in an aqueous solution containing fuel and product species, which react on a catalytic hemisphere of the particles. A schematic representation of the system of interest is shown in figure 1. The study of this system is important not only because it has experimental realisations, but also because the mathematical deduction of the active stress tensor in dilute suspensions made of synthetic Janus particles can provide guiding lines to infer the active stress tensor in more complex forms of active matter such as active gels, suspensions of bacteria or other unicellular microorganisms, and biological tissues like muscles.
Our goal is thus to deduce the stress tensor of the active suspension of Janus particles under dilute conditions, following the approach developed by Gaspard & Kapral (Reference Gaspard and Kapral2020) to include all the terms that are linear in the thermodynamic forces, i.e. the affinity of chemical reactions and the gradients of the velocity and concentration fields (de Groot & Mazur Reference de Groot and Mazur1984). The novelty of this approach is that the description is based on the distribution function of the position and the orientation of the Janus particles, providing a clear understanding of the ways the couplings between the fluid flow, the diffusion of molecular species and their reactions may depend on the orientation of the particles. Afterwards, the macroscopic theories formulated in terms of the macrofields of polar or nematic ordering (Hatwalne et al. Reference Hatwalne, Ramaswamy, Rao and Simha2004; Kruse et al. Reference Kruse, Joanny, Jülicher, Prost and Sekimoto2005; Ramaswamy Reference Ramaswamy2010; Marchetti et al. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Aditi Simha2013; Cates & Tjhung Reference Cates and Tjhung2018; Jülicher et al. Reference Jülicher, Grill and Salbreux2018) can be deduced therefrom. The mesoscopic and macroscopic levels of description of active suspensions can thus be bridged together within the framework of this approach.
After this introduction, the plan of the paper is the following. The chemohydrodynamics of dilute active suspensions is presented in § 2. Therein, these suspensions are described from the scale of the micrometric colloidal particles up to the macroscale, where spatial averages are taken over the statistical distribution of the position and orientation of the particles. We consider hydrodynamically spherical Janus particles made of two hemispheres: one is catalytic for the chemical reaction powering the activity and the other is chemically inert. Each hemisphere and each chemical species is characterised by its own diffusiophoretic constant. The Navier slip length is supposed to be uniform on the whole sphere of the particle surface. We assume that there is no external force field exerted on the suspension, so that the particles are force free and torque free. Furthermore, the suspension is assumed to be dilute, at low Reynolds and Péclet numbers, and in the quasi-stationary approximation. Under these conditions, the equations for the velocity, the pressure and the concentration fields can be linearised and decoupled from each other, reducing to Stokes’ and Laplace’s equations, respectively. Therefore, the disturbances of the velocity and pressure fields can be obtained around isolated and individual particles moving in a background with velocity and concentration gradients and, subsequently, the effective stress tensor can be deduced using the aforementioned methods. The stationary fields around a single isolated Janus particle are obtained in § 3 by expanding the concentration and pressure fields into series of spherical harmonics and the velocity field into series of vectorial spherical harmonics. The details of the calculations for the one-particle profiles of molecular concentrations are given in Appendix A. After a summary of the definitions and relevant properties of the vectorial spherical harmonics in Appendix B, the calculations of the velocity field around a single isolated Janus particle are detailed in Appendix C. The stress tensor of the dilute active suspension is deduced in § 4 using the calculations of Appendix D based on the leading contributions of the velocity and pressure disturbances at large distances from the particle. For the hydrodynamically spherical particles we consider, the stress tensor is shown to have contributions involving not only the molecular concentrations but also the concentration gradients. The case where there is no chemical reaction is discussed in § 5. Section 6 is devoted to the case of Janus particles with uniform diffusiophoretic constants in the presence of reaction. The results obtained by Gaspard & Kapral (Reference Gaspard and Kapral2020) for the rectilinear and angular velocities are adapted in Appendix E to the case of the Janus particles that are discussed in § 6. Section 7 gives a summary and perspectives.
2. Chemohydrodynamics of dilute active suspensions
For the purpose of deriving the stress tensor for a dilute suspension of active Janus particles, the mathematical formulation is here set up based on chemohydrodynamics, which rules the motion of the suspension containing force-free and torque-free particles.
2.1. Description of the system
The suspension is composed of many spherical Janus particles of micrometric radius
$R$
immersed in an incompressible aqueous solution containing the solute molecular species
$\textrm {A}$
and
$\textrm {B}$
, which are mixed with a solvent molecular species
$\textrm {S}$
. The Janus particles have catalytic and non-catalytic hemispheres, which are denoted
$\textrm {c}$
and
$\textrm {n}$
, respectively. The reversible chemical reaction
happens on the catalytic hemisphere with rate constants per unit area
$\kappa _\pm ^{ {c}}$
,
$\textrm {C}$
denoting the colloidal particles as third solute species in addition to the molecular species
$\textrm {A}$
and
$\textrm {B}$
. The solvent species
$\textrm {S}$
is assumed to be chemically inert. The Janus particles are propelled by self-diffusiophoresis caused by the local concentration profiles of
$\textrm {A}$
and
$\textrm {B}$
generated around each particle by the reaction (2.1). There is no reaction in the bulk of the solution.
The suspension is described in terms of the velocity
$\boldsymbol{v}(\boldsymbol{r},t)$
of the fluid and the molecular concentrations
$c_{{A}}(\boldsymbol{r},t)$
and
$c_{{B}}(\boldsymbol{r},t)$
, which are fields defined in space
$\boldsymbol{r}=(x,y,z)$
outside the particles and at time
$t$
. In addition, the suspension contains
$N_{{C}}$
colloidal particles, which are labelled with the index
$c=1,2,\ldots ,N_{{C}}$
. Each Janus particle has a position
$\boldsymbol{R}_c$
, an orientation given by the unit vector
$\boldsymbol{u}_c$
with direction along the symmetry axis and pointing towards its catalytic hemisphere, a rectilinear velocity
$\boldsymbol{V}_c$
, and an angular velocity
$\boldsymbol{\varOmega }_c$
. The orientation of each particle is equivalently described by the polar and azimuthal angles
$(\theta _c,\varphi _c)$
, such that
$\boldsymbol{u}_c=\sin \theta _c \cos \varphi _c \, \boldsymbol{1}_x+\sin \theta _c \sin \varphi _c \, \boldsymbol{1}_y+\cos \theta _c \, \boldsymbol{1}_z$
, where
$\{\boldsymbol{1}_x, \boldsymbol{1}_y, \boldsymbol{1}_z\}$
is a Cartesian basis of orthogonal unit vectors of the tridimensional space. The velocity field inside each Janus particle is thus given by
The velocity and the concentration fields satisfy boundary conditions at the interface
$\Vert \boldsymbol{r}-\boldsymbol{R}_c(t)\Vert =R$
between the particles
$c=1,2,\ldots ,N_{{C}}$
and the surrounding fluid.
2.2. Equations for the velocity field of the fluid
Since the fluid is incompressible and the colloidal particles are rigid, the velocity field of the suspension satisfies the condition
Furthermore, the velocity field obeys the Navier–Stokes equation
where
$\rho$
is the mass density of the fluid and
${\boldsymbol{\sigma }}=(\sigma _{\textit{ij}})$
is the stress tensor. The latter is related to the pressure tensor
$\boldsymbol{\unicode{x1D64B}}=(P_{\textit{ij}})$
by
${\boldsymbol{\sigma }}=-{\boldsymbol{\unicode{x1D64B}}}$
. The components of the stress tensor are given by
\begin{align} \sigma _{\textit{ij}} = -p \, \delta _{\textit{ij}} +\eta _0 \left (\boldsymbol{\nabla} _i v_{\!j}+\boldsymbol{\nabla} _{\!j} v_i\right ) + \sum _{c=1}^{N_{{C}}} \sigma _{\textit{ij}}^{(c)} , \end{align}
where
$p$
is the pressure,
$\eta _0$
the shear viscosity of the fluid and
$\sigma _{\textit{ij}}^{(c)}$
the contribution of the
$c{\mbox{th}}$
colloidal particle to the stress tensor. The latter can be expressed as a distribution that is singular at the interface between the particle and the fluid and that is determined by the boundary conditions on the velocity field (Mazur & Bedeaux Reference Mazur and Bedeaux1974; Albano, Bedeaux & Mazur Reference Albano, Bedeaux and Mazur1975; Gaspard & Kapral Reference Gaspard and Kapral2018, Reference Gaspard and Kapral2019). These boundary conditions include the effects of sliding friction due to possible partial slip between the fluid and the particle, and by diffusiophoresis coupling the fluid velocity with the concentrations of the molecular species. The surface
$\Vert \boldsymbol{r}-\boldsymbol{R}_c(t)\Vert =R$
of the
$c{\mbox{th}}$
particle is denoted
$\varSigma _c(t)$
. The unit vector normal to this surface and oriented towards the fluid is denoted
$\boldsymbol{n}_c=R^{-1}(\boldsymbol{r}-\boldsymbol{R}_c)$
for
$\boldsymbol{r}\in \varSigma _c(t)$
. In the direction normal to the interface, the boundary condition reads
In the two directions tangential to the interface, they are given by
\begin{align} \left ({\boldsymbol{\unicode{x1D644}}}-\boldsymbol{n}_c\boldsymbol{n}_c\right ) \boldsymbol{\cdot }\left [\boldsymbol{v} -\boldsymbol{v}_{\textit{solid}}^{(c)} - b\big (\boldsymbol{\nabla }\boldsymbol{v}+\boldsymbol{\nabla }\boldsymbol{v}^{{T}}\big )\boldsymbol{\cdot }\boldsymbol{n}_c +\sum _\varkappa b_\varkappa \boldsymbol{\nabla }c_\varkappa \right ]_{\varSigma _c} = 0 , \end{align}
where
${\boldsymbol{\unicode{x1D644}}}=(\delta _{\textit{ij}})$
is the tridimensional identity matrix,
$b$
is the Navier slip length characterising the possible partial slip between the particle and the fluid, and
$b_\varkappa$
is the diffusiophoretic constant of the molecular species
$\varkappa \in \{\textrm {A},\textrm {B}\}$
. Here,
$(\boldsymbol{\nabla }\boldsymbol{v})^{{T}}$
denotes the transpose of the matrix
$\boldsymbol{\nabla }\boldsymbol{v}$
.
The Navier slip length
$b$
is assumed to be uniform on the whole interface between the spherical particle and the fluid. The diffusiophoretic constants are supposed to take different values on each hemisphere
where
$H^{h}(\theta _c)$
denotes the Heaviside function that is equal to one on hemisphere
$h$
and zero otherwise. The catalytic hemisphere is taken as
$0\leqslant \theta _c\leqslant ({\pi }/{2})$
and the non-catalytic hemisphere as
$({\pi }/{2})\lt \theta _c\leqslant \pi$
in the frame where the polar axis is oriented in the direction of the unit vector
$\boldsymbol{u}_c$
.
As a remark, we should note that the boundary conditions (2.7) result from the thin-interaction-layer approximation (Anderson Reference Anderson1989; Michelin & Lauga Reference Michelin and Lauga2014), holding if the diffusiophoretic lengths
$\lambda _\varkappa \equiv \sqrt {b_\varkappa \eta _0/k_{{B}}T}$
(
$k_{{B}}$
denotes the Boltzmann constant) as well as the Navier slip length
$b$
are much smaller than the particle radius
$R$
. These interfacial characteristic lengths are in the range of
$0.1$
–
$50$
nm (Joly, Ybert & Bocquet Reference Joly, Ybert and Bocquet2006; Bouzigues, Tabeling & Bocquet Reference Bouzigues, Tabeling and Bocquet2008; Li et al. Reference Li, D’eramo, Lee, Monti, Yonger, Tabeling, Chollet, Bresson and Tran2015; Ibrahim et al. Reference Ibrahim, Golestanian and Liverpool2017; Collis et al. Reference Collis, Olcum, Chakraborty, Manalis and Sader2021), so that the thin-interaction-layer approximation is valid if the particles are micrometric.
2.3. Equations for the molecular concentrations
The concentrations of the molecular species
$\varkappa \in \{\textrm {A},\textrm {B}\}$
satisfy the following advection–diffusion–reaction equations
\begin{align} \left (\partial _t + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) c_\varkappa = D_\varkappa {{\nabla} }^2 c_\varkappa + \sum _{c=1}^{N_{{C}}} \nu _\varkappa w^{(c)} , \end{align}
where
$D_\varkappa$
is the associated molecular diffusivity,
$\nu _\varkappa$
is the stoichiometric coefficient of species
$\varkappa$
in the reaction
$\textrm {A}\to \textrm {B}$
of (2.1) (
$\nu _{{A}}=-1$
and
$\nu _{{B}}=+1$
) and
$w^{(c)}$
is the contribution of the
$c{\mbox{th}}$
particle to the reaction rate density. (If several reactions occurred,
$\nu _\varkappa w^{(c)}$
should be replaced by
$\sum _r \nu _{\varkappa r} w_r^{(c)}$
, where
$\nu _{\varkappa r}$
is the stoichiometric coefficient of species
$\varkappa$
in reaction
$r$
and
$w_r^{(c)}$
is the corresponding rate density.) This contribution is also singular at the interface between the Janus particle and the fluid. It is fixed by the boundary conditions
where
are the rate constants per unit area of the forward and reverse surface reactions in (2.1).
The overall reaction rate of the
$c{\mbox{th}}$
Janus particle is defined as the number of reactive events
$\textrm {A}\to \textrm {B}$
per unit time taking place on the catalytic hemisphere of the particle, A being considered as the fuel and B as the product. This overall reaction rate is given by the incoming flux of fuel molecules A (or the outgoing flux of product molecules B) at the surface of the Janus particle,
as expressed in terms of the surface area element
$\textrm {d}\boldsymbol{\varSigma }=\boldsymbol{n}_c \textrm { d}\varSigma$
and the rate constants in (2.11). As a consequence, the total number of solute molecules,
$N_0\equiv N_{{A}}+N_{{B}}$
with
$N_\varkappa \equiv \int _{\mathcal V} c_\varkappa \, {\rm d}^3r$
, is a conserved quantity for the reaction (2.1), i.e.
$\textrm {d}N_0/\textrm {d}t=0$
.
We note that the suspension is active if the molecular concentrations are not in their Guldberg–Waage equilibrium ratio, i.e. if
$\kappa _+^{{c}}\, c_{{A}} \ne \kappa _-^{{c}} \, c_{{B}}$
. Otherwise, the suspension is said to be passive. Even if the suspension is at chemical equilibrium, i.e. if
$\kappa _+^{{c}}\, c_{{A}} = \kappa _-^{{c}} \, c_{{B}}$
, the system may evolve out of equilibrium under the effects of non-zero gradients of velocity or molecular concentrations. The suspension is thus active because the chemical reaction is out of equilibrium.
2.4. Force and torque exerted by the fluid on each particle
Once the velocity field of the fluid is known, the force and the torque exerted by the fluid on each colloidal particle are respectively given by
where
$\boldsymbol{\sigma }^{({fl})}=-p\, \boldsymbol{\unicode{x1D644}}+\eta _0(\boldsymbol{\nabla }\boldsymbol{v}+\boldsymbol{\nabla }\boldsymbol{v}^{{T}})$
is the stress tensor of the fluid around the particle and
$c=1,2,\ldots ,N_{{C}}$
.
In the absence of an external force field, the particles are force free and torque free, which determines their rectilinear and angular velocities with respect to the fluid flow. In this way, the equations of motion can be obtained for the position
$\boldsymbol{R}_c(t)$
and the orientation
$\boldsymbol{u}_c(t)$
of each particle
$c=1,2,\ldots ,N_{{C}}$
.
2.5. Distribution function of the Janus particles
The distribution of the position and orientation of the Janus particles is defined by

where
$\delta ^3[\boldsymbol{r}-\boldsymbol{R}_c(t)]$
is the tridimensional Dirac distribution centred at the position
$\boldsymbol{R}_c$
of the
$c{\mbox{th}}$
particle and
$\delta ^2[\boldsymbol{u}-\boldsymbol{u}_c(t)]$
stands for
$\delta (\theta -\theta _c)\delta (\varphi -\varphi _c)/\sin \theta _c$
. The distribution (2.15) is normalised to the invariant total number of particles according to
with
${\rm d}^2u = \textrm {d}\varOmega =\sin \theta\, \textrm {d}\theta \,\textrm {d}\varphi$
.
The suspension may be described as a continuous medium over length scales larger than the size of coarse-graining cells
${\mathcal V}_{\boldsymbol{r}}$
of volume
$V_{\boldsymbol{r}}$
and centred around the positions
$\boldsymbol{r}$
, if the statistical properties of the suspension may be assumed to remain quasi-uniform in these cells. The cells
${\mathcal V}_{\boldsymbol{r}}$
should contain a sufficiently large number of particles to justify a statistical formulation. Accordingly, we may introduce the spatial averages of the quantities of interest over the cells
${\mathcal V}_{\boldsymbol{r}}$
as
In particular, the mean distribution function of the particles may thus be evaluated as

such that
for any quantity
$(\boldsymbol{\cdot })$
depending on the coordinates
$\{\boldsymbol{R}_c,\boldsymbol{u}_c\}_{c=1}^{N_{ { C}}}$
of the Janus particles. We note that we may also consider a statistical ensemble of configurations for the Janus particles and the corresponding statistical average
$\langle \boldsymbol{\cdot }\rangle$
, which should be equal to the spatial average,
, if ergodicity is held (Batchelor Reference Batchelor1970).
The density of colloidal particles is given by
In addition, we may introduce the moments of the distribution of the orientation of the Janus particles as
\begin{align} r_{\textit{ijk}}(\boldsymbol{r},t) &\equiv \int {\rm d}^2u \, f_{{C}}(\boldsymbol{r},\boldsymbol{u},t) \, R_{\textit{ijk}}(\boldsymbol{u}) \nonumber \\&\qquad \mbox{with}\quad R_{\textit{ijk}}(\boldsymbol{u}) \equiv u_i \, u_{\!j} \, u_k - \frac {1}{5} \, u_i \, \delta _{jk} - \frac {1}{5} \, u_{\!j} \, \delta _{ik} - \frac {1}{5} \, u_k \, \delta _{\textit{ij}} , \end{align}
such that the distribution function can be decomposed as
\begin{align} f_{{C}}(\boldsymbol{r},\boldsymbol{u},t) &= \frac {1}{4\pi } \, n_{{C}}(\boldsymbol{r},t) + \frac {3}{4\pi } \, p_i(\boldsymbol{r},t) \, u_i + \frac {15}{8\pi } \, q_{\textit{ij}}(\boldsymbol{r},t) \, Q_{\textit{ij}}(\boldsymbol{u})\nonumber\\&\quad + \frac {35}{8\pi } \, r_{\textit{ijk}}(\boldsymbol{r},t) \, R_{\textit{ijk}}(\boldsymbol{u}) + {\cdots} \end{align}
with the convention of summation over repeated indices (Turzi Reference Turzi2011; Saintillan & Shelley Reference Saintillan and Shelley2013). We note that the tensors
$q_{\textit{ij}}$
and
$Q_{\textit{ij}}$
are symmetric and traceless, while the tensors
$r_{\textit{ijk}}$
and
$R_{\textit{ijk}}$
are totally symmetric in the indices and traceless in two indices, i.e. satisfying
$r_{iik}=0$
and
$R_{iik}=0$
, as well as similar identities obtained by the permutation of indices.
As an alternative to the moments (2.20), (2.21), (2.22),
$\ldots$
, we may introduce the local mean values
We note that the mean distribution function (2.17) obeys a local conservation equation generalised to the configuration space of the position and orientation of the particles:
where
$\boldsymbol{J}_{ {Ct}}$
and
$\boldsymbol{J}_{{Cr}}$
denote the translational and rotational components of the colloidal current density, and
$\textrm {div}_{{r}}$
is the divergence in the orientation subspace (Brenner Reference Brenner1965). Equation (2.25) is often referred to as the Smoluchowski equation if the colloidal current density includes the diffusion terms (Saintillan & Shelley Reference Saintillan and Shelley2015).
2.6. Dilute suspensions
The suspension is dilute if the mean distance
${\mathcal R}_{{C}}$
between nearest neighbouring particles is much larger than their radius
$R$
:
In this case, the mechanical interactions between the colloidal particles can be neglected and the suspension behaves as a collection of largely separated and quasi-isolated particles immersed in the coarse-grained velocity and concentration fields.
Since the cells
${\mathcal V}_{\boldsymbol{r}}$
used for coarse graining contain many colloidal particles, they can be decomposed as
${\mathcal V}_{\boldsymbol{r}}=\bigcup _{c:\,\boldsymbol{R}_c\in {\mathcal V}_{\boldsymbol{r}}}{\mathcal V}_{\boldsymbol{r}}^{(c)}$
into smaller domains
${\mathcal V}_{\boldsymbol{r}}^{(c)}$
, each containing a single particle, for instance, using a Voronoi decomposition. The size of these domains is of the order of the mean distance between nearest neighbours,
${\mathcal R}_{{C}}\sim n_{{C}}^{-1/3}$
, satisfying the condition (2.26) in a dilute suspension. As a consequence, the rim of each domain
${\mathcal V}_{\boldsymbol{r}}^{(c)}$
is at large distances from the Janus particle located at its centre. Inside each domain
${\mathcal V}_{\boldsymbol{r}}^{(c)}$
, the profiles of the velocity, pressure and concentration fields are expected to take the following forms:
as expressed in a frame centred on the position of each particle, i.e. such that
$\boldsymbol{R}_c=0$
. In (2.27),
${\boldsymbol{\unicode{x1D642}}}=(G_{\textit{ij}})$
denotes a possible tensor of velocity gradients holding on average across the domain
${\mathcal V}_{\boldsymbol{r}}^{(c)}$
. Since the fluid is incompressible, this tensor is traceless,
$\textrm {tr}\,{\boldsymbol{\unicode{x1D642}}}=G_{ii}=0$
. Moreover,
$\boldsymbol{v}^{\circ }$
is a constant value for the fluid velocity and
$\Delta \boldsymbol{v}$
is the disturbance velocity due to the presence of the Janus particle in the centre of the domain
${\mathcal V}_{\boldsymbol{r}}^{(c)}$
. This disturbance velocity vanishes at large distances from the particle. Similarly, in (2.28) and (2.29),
$p^{\circ }$
and
$c_\varkappa ^{\circ }$
denote constant values, and
$\Delta p$
and
$\Delta c_\varkappa$
the disturbance pressure and concentrations, vanishing away from the particle. In (2.29),
$\boldsymbol{g}_\varkappa$
denotes a possible gradient in the profile of the molecular concentration
$c_\varkappa$
, holding on average across the domain
${\mathcal V}_{\boldsymbol{r}}^{(c)}$
.
2.7. Linearisation of the equations of motion and quasi-stationarity
The equations of chemohydrodynamics are linearised around each particle in the dilute-suspension limit (2.26) to obtain the disturbances of the velocity, pressure and concentration fields, i.e.
$\Delta \boldsymbol{v}$
,
$\Delta p$
and
$\Delta c_\varkappa$
in (2.27), (2.28) and (2.29), respectively.
For this purpose, we suppose that the flow is laminar with low Reynolds numbers,
and low Péclet numbers,
with respect to the speed of the particles
$V_c\equiv \Vert \boldsymbol{V}_c\Vert$
and the characteristic speed associated with the velocity gradients in shear flows
$V_{\boldsymbol{\unicode{x1D642}}}=R\Vert {\boldsymbol{\unicode{x1D642}}}\Vert$
. (The norm of the tensor of velocity gradients can be defined as
$\Vert {\boldsymbol{\unicode{x1D642}}}\Vert \equiv \sup _{\Vert \boldsymbol{n}\Vert =1}\Vert \boldsymbol{\unicode{x1D642}}\boldsymbol{\cdot }\boldsymbol{n}\Vert$
.) The conditions with respect to the particle speed
$V_c$
are satisfied in aqueous suspensions, where the active particles are micrometric (
$R\sim 1 \,\unicode{x03BC}$
m) and have speeds of the order of
$V_c \sim 10\, \unicode{x03BC}$
m s−1, so that
$Re\sim 10^{-5}$
for
$\eta _0/\rho \simeq 10^{-6}\,$
m
$^2$
s−1, and
$Pe_\varkappa \sim 10^{-2}$
for
$D_\varkappa \sim 10^{-9}\,$
m
$^2$
s−1. The conditions with respect to the shear-flow characteristic speed
$V_{\boldsymbol{\unicode{x1D642}}}$
are fulfilled if the velocity gradients are small enough, satisfying
$\Vert {\boldsymbol{\unicode{x1D642}}}\Vert \ll \eta _0/(\rho R^2)\sim 10^6\,$
s
$^{-1}$
and
$\Vert {\boldsymbol{\unicode{x1D642}}}\Vert \ll D_\varkappa /R^2\sim 10^3\,$
s
$^{-1}$
.
Under these conditions, the Navier–Stokes equation (2.4) reduces to the linear Stokes equation, and the advection–diffusion–reaction equation (2.9) can also be linearised by neglecting
$\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }(\boldsymbol{\cdot })$
in front of
$\partial _t(\boldsymbol{\cdot })$
. With these linearisations, the equations for the molecular concentrations are decoupled from the equations for the velocity field.
Furthermore, the velocity and concentration fields are assumed to vary over time scales that are longer than those associated with viscosity and molecular diffusion,
$t_{{visc}}\sim \rho R^2/\eta _0$
and
$t_{{ diff}}\sim R^2/D_\varkappa$
, respectively. As a consequence, we may consider stationary solutions for the linearised (2.4) and (2.9) to determine the profiles of the fields around each Janus particle. Therefore, these equations reduce to
which should be solved with the boundary conditions (2.6), (2.7) and (2.10) for a single colloidal particle located at
$\boldsymbol{r}=0$
, and
at large distances from the particle, i.e. for
$\Vert \boldsymbol{r}\Vert \gg R$
.
These approximations are justified to the extent that the stress tensor is limited to the contributions that are linear in the gradients of the velocity and concentration fields. If the conditions (2.30) and (2.31) were not satisfied for either the particle speed
$V_c$
or the shear-flow characteristic speed
$V_{\boldsymbol{\unicode{x1D642}}}$
, corrections of higher order in the gradients would be required beyond the linear contributions that are calculated here.
2.8. Effective stress tensor of the suspension
We suppose that the particles are force free, i.e. the forces (2.13) are equal to zero,
$\boldsymbol{F}_c=0$
. Under these circumstances, the disturbance velocity and pressure induced by each particle behave as
for
$r=\Vert \boldsymbol{r}\Vert \to \infty$
. In such dilute suspensions, the effective stress tensor can be calculated as follows using the spatial average (2.16):
\begin{align} \overline {\sigma }_{\textit{ij}} &= -\overline {p}\, \delta _{\textit{ij}} + \eta _0\, \overline {\boldsymbol{\nabla} _i v_{\!j}+\boldsymbol{\nabla} _{\!j} v_i} \nonumber\\&\quad + \frac {1}{V_{\boldsymbol{r}}} \sum _{c:\,\boldsymbol{R}_c\in {\mathcal V}_{\boldsymbol{r}}} \int _{{\mathcal V}_{\boldsymbol{r}}^{(c)}} \left [ \Delta \sigma _{ik}^{( {fl})} \, r_{\!j} \, \textrm {d}\varSigma _k - \eta _0 \left ( \Delta v_i \, \textrm {d}\varSigma _{\!j} + \Delta v_{\!j} \, \textrm {d}\varSigma _i \right )\right ]\! , \end{align}
where
$\Delta \sigma _{ik}^{( {fl})} = -\Delta p \, \delta _{ik} +\eta _0 ( \boldsymbol{\nabla} _i \Delta v_k + \boldsymbol{\nabla} _k \Delta v_i )$
is the disturbance in the stress tensor of the fluid due to the presence of a particle at the centre of the cell
${\mathcal V}_{\boldsymbol{r}}^{(c)}$
and evaluated at large distances from this centre, after having solved the one-particle problem (Landau & Lifshitz Reference Landau and Lifshitz1959; Batchelor Reference Batchelor1970). In (2.38), the spatial averages of the velocity gradients correspond to the velocity gradients in the domains
${\mathcal V}_{\boldsymbol{r}}^{(c)}$
as introduced in (2.27) according to
$\overline {\boldsymbol{\nabla} _i v_{\!j}}=G_{\textit{ji}}$
and
$\overline {\boldsymbol{\nabla} _{\!j} v_i}=G_{\textit{ij}}$
.
Furthermore, the particles are assumed to be torque free, i.e. the torques in (2.14) are equal to zero,
$\boldsymbol{T}_c=0$
. Consequently, the stress tensor of the suspension remains symmetric,
$\sigma _{\textit{ij}}=\sigma _{\textit{ji}}$
, as shown later.
The calculation of the stress tensor will be carried out for the dilute active suspension in the following sections.
3. Stationary fields around a single Janus particle
3.1. Molecular concentrations
The first step of the problem is to obtain the concentration profiles, because they are decoupled from the velocity field as a consequence of the low-Péclet-number limits (2.31) and the facts that the resulting (2.33) for the concentrations and the boundary conditions (2.10) and (2.35) do not depend on the velocity field. The calculation is presented in Appendix A by assuming that the Janus particle is oriented in the direction
$\boldsymbol{u}=\boldsymbol{1}_z$
. The quantity entering the boundary conditions (2.7) on the tangential components of the velocity field is thus given by
\begin{align} \left ( \sum _\varkappa b_\varkappa \boldsymbol{\nabla } c_\varkappa \right )_R = \frac {3}{2} \left ( \boldsymbol{\unicode{x1D644}} - \boldsymbol{n}\boldsymbol{n}\right ) \boldsymbol{\cdot }\left ( b_{{A}} \boldsymbol{g}_{{ A}} + b_{{B}} \boldsymbol{g}_{{B}}\right ) + R \left (\frac {b_{{B}}}{D_{{B}}} - \frac {b_{{A}}}{D_{{A}}}\right ) (\boldsymbol{\nabla } Z)_R , \end{align}
where the function
$Z$
obeys the equation
${\nabla} ^2Z=0$
with the boundary conditions
and
$(\boldsymbol{\nabla } Z)_\infty =0$
, as expressed in terms of the dimensionless Damköhler number
the reaction rate
and its gradient
The analytical expression of the function
$Z$
is given in Appendix A.
In the right-hand side of (3.1), the first contribution comes from simple diffusiophoresis in the gradients of concentrations and the second is due to self-diffusiophoresis generated by the reaction. In the absence of reaction (i.e. if
$\kappa _\pm ^{{c}}=0$
) or if the reaction is at equilibrium (i.e. if
$\kappa _+^{{c}} \, c_{{A}}=\kappa _-^{{c}} \, c_{{B}}$
), the second contribution is equal to zero.
The system evolves in the reaction-limited regime if
${\textit {Da}}\ll 1$
, i.e. if
$R\kappa _+^{{c}} \ll D_{{ A}}$
and
$R\kappa _-^{{c}} \ll D_{{B}}$
; or in the diffusion-controlled regime if
${\textit {Da}}\gg 1$
, i.e. if
$R\kappa _+^{{c}} \gg D_{{A}}$
or
$R\kappa _-^{{c}} \gg D_{{B}}$
.
3.2. Pressure and velocity fields
Taking the divergence of Stokes’ equation (2.32) and using the incompressibility condition (2.3), we obtain the Laplace equation
${\nabla} ^2p=0$
for the pressure. Around a single Janus particle located at
$\boldsymbol{r}=0$
, the pressure is expanded into spherical harmonics as
where
$l\geqslant \vert m \vert \geqslant 0$
and
This linear second-order homogeneous ordinary differential equation is solved by considering the following linear superposition:
with some coefficients
$a_{lm}$
and
$b_{lm}$
to be determined.
The vector field of the fluid velocity around the Janus particle can be expanded as
into the vectorial spherical harmonics that are defined in Appendix B (Barrera, Estévez & Giraldo Reference Barrera, Estévez and Giraldo1985).
In Appendix C, Stokes’ equation (2.32) and the incompressibility condition (2.3) are solved with the boundary conditions (2.6), (2.7) and (2.34) to obtain the radial functions
$f_{lm}(r)$
,
$g_{lm}(r)$
and
$h_{lm}(r)$
of the velocity field around a single Janus particle. These radial functions are expressed as powers of the radial distance
$r$
with coefficients given in terms of the velocity
$v_i^{\circ }$
and velocity gradients
$G_{\textit{ij}}$
, the particle rectilinear velocity
$V_i$
and angular velocity
$\varOmega _i$
, and the molecular concentrations
$c_\varkappa ^{\circ }$
and concentration gradients
$\boldsymbol{g}_\varkappa$
involved because of diffusiophoresis and the chemical reaction. The analytical expressions of the radial functions are calculated in Appendix C.
Finally, we find that the disturbance velocity and pressure can be expressed as
where the matrix of this stresslet has three parts:
given by (C92), (C101) and (C121), respectively. The first part is proportional to the traceless tensor of velocity gradients, the second one is due to simple diffusiophoresis and the third one to self-diffusiophoresis generated by the reaction. These three parts contribute to the stress tensor of the active suspension according to (2.38), as shown in the following section.
4. Stress tensor of the dilute active suspension
4.1. General expression
According to (2.38) and (2.18), the stress tensor of the dilute active suspension can be expressed as
with
in terms of the disturbance velocity and pressure given by (3.10). The calculation of (4.2) for the part of the stress tensor due to the colloidal particles is carried out in Appendix D. This calculation shows that
in terms of the stresslet matrix (3.11). The three parts of this matrix lead to three corresponding contributions to the stress tensor,
which are given in the following subsections using the results of the calculations carried out in Appendix C.
4.2. Contribution to the shear viscosity
Inserting the stresslet matrix given in (C92) into (4.3) and using the property that
$\textrm {tr}\, \boldsymbol{\unicode{x1D63C}}^{({v})}=A_{kk}^{({v})}=0$
show that
\begin{align} \Delta \sigma _{\textit{ij}}^{({v})} = \frac {5}{2} \, \frac {\displaystyle 1 + \frac {2b}{R}}{\displaystyle 1 + \frac {5b}{R}} \, \phi \, \eta _0 \left ( G_{\textit{ij}}+G_{\textit{ji}}\right )\! , \qquad \mbox{where}\quad \phi \equiv \frac {4\pi }{3}\, R^3 \, n_{{C}} \end{align}
is the volume fraction occupied by the Janus particles in the suspension. Since
$G_{\textit{ij}}=\overline {\boldsymbol{\nabla} _{\!j} v_i}$
, this contribution to the stress tensor (4.1) implies that the effective shear viscosity of the suspension is given by
\begin{align} \eta = \eta _0 \left (1 + \frac {5}{2} \, \frac {\displaystyle 1 + \frac {2b}{R}}{\displaystyle 1 + \frac {5b}{R}} \, \phi \right ) \end{align}
at leading order in the colloidal fraction, which is a well-known result since Einstein’s pioneering work (Einstein Reference Einstein1906; Happel & Brenner Reference Happel and Brenner1983; Luo & Pozrikidis Reference Luo and Pozrikidis2008; Premlata & Wei Reference Premlata and Wei2022). In the limit
$b\ll R$
, Einstein’s formula
$\eta =\eta _0(1+5\phi /2)$
is recovered. Corrections of order
$\phi ^2$
or higher could be calculated, in particular, with the methods of Bedeaux et al. (Reference Bedeaux, Kapral and Mazur1977) and Happel & Brenner (Reference Happel and Brenner1983).
We note that, as expected for hydrodynamically spherical particles with a uniform Navier slip length, the non-equilibrium driving forces given by the velocity gradients
$G_{\textit{ij}}$
generate only the contribution (4.5). Moreover, for the colloidal particles we consider, there are no non-Newtonian modifications like
$u_iu_kG_{jk}$
,
$u_{\!j}u_kG_{ik}$
or
$Q_{\textit{ij}}u_ku_lG_{kl}$
(Ramachandran & Khair Reference Ramachandran and Khair2009; Saintillan & Shelley Reference Saintillan and Shelley2013; Saintillan Reference Saintillan2018; Premlata & Wei Reference Premlata and Wei2022).
4.3. Contribution from simple diffusiophoresis
Here, we present the contribution from simple diffusiophoresis given by the stresslet matrix (C101). Inserting it into (4.3) leads to
\begin{align} \Delta \sigma _{\textit{ij}}^{({d})} = \frac {15\pi }{8} \, \frac {R^2 \eta _0}{\displaystyle 1 + \frac {5b}{R}} \int {\rm d}^2u \, f_{{C}} \sum _\varkappa \Delta b_\varkappa \left [ g_{\varkappa i} \, u_{\!j} + g_{\varkappa j} \, u_i+ \boldsymbol{g}_\varkappa \boldsymbol{\cdot }\boldsymbol{u} \left ( u_i \, u_{\!j} - \delta _{\textit{ij}} \right )\right ]\! , \end{align}
where
$\boldsymbol{g}_\varkappa =\overline {\boldsymbol{\nabla } c_\varkappa }$
are the molecular concentration gradients defined on the scale of the coarse-grained description and
$\Delta b_\varkappa = b_\varkappa ^{{c}} - b_\varkappa ^{{n}}$
are the differences of diffusiophoresis constants between the catalytic and non-catalytic hemispheres of the Janus particles. Therefore, this contribution is not equal to zero if there is heterogeneity in the diffusiophoretic effect at the surface of the particles. We note that this contribution to the stress tensor is symmetric,
$\Delta \sigma _{\textit{ij}}^{({ d})}=\Delta \sigma _{\textit{ji}}^{({d})}$
, and traceless,
$\textrm {tr}\, \Delta \boldsymbol{\sigma }^{({ d})}=\Delta \sigma _{ii}^{({d})}=0$
, as expected for torque-free rigid Janus particles.
Using the decomposition (2.23) of the distribution function of the Janus particles into the moments of their orientation (2.20)–(2.22), we find that
\begin{align} \Delta \sigma _{\textit{ij}}^{({d})} = \frac {9\pi }{4} \, \frac {R^2 \eta _0}{\displaystyle 1 + \frac {5b}{R}} \sum _\varkappa \Delta b_\varkappa \left ( g_{\varkappa i} \, p_{\!j} + g_{\varkappa j} \, p_i - \frac {2}{3} \, \boldsymbol{g}_\varkappa \boldsymbol{\cdot }\boldsymbol{p} \, \delta _{\textit{ij}} + \frac {5}{6} \, g_{\varkappa k} \, r_{\textit{ijk}} \right )\! , \end{align}
which is expressed in terms of the concentration gradients
$g_{\varkappa k}=\overline {\boldsymbol{\nabla} _k c_\varkappa }$
of the molecular species
$\varkappa \in \{\textrm {A},\textrm {B}\}$
, the vector field (2.20) and the rank-three tensor (2.22) of the particle orientation. Here, the non-equilibrium driving forces are the concentration gradients of the solute molecular species. This coupling between the stress tensor and the concentration gradients is allowed by Curie’s symmetry principle because the colloidal Janus particles have an orientation characterised by the polar-order vector (2.20) and the octupolar-order tensor (2.22). The contribution (4.8) would vanish in the absence of orientational order, i.e. if
$p_i=0$
and
$r_{\textit{ijk}}=0$
. We also note that the contribution (4.8) to the stress tensor can induce non-Newtonian effects associated with non-vanishing normal stress differences
$N_1\equiv \sigma _{11}-\sigma _{22}\ne 0$
or
$N_2\equiv \sigma _{22}-\sigma _{33}\ne 0$
, here depending on the concentration gradients.
4.4. Contribution from reaction
Finally, the contribution from the reaction can be obtained by inserting the stresslet matrix (C121) into the formula (4.3), giving
\begin{align} \Delta \sigma _{\textit{ij}}^{({r})} &= \frac {4\pi R^2 \eta _0}{\displaystyle 1 + \frac {5b}{R}} \int \textrm {d}^2u \, f_{{C}} \,\Bigg [\left (\alpha _\varLambda \, \varrho + \frac {3R}{2} \, \beta _\varLambda \, \boldsymbol{g}_{\varrho }\boldsymbol{\cdot }\boldsymbol{u}\right ) \left ( 3 \, u_i \, u_{\!j} -\delta _{\textit{ij}} \right ) \nonumber \\ &\quad + \frac {3R}{2} \, \gamma _\varLambda \left ( g_{\varrho i} \, u_{\!j}+ g_{\varrho j}\, u_i - 2 \, \boldsymbol{g}_\varrho \boldsymbol{\cdot }\boldsymbol{u} \, u_i \, u_{\!j} \right )\Bigg ] \end{align}
in terms of the reaction rate (3.4) with the average molecular concentrations
$c_\varkappa ^{\circ }=\overline {c_\varkappa }$
, its gradient (3.5) with
$\boldsymbol{g}_\varkappa =\overline {\boldsymbol{\nabla } c_\varkappa }$
and the coefficients (C116)–(C118). This contribution due to the self-diffusiophoresis generated by the reaction is also symmetric,
$\Delta \sigma _{\textit{ij}}^{({r})}=\Delta \sigma _{\textit{ji}}^{({r})}$
, and traceless,
$\textrm {tr}\, \Delta \boldsymbol{\sigma }^{({r})}=\Delta \sigma _{ii}^{({r})}=0$
.
Using the decomposition (2.23), this contribution can be expressed as follows:
\begin{align} \Delta \sigma _{\textit{ij}}^{({r})} &= \frac {12 \pi R^2 \eta _0}{\displaystyle 1 + \frac {5b}{R}} \, \Bigg [ \alpha _\varLambda \, \varrho \, q_{\textit{ij}} + \frac {3}{10} \left (\beta _\varLambda + \gamma _\varLambda \right ) R \left ( g_{\varrho i} \, p_{\!j} + g_{\varrho j} \, p_i - \frac {2}{3} \, \boldsymbol{g}_{\varrho } \boldsymbol{\cdot }\boldsymbol{p} \, \delta _{\textit{ij}} \right ) \nonumber \\ &\quad + \left (\frac {3}{2} \, \beta _\varLambda -\gamma _\varLambda \right ) R \, g_{\varrho k} \, r_{\textit{ijk}} \Bigg ] \end{align}
in terms of the orientation tensors of rank one (2.20), rank two (2.21) and rank three (2.22).
We note that the whole contribution (4.10) is equal to zero if the suspension is passive or at chemical equilibrium with
$\varrho =0$
and, thus,
$\boldsymbol{g}_\varrho =\boldsymbol{\nabla }\varrho =0$
. Otherwise, in active suspensions, this contribution can generate non-Newtonian effects characterised by non-vanishing normal stress differences
$N_1\ne 0$
or
$N_2\ne 0$
, as for the contribution (4.8) but here in the presence of the reaction. The first term going as
$\varrho q_{\textit{ij}}$
is proportional to the quadrupolar-order tensor (2.21). Its existence has been anticipated by Simha & Ramaswamy (Reference Simha and Ramaswamy2002) and further considered, e.g. by Hatwalne et al. (Reference Hatwalne, Ramaswamy, Rao and Simha2004) and Ramaswamy (Reference Ramaswamy2010). The other terms involve the concentration gradients and they would be equal to zero if the molecular concentrations were uniform. The second term depends on the polar-order vector (2.20) and the third term on the octupolar-order tensor (2.22).
The contribution (4.10) can be expressed in terms of the affinity of the chemical reaction (2.1), which is given as
${\mathcal A}\equiv (\mu _{{A}}-\mu _{{B}})/(k_{{ B}}T)=\ln [(\kappa _+^{{c}} \overline {c_{{A}}})/(\kappa _-^{{c}} \overline {c_{{B}}})]$
in terms of the chemical potentials
$\mu _\varkappa =\mu _{\varkappa 0}+k_{{B}}T \ln (\overline {c_\varkappa }/c_0)$
of the molecular species involved in the reaction, where
$c_0$
is a reference concentration value of one mole per litre. Accordingly, the reaction rate (3.4) can be written into the following equivalent form:
$\varrho =\kappa _+^{{c}} \overline {c_{{A}}} (1-\textrm { e}^{-{\mathcal A}})=\kappa _+^{{c}} \overline {c_{{A}}}{\mathcal A}+O({\mathcal A}^2)$
. In a regime close to chemical equilibrium where
${\mathcal A}\ll 1$
, the reaction rate is thus proportional to the chemical affinity, which is the non-equilibrium driving force of the reaction (de Groot & Mazur Reference de Groot and Mazur1984).
5. Case of chemically inert colloids with non-uniform diffusiophoretic constants
In the absence of chemical reaction, i.e. if the rate constants are equal to zero,
$\kappa _\pm ^{{c}}=0$
, the suspension is no longer active and the contribution (4.10) vanishes because
$\varrho =0$
and
$\boldsymbol{g}_\varrho =0$
. In such a passive suspension, the stress tensor may still depend on the molecular concentration gradients due to the contribution (4.8) if the diffusiophoretic constants are non-uniform on the surface of the colloidal particles, i.e. if
$\Delta b_\varkappa \ne 0$
.
The contribution (4.8) can be rewritten as
\begin{align} \Delta \sigma _{\textit{ij}}^{({d})} = \frac {27}{16 R} \, \frac {\phi \, \eta _0}{\displaystyle 1 + \frac {5b}{R}} \sum _\varkappa \Delta b_\varkappa \left ( g_{\varkappa i} \langle u_{\!j}\rangle + g_{\varkappa j} \langle u_i\rangle - \frac {2}{3} \, \boldsymbol{g}_\varkappa \boldsymbol{\cdot }\langle \boldsymbol{u}\rangle \, \delta _{\textit{ij}} + \frac {5}{6} \, g_{\varkappa k} \langle R_{\textit{ijk}} \rangle \right )\! , \end{align}
in terms of the local mean values (2.24) and the volume fraction
$\phi$
introduced in (4.5). This contribution is thus also proportional to the volume fraction as the modification of the shear viscosity (4.6).
If the particles were fully polarised in the
$x$
-direction with
$\boldsymbol{u}=(1,0,0)$
at some instant of the time evolution, the contribution (4.8) to the stress tensor would be given by
\begin{align} \Delta \boldsymbol{\sigma }^{({d})} = \left [\Delta \sigma _{\textit{ij}}^{({d})}\right ] = \frac {45}{32 R} \, \frac {\phi \, \eta _0}{\displaystyle 1 + \frac {5b}{R}} \sum _\varkappa \Delta b_\varkappa \left [ \begin{array}{ccc} 2 \, g_{\varkappa x} & g_{\varkappa y} & g_{\varkappa z} \\ g_{\varkappa y} & - g_{\varkappa x} & 0 \\ g_{\varkappa z} & 0 & - g_{\varkappa x} \\ \end{array} \right ]\! . \end{align}
The condition of incompressibility is satisfied, because
$\textrm {tr}\, \Delta \boldsymbol{\sigma }^{( { d})} =0$
. As a consequence, the pressure would change by
$\Delta P_{xx}^{({d})}=-\Delta \sigma _{xx}^{({d})}$
in the
$x$
-direction and by
$\Delta P_{yy}^{( { d})}=\Delta P_{zz}^{({d})}=({1}/{2})\, \Delta \sigma _{xx}^{({d})}$
in the orthogonal directions. To evaluate the contribution (5.2), we may assume that the Navier slip length is much smaller than the particle radius (i.e.
$b\simeq 0$
). (The condition
$b\ll R$
is indeed satisfied for common solid surfaces in contact with aqueous solutions (Joly et al. Reference Joly, Ybert and Bocquet2006; Bouzigues et al. Reference Bouzigues, Tabeling and Bocquet2008; Li et al. Reference Li, D’eramo, Lee, Monti, Yonger, Tabeling, Chollet, Bresson and Tran2015; Collis et al. Reference Collis, Olcum, Chakraborty, Manalis and Sader2021).) Supposing furthermore that the particles have a radius
$R=1\,\unicode{x03BC}$
m, occupy a volume fraction
$\phi =0.1$
in an aqueous solution of viscosity
$\eta _0=10^{-3}\,$
Pa s, and that typical experimental values are taken for the concentrations
$c_\varkappa \simeq 2\times 10^{27}\,$
m
$^{-3}$
, their relative gradients
$\Vert \boldsymbol{\nabla }c_\varkappa \Vert /c_\varkappa \sim 10^5\,$
m
$^{-1}$
and the diffusiophoretic constants such that
$\Delta b_\varkappa \sim 2\times 10^{-38}\,$
m
$^5$
s−1 (Ibrahim et al. Reference Ibrahim, Golestanian and Liverpool2017), the change of pressure would be of the order of
$\vert \Delta P_{xx}^{({d})}\vert \sim 10^{-3}$
Pa and it would depend on the magnitude and the direction of the concentration gradients.
6. Case of chemically active colloids with uniform diffusiophoretic constants
Here, we consider the case where the diffusiophoretic constants (2.8) take the same values on both hemispheres, i.e.
$\Delta b_\varkappa =0$
for
$\varkappa \in \{\textrm {A},\textrm {B}\}$
, so that the contribution (4.8) from simple diffusiophoresis is equal to zero. Moreover, in this case, the coefficients (C116), (C117) and (C118) of the contribution (4.10) from the reaction reduce to
according to (C104). The coefficients
$a_2$
,
$b_2$
and
$c_2$
depend on the Damköhler number (3.3), in a way which can be calculated using (A16)–(A18) for
$l=2$
. Next, they may be substituted into (4.10). Furthermore, the reaction rate (3.4) and its gradient (3.5) can be expressed in terms of the particle velocity generated by the reaction, which is known to take the following value (Gaspard & Kapral Reference Gaspard and Kapral2020) (this velocity was denoted
$V_{{sd}}\equiv V^{( {r})}$
by Gaspard & Kapral (Reference Gaspard and Kapral2020) to refer to its self-diffusiophoretic origin):
\begin{align} V^{({r})} = \frac {2 a_1}{\displaystyle 3\left (1+\frac {2b}{R}\right )} \, \varLambda \varrho , \end{align}
as summarised in Appendix E. An active particle is said to be a pusher if
$\varLambda \varrho \gt 0$
and a puller if
$\varLambda \varrho \lt 0$
(Reigh et al. Reference Reigh, Huang, Schofield and Kapral2017; Gaspard & Kapral Reference Gaspard and Kapral2019). (The notation
$\varUpsilon \equiv \varLambda \varrho$
is used by Gaspard & Kapral (Reference Gaspard and Kapral2019).)
Therefore, using the velocity (6.2) and the local mean values (2.24), the contribution (4.10) from the reaction to the stress tensor can be equivalently written as
\begin{align} \Delta \sigma _{\textit{ij}}^{({r})} &= \frac {27}{2 a_1} \, \frac {\displaystyle 1 + \frac {2b}{R}}{\displaystyle 1 + \frac {5b}{R}} \, \phi \, \eta _0 \, \bigg [\ a_2 \, \frac {V^{({r})}}{R} \, \langle Q_{\textit{ij}}\rangle \nonumber \\ &\quad + \frac {3}{10} \, (b_2 + c_2) \left ( \boldsymbol{\nabla} _i V^{({r})} \, \langle u_{\!j}\rangle + \boldsymbol{\nabla} _{\!j} V^{({r})} \, \langle u_i\rangle - \frac {2}{3} \, \boldsymbol{\nabla } V^{(\textrm { r})} \boldsymbol{\cdot }\langle \boldsymbol{u}\rangle \, \delta _{\textit{ij}} \right ) \nonumber \\&\quad + \left (\frac {3}{2} \, b_2 - c_2\right ) \boldsymbol{\nabla} _k V^{({r})} \, \langle R_{\textit{ijk}}\rangle \bigg ] , \end{align}
where
$\phi$
is the volume fraction of the Janus particles in the suspension, as introduced in (4.5). The first term with the coefficient
$a_2$
is proportional to the nematic-order parameter
$\langle Q_{\textit{ij}}\rangle$
and to the particle velocity (6.2), which is the non-equilibrium driving force caused by the reaction. The other terms with the coefficients
$b_2$
and
$c_2$
are proportional to the velocity gradient
\begin{align} \boldsymbol{\nabla }V^{ {(r)}} = \frac {2 a_1}{\displaystyle 3\left (1+\frac {2b}{R}\right )} \, \varLambda \boldsymbol{g}_\varrho = \frac {2 a_1 \varLambda }{\displaystyle 3\left (1+\frac {2b}{R}\right )} \left ( \kappa _+^{ {c}} \, \boldsymbol{\nabla } c_{{A}} - \kappa _-^{ {c}} \, \boldsymbol{\nabla } c_{{ B}} \right )\! , \end{align}
which is a further non-equilibrium driving force depending on the gradients of the molecular concentrations
$\boldsymbol{g}_{{A}}=\boldsymbol{\nabla } c_{{A}}$
and
$\boldsymbol{g}_{{B}}=\boldsymbol{\nabla } c_{{B}}$
. These other terms involve the polar- and octupolar-order parameters,
$\langle u_i\rangle$
and
$\langle R_{\textit{ijk}}\rangle$
, respectively.
We note that, under the same conditions, the overall reaction rate of a single Janus particle is given by (A25) in terms of the coefficients
$a_0$
and
$b_0$
; and its vectorial velocity by (E4) in terms of the coefficients
$a_1$
,
$b_1$
and
$c_1$
(Gaspard & Kapral Reference Gaspard and Kapral2020).
All the dimensionless coefficients
$a_0$
,
$a_1$
,
$a_2$
,
$b_0$
,
$b_1$
,
$b_2$
,
$c_1$
and
$c_2$
, which determine the macroscopic behaviour of the active suspension, can be calculated with (A16)–(A18) using Mathematica (Wolfram Reference Wolfram1988). Their values are given in table 1 and plotted in figures 2–4, as a function of the Damköhler number (3.3). In the diffusion-controlled regime, they decrease as
${\textit {Da}}^{-1}$
for
${\textit {Da}} \gg 1$
. We note that the coefficient
$a_2$
is negative and equal to zero at
${\textit {Da}}=0$
. Moreover, we see in figure 2 that the coefficient
$a_2$
goes to zero for
${\textit {Da}}\to 0$
in a way that is proportional to
$\textit {Da}$
. The fact that
$a_2({\textit {Da}}=0)=0$
is a consequence of the assumption that the catalytic and inert parts of the Janus particles are hemispheric. This would not be the case for Janus particles with non-symmetric catalytic and inert caps, for which
$a_2({\textit {Da}}=0)\ne 0$
(Michelin & Lauga Reference Michelin and Lauga2014). For hemispheric catalytic and inert caps, the other coefficients
$a_0$
,
$a_1$
,
$b_0$
,
$b_1$
,
$b_2$
,
$c_1$
and
$c_2$
take non-zero values at
${\textit {Da}}=0$
and, thus, in the reaction-limited regime for
${\textit {Da}}\ll 1$
, as shown in table 1 and figures 2–4.
Values of the coefficients
$a_0$
,
$a_1$
,
$a_2$
,
$b_0$
,
$b_1$
,
$b_2$
,
$c_1$
and
$c_2$
versus the Damköhler number (3.3). The coefficients
$a_0$
and
$b_0$
determine the overall reaction rate (A25);
$a_1$
,
$b_1$
and
$c_1$
the particle velocity (E4); and
$a_2$
,
$b_2$
and
$c_2$
the contribution (6.3) to the stress tensor. They are computed using (A16)–(A18) truncated to
$l^{\prime}_{{max}}=100$
.

Plots of the values of the coefficients
$a_0$
,
$a_1$
and
$a_2$
versus the Damköhler number (3.3). The coefficient
$a_0$
is involved in the overall reaction rate (A25),
$a_1$
in the particle velocity (E4) and
$a_2$
in the contribution (6.3) to the stress tensor. They determine the dependence of these quantities on the reaction rate
$\varrho$
. They are computed using (A16)–(A18) truncated to
$l^{\prime}_{{max}}=100$
.

Plots of the values of the coefficients
$b_0$
,
$b_1$
and
$b_2$
versus the Damköhler number (3.3). The coefficient
$b_0$
is involved in the overall reaction rate (A25),
$b_1$
in the particle velocity (E4) and
$b_2$
in the contribution (6.3) to the stress tensor. They determine the dependences of these quantities on the gradient
$\boldsymbol{g}_\varrho =\boldsymbol{\nabla }\varrho$
. They are computed using (A16)–(A18) truncated to
$l^{\prime}_{{max}}=100$
.

Plots of the values of the coefficients
$c_1$
and
$c_2$
versus the Damköhler number (3.3). The coefficient
$c_1$
is involved in the particle velocity (E4) and
$c_2$
in the contribution (6.3) to the stress tensor. They determine the dependences of these quantities on the gradient
$\boldsymbol{g}_\varrho =\boldsymbol{\nabla }\varrho$
. They are computed using (A16)–(A18) truncated to
$l^{\prime}_{{max}}=100$
.

The contribution (6.3) from the reaction can be compared with the contribution (4.5) to shear viscosity. They would have comparable values if the gradients
$G_{\textit{ij}}$
of the fluid velocity would be of the same order of magnitude as the particle velocity (6.2) divided by the particle radius
$R$
or its gradient components
$\boldsymbol{\nabla} _i V^{({r})}$
, and if the order parameters
$\langle Q_{\textit{ij}}\rangle$
or
$\langle u_i\rangle$
of the particle orientation were maximal. Generally speaking, (6.3) is the analogue of the force–velocity relation characterising active systems such as muscles.
To get an evaluation of the contribution (6.3), we assume that the terms with the components of the gradient
$\boldsymbol{\nabla }V^{({r})}$
are negligible and, again, that the Navier slip length is much smaller than the particle radius (i.e.
$b\simeq 0$
). In this case, we have that
because the particle velocity (6.2) is given by
$V^{({r})}\simeq (2/3) a_1 \varLambda \varrho$
in the same approximation.
Now, if the particles were fully polarised in the
$x$
-direction with
$\boldsymbol{u}=(1,0,0)$
, the activity would change the pressure by
$\Delta P_{xx}^{({r})}=-\Delta \sigma _{xx}^{({r})}\simeq -6 a_2 \phi \eta _0 R^{-1}\varLambda \varrho$
in that direction and by
$\Delta P_{yy}^{({r})}=\Delta P_{zz}^{(\textrm { r})}=-({1}/{2})\, \Delta P_{xx}^{({r})}$
in the orthogonal directions to maintain the incompressibility. Since
$a_2 \lt 0$
, the medium is extensile with
$\Delta P_{xx}^{({r})}\gt 0$
if the active particles are pushers and the medium is contractile with
$\Delta P_{xx}^{({r})}\lt 0$
if the active particles are pullers, which is consistent with expectations (Hatwalne et al. Reference Hatwalne, Ramaswamy, Rao and Simha2004; Ramaswamy Reference Ramaswamy2010). In an active suspension of Damköhler number equal to
${\textit {Da}}=10$
composed of particles of radius
$R=1\,\unicode{x03BC}$
m, moving with the speed
$V^{({r})}=10\,\unicode{x03BC}$
m s−1 in an aqueous solution of viscosity
$\eta _0=10^{-3}\,$
Pa s and occupying a volume fraction
$\phi =0.1$
, the activity would generate a change of pressure of about
$\Delta P_{xx}^{({r})} \simeq 1.7\times 10^{-3}$
Pa. This effect is proportional to the particle velocity
$V^{({r})}$
and to their volume fraction
$\phi$
.
Furthermore, we point out that the contribution (6.3) from the reaction to the stress tensor can generate the non-Newtonian effects aforementioned in § 4.4. However, the activity of the reaction does not modify the effective shear viscosity with respect to its known expression (4.6), because the particles are here assumed to be hydrodynamically spherical.
To recapitulate, note that, in the case where the colloids are chemically active and, moreover, have non-uniform diffusiophoretic constants (i.e.
$\Delta b_\varkappa \ne 0$
), the stress tensor includes the two contributions (4.8) and (4.10) in addition to the standard contribution (4.5) expected for hydrodynamically spherical particles.
7. Summary and perspectives
In this paper, the method pioneered by Einstein (Reference Einstein1906) (Landau & Lifshitz Reference Landau and Lifshitz1959; Batchelor Reference Batchelor1970) to calculate the stress tensor was extended from passive to active dilute suspensions composed of Janus colloidal particles that are propelled by self-diffusiophoresis and powered with chemical energy. In the dilute-system limit, the particles may be considered to be far apart and isolated from each other, so that their mutual interactions can be neglected. Under such circumstances, the stress tensor of the suspension can be deduced from the velocity and pressure of the fluid around a single Janus particle moving in coarse-grained background fields.
The particles were supposed to be spherical and made of two hemispheres: a catalytic hemisphere, where the chemical reaction powering the activity occurs, and the other one being non-catalytic. The propulsion of the particles was assumed to be induced by diffusiophoresis locally coupling the fluid velocity with the molecular concentration gradients generated by the surface chemical reaction on the catalytic hemisphere of the Janus particles. Accordingly, the molecular concentrations were first obtained before calculating the one-particle velocity field. The latter was expanded into vectorial spherical harmonics and corresponding radial functions, which were calculated by solving the steady Stokes equation and the condition of incompressibility. The coefficients of the radial functions were determined using the boundary conditions at the surface of the Janus particle and far away from the particle, together with the force-free and torque-free conditions.
As a consequence, the velocity and pressure disturbance fields of a single particle behave as
$\Delta \boldsymbol{v}=O(r^{-2})$
and
$\Delta p=O(r^{-3})$
, leading to three contributions to the stress tensor due to the presence of Janus particles in the suspension. In general, these contributions depend on the orientation of the Janus particles and their distribution. All these contributions to the stress tensor are traceless and symmetric in consistency with the incompressibility and torque-free conditions, respectively.
The first contribution leads to an enhancement of the shear viscosity, which is proportional to the volume fraction of suspended particles. This contribution is well known since the classic work of Einstein (Reference Einstein1906). Because the velocity field was here assumed to satisfy partial slip boundary conditions, the effective shear viscosity was obtained with its dependence on the Navier slip length, generalising Einstein’s result for no-slip spheres (Happel & Brenner Reference Happel and Brenner1983; Luo & Pozrikidis Reference Luo and Pozrikidis2008; Premlata & Wei Reference Premlata and Wei2022).
The second contribution is due to the simple diffusiophoresis of the particles and it depends on the gradients of molecular concentrations. It would persist in the absence of chemical reactions, i.e. if the suspension were passive. However, this contribution only exists if the diffusiophoretic constants take different values on the materials composing the two hemispheres of the Janus particles.
The third contribution constitutes the active part of the stress tensor, because it is proportional to the reaction rate and its gradient, which are equal to zero at chemical thermodynamic equilibrium. This contribution from the activity of the suspension also depends on the orientation of the Janus particles. Since the stress tensor is a rank-two tensor and the reaction rate density is a scalar, the corresponding term is proportional to the tensor characterising the nematic order of the Janus particles. Further terms coupling the gradient of the reaction rate with the vector of polar order and the rank-three tensor of octupolar order are also identified. These additional terms depend on the gradients of the molecular concentrations, as is the case for the second contribution due to simple diffusiophoresis, but now caused by the reaction.
These three contributions to the stress tensor are obtained in the limits of a dilute suspension (2.26), low Reynolds numbers (2.30), low Péclet numbers (2.31) and quasi-static responses in the low-frequency range where
$\omega \ll \eta _0/(\rho R^2)$
and
$\omega \ll D_\varkappa /R^2$
. In this respect, they are proportional to the volume fraction occupied by the suspended Janus particles, they are linear in the gradients of the velocity and concentration fields, and they are independent of time or frequency. Nevertheless, the systematic calculation achieved in this paper can predict the magnitude of these effects in terms of the microstructural properties of active suspensions for arbitrary values of the Damköhler number (3.3). These quantitative predictions can be tested experimentally in active suspensions where the molecular concentrations and their gradients are controlled. Furthermore, these effects are expected to become larger in denser suspensions. In this regard, the contributions of higher powers in the volume fraction could be obtained by extending the systematic methods developed by Bedeaux et al. (Reference Bedeaux, Kapral and Mazur1977) from passive to active suspensions. These methods could also be used to go beyond the quasi-static approximation and possibly obtain contributions to the stress tensor with time-dependent responses like for viscoelasticity (Zia & Brady Reference Zia and Brady2015). Furthermore, there exist methods to calculate the possible contributions from higher order in the Péclet numbers (Frankel & Khair Reference Frankel and Khair2014; Yariv & Kaynan Reference Yariv and Kaynan2017) and from their scaling behaviour at large Péclet numbers (Yariv & Michelin Reference Yariv and Michelin2015). For known micrometric active particles, the Reynolds numbers are significantly smaller than the Péclet numbers, so that the low-Reynolds-number limit should provide satisfactory predictions in the current circumstances.
To fulfil the programme of fluid mechanics for the active suspension of interest, the effective stress tensor we calculated should be inserted into the Navier–Stokes equation for the coarse-grained velocity field, which should next be solved to predict the flow of the suspension. In this regard, the Navier–Stokes equation must be supplemented with the advection–diffusion–reaction equations for the coarse-grained molecular concentrations and also the generalised local conservation (2.25) for the mean distribution function
$f_{{C}}(\boldsymbol{r},\boldsymbol{u},t)$
of the Janus particles. All these equations form a closed set of coupled partial differential equations to be solved with suitable boundary conditions to determine the time evolution of the active suspension of interest on macroscales larger than the mean distance between nearest neighbouring particles.
Similar programmes can also be carried out for active particles with other geometries like ellipsoidal active particles (Saintillan & Shelley Reference Saintillan and Shelley2015; Saintillan Reference Saintillan2018), spheroidal particles with a varying size of the catalytic patch (Popescu et al. Reference Popescu, Dietrich, Tasinkevych and Ralston2010; Michelin & Lauga Reference Michelin and Lauga2014) or sphere dimers with a catalytic sphere rigidly bound to a non-catalytic one (Reigh & Kapral Reference Reigh and Kapral2015). Already, if the colloidal particles are not hydrodynamically spherical because of a non-uniform Navier slip length on their surface, the effective stress tensor should include extra terms going as
$u_iu_kG_{jk}$
,
$u_{\!j}u_kG_{ik}$
or
$Q_{\textit{ij}}u_ku_lG_{kl}$
, which lead to non-Newtonian effects (Ramachandran & Khair Reference Ramachandran and Khair2009; Saintillan & Shelley Reference Saintillan and Shelley2013; Premlata & Wei Reference Premlata and Wei2022) and to corrections proportional to the nematic-order tensor
$q_{\textit{ij}}$
in the case of active particles like stick-slip squirmers (Yang et al. Reference Yang, Chen, Premlata and Wei2023). Furthermore, for active rodlike particles, similar contributions have been established leading to important changes in shear viscosity due to the activity of the suspension (Saintillan Reference Saintillan2010; Saintillan & Shelley Reference Saintillan and Shelley2013, Reference Saintillan and Shelley2015; López et al. Reference López, Gachelin, Douarche, Auradou and Clément2015; Saintillan Reference Saintillan2018).
In a more general perspective, the results based on chemohydrodynamics provide a fundamental justification to the analytic expressions commonly assumed for the stress tensor in the theories describing suspensions as continuous media. They also show that the study of synthetic active matter can open interesting perspectives to understand analogue phenomena in biological systems (Koch & Subramanian Reference Koch and Subramanian2011; Saintillan Reference Saintillan2018). In particular, the sarcomers composing muscles have a head-to-head structure with nematic ordering in the muscular fibres (Günther & Kruse Reference Günther and Kruse2007). In this analogy, the force exerted by muscles can be understood in terms of an active stress tensor similar to the one that is here obtained based on the physico-chemical principles. Along these lines, the present study may also lead to a better understanding of the active stress tensor for the mathematical description of biological tissues.
The present calculations leave open the issue of external forces or external torques exerted on the suspended particles. The presence of external forces should add a corresponding force density in the coarse-grained Navier–Stokes equation, while external torques should contribute by an antisymmetric part to the effective stress tensor (Landau & Lifshitz Reference Landau and Lifshitz1959; Batchelor Reference Batchelor1970).
Finally, we note that the mesoscopic approach developed in this paper is based on the knowledge of the constitutive properties of the solution and the fluid–solid interfaces. In turn, these properties should be determined by the microscopic dynamics of the atoms and molecules composing the system. Therefore, the fully microscopic derivation of the stress tensor in active suspensions will require first-principles approaches based on non-equilibrium statistical mechanics (Robertson et al. Reference Robertson, Schofield, Gaspard and Kapral2020, Reference Robertson, Schofield and Kapral2024; Gaspard Reference Gaspard2022).
Acknowledgements
The author thanks Raymond Kapral for many inspiring discussions, Joël Mabillard for reading the preprint and sharing useful comments, and David Gaspard for expertise support during this work.
Funding
This research was supported by the Université Libre de Bruxelles (ULB).
Declaration of interests
The author declares no competing interests.
Data availability statement
The numerical data reported in this study are available from the corresponding author upon reasonable request.
Appendix A. Stationary one-particle profiles of molecular concentrations
In this appendix, (2.33) with the boundary conditions (2.10) and (2.35) are solved to obtain the stationary molecular concentrations around one Janus particle located at
$\boldsymbol{r}=0$
and oriented vertically with
$\boldsymbol{u}=\boldsymbol{1}_z$
. The geometry of the particle and the spherical coordinates are shown in figure 3. The methods of Gaspard & Kapral (Reference Gaspard and Kapral2020) are used in the following calculation.
Since
$\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }=\partial _r$
because
$\boldsymbol{n}=\boldsymbol{1}_r=\boldsymbol{r}/r$
is the unit vector in the radial direction, the concentrations of molecular species
$\textrm {A}$
and
$\textrm {B}$
satisfy the following equations and boundary conditions:
where
$H^{{c}}(\theta )=1$
on the catalytic hemisphere and
$H^{ {c}}(\theta )=0$
on the non-catalytic one.
The molecular concentrations can be decomposed as
with
into the following functions:
obeying
where
We note that the function
$\phi$
introduced in (A5) is associated with the conservation of the total number
$N_0=N_{{A}}+N_{{B}}$
of solute molecules involved in the chemical reaction (2.1). The other function
$\psi$
is defined in such a way to satisfy (A7) that are decoupled from (A6) for
$\phi$
. In the theory of chemical reaction networks, there exist general methods to identify the conserved quantities, but the other quantities may have to satisfy nonlinear boundary conditions in the presence of autocatalysis or cross-catalysis (Gaspard Reference Gaspard2022).
Schematic representation of one Janus particle and the used spherical coordinates.

The solutions of the two problems (A6) and (A7) are given by
with
$\phi ^\circ \equiv \ell (D_{{A}}\, c_{{A}}^\circ + D_{{B}}\, c_{{B}}^\circ )$
and
in terms of the reaction rate (3.4), its gradient (3.5) and the function
$Z(r,\theta ,\varphi )$
that is the solution of the problem (3.2) as a consequence of (A7). (The notation used by Gaspard & Kapral (Reference Gaspard and Kapral2020) is
$(\varPhi )_{2020}=\phi$
,
$(\varPsi )_{2020}=\psi$
,
$(\varPsi _0)_{2020}=\ell ^2\varrho$
and
$({\mathcal F})_{2020}=\varrho ^{-1}Z$
.) The latter can be expanded as
\begin{align} Z(r,\theta ,\varphi ) = \sum _{lm} z_{lm} \left (\frac {R}{r}\right )^{l+1} Y_{lm}(\theta ,\varphi ) \end{align}
into the spherical harmonics defined by
\begin{eqnarray} && Y_{lm}(\theta ,\varphi ) =(-1)^m \sqrt {\frac {(2l+1)(l-m)!}{4\pi (l+m)!}}\, P_l^m(\cos \theta ) \, \textrm { e}^{ {i} m \varphi } \quad \mbox{and}\qquad \nonumber \\ && Y_{l,-m}(\theta ,\varphi ) = (-1)^m \,Y_{lm}^*(\theta ,\varphi ) \quad \mbox{for}\quad m\geqslant 0 \end{eqnarray}
in terms of the associated Legendre functions
$P_l^{\vert m \vert }(\xi )=(1-\xi ^2)^{\vert m \vert /2} ({ {\rm d}^{\vert m \vert }}/{\textrm {d}\xi ^{\vert m \vert }}) P_l(\xi )$
, where
$P_l(\xi )$
is the Legendre polynomial of degree
$l \geqslant \vert m \vert \geqslant 0$
(Bransden & Joachain Reference Bransden and Joachain2000).
The coefficients
$z_{lm}$
can be calculated by inserting the expansion (A11) into the boundary conditions of the problem (3.2) and by using the orthonormality of the spherical harmonics,
$\int Y_{lm}^* Y_{l'm'} \textrm {d}\varOmega =\delta _{ll'} \delta _{mm'}$
, leading to
where the coefficients are given by (the present quantities
$c_l$
,
${\mathcal C}_l$
and
$N_{ll'}$
are related to those used by Gaspard & Kapral (Reference Gaspard and Kapral2020) according to
$c_l=(c_l)_{2020}\sqrt {l(l+1)/2}$
,
${\mathcal C}_l=({\mathcal C}_l)_{2020}/\sqrt {2l(l+1)}$
and
$N_{ll'}=(N_{ll'})_{2020}/\sqrt {l(l+1)l'(l'+1)}$
. The other quantities are the same. This change in the definition of the coefficients
$c_l$
allows us to have similar expansions for
$\gamma _\varLambda$
in (C118) as for
$\alpha _\varLambda$
and
$\beta _\varLambda$
in (C116) and (C117).)
with
\begin{align} M_{ll'} &\equiv 2 \, \frac {l+1}{2l+1} \, \delta _{ll'} + {\textit {Da}} \, \int _0^1 \textrm {d}\xi \, P_l(\xi )\, P_{l'}(\xi ) \, \quad \mbox{and}\qquad \nonumber \\ N_{ll'} &\equiv \frac {2}{l(2l+1)} \, \delta _{ll'} + {\textit {Da}} \, \int _0^1 \textrm {d}\xi \, \frac {P_l^1(\xi )}{l(l+1)} \, \frac {P_{l'}^1(\xi )}{l'(l'+1)} . \\[-2pt] \nonumber \end{align}
We note that the function
$Z$
can also be decomposed as
into the functions
\begin{align} K &\equiv \sum _{l=0}^{\infty } a_l \, P_l(\cos \theta ) \left (\frac {R}{r}\right )^{l+1} , \\[-9pt] \nonumber \end{align}
\begin{align} L_z &\equiv \sum _{l=0}^{\infty } b_l \, P_l(\cos \theta ) \left (\frac {R}{r}\right )^{l+1} , \\[-9pt] \nonumber \end{align}
\begin{align} \frac {L_x}{\cos \varphi } &= \frac {L_y}{\sin \varphi } \equiv 2 \sum _{l=1}^{\infty } c_l \, \frac {P_l^1(\cos \theta )}{l(l+1)} \left (\frac {R}{r}\right )^{l+1} . \\[9pt] \nonumber \end{align}
The knowledge of these functions allows us to calculate the overall reaction rate (2.12) for a single Janus particle. Indeed, the molecular concentrations can be expressed by (A3) in terms of the functions defined in (A5). Using the boundary condition
$(\partial _r \phi )_R=0$
and the relation (A10), we obtain
Next, the expansion (A11) shows that
leading to the result
$W = \sqrt {4\pi } \, R^2 \, z_{00}$
, since
$\textrm {d}\varSigma = R^2 \,\textrm {d}\varOmega$
with
$\textrm {d}\varOmega =\sin \theta \,\textrm {d}\theta \,\textrm {d}\varphi$
and
$\int Y_{lm} \, \textrm {d}\varOmega = \sqrt {4\pi } \, \delta _{l0} \, \delta _{m0}$
. Using (A13) for the coefficient
$z_{00}$
, the overall reaction rate of a single Janus particle is finally given by
in terms of the reaction rate (3.4) and its gradient (3.5) (Gaspard & Kapral Reference Gaspard and Kapral2020). The second term involving the gradients of the molecular concentrations is negligible in front of the first under the conditions
$(R\Vert \boldsymbol{\nabla } c_\varkappa \Vert /c_\varkappa ) \ll 1$
, i.e. when the characteristic length scales of the gradients are larger than the radius
$R$
of the Janus particles. In the case where the gradients can be considered to be equal to zero (so that
$\boldsymbol{g}_\varrho =0$
), the expression (A25) for the overall reaction rate reduces to its first term
$W= 4\pi R^2 a_0 \varrho$
(Gaspard & Kapral Reference Gaspard and Kapral2018, Reference Gaspard and Kapral2019).
Appendix B. Vectorial spherical harmonics and their properties
The vectorial spherical harmonics are functions of the polar and azimuthal angles
$(\theta ,\varphi )$
. They are defined as follows in terms of the spherical harmonics
$Y_{lm}$
(Barrera et al. Reference Barrera, Estévez and Giraldo1985):
for
$l=0,1,2,\ldots$
and
$m=-l,-l+1,\ldots ,l-1,l$
, where
$\boldsymbol{1}_r\equiv \boldsymbol{r}/r$
. In the orthonormal basis of the units vectors
$\{\boldsymbol{1}_r,\boldsymbol{1}_\theta ,\boldsymbol{1}_\varphi \}$
associated with the spherical coordinates, the latter two vectorial spherical harmonics can be expressed as
showing that the three kinds of vectorial spherical harmonics are orthogonal to each other according to
Since
$Y_{00}=1/\sqrt {4\pi }$
, we have that
$\boldsymbol{\varPsi }_{00}=0$
and
$\boldsymbol{\varPhi }_{00}=0$
.
Furthermore, the vectorial spherical harmonics satisfy the relations
where
The divergences, the rotationals and the Laplacians of vectorial spherical harmonics multiplied by some radial function are given by the following formulae (Barrera et al. Reference Barrera, Estévez and Giraldo1985):
Appendix C. Stationary one-particle profile of the velocity
The expansion (3.9) of the velocity field inside and around a single Janus particle into vectorial spherical harmonics is determined by using successively the expression (2.2) for the velocity inside the particle, Stokes’ equation (2.32), the incompressibility condition (2.3), the boundary conditions (2.34) at large distances, the boundary condition (2.6) at
$r=R$
in the direction normal to the surface and the boundary conditions (2.7) at
$r=R$
in the two directions tangential to the surface. According to (B7)–(B9), the radial functions of the expansion (3.9) are in principle given by
C.1. Velocity inside the solid particle
Since the centre of the particle is located at
$\boldsymbol{r}=\boldsymbol{R}=0$
and moves with the rectilinear velocity
$\boldsymbol{V}$
, while the particle rotates with the angular velocity
$\boldsymbol{\varOmega }$
, the velocity field (2.2) inside the particle is given by
for
$r=\Vert \boldsymbol{r}\Vert \lt R$
. Using (C1)–(C3), this vector field can be expanded into vectorial spherical harmonics according to
\begin{align} \boldsymbol{v}_{\textit{solid}}= \sum _{m=-1}^{+1} \big ( V_{1m} \, \boldsymbol{\varUpsilon }_{1m} +V_{1m} \, \boldsymbol{\varPsi }_{1m} - r \, \varOmega _{1m} \, \boldsymbol{\varPhi }_{1m} \big ) \end{align}
with
We note that the particle velocity
$\boldsymbol{V}$
and its angular velocity
$\boldsymbol{\varOmega }$
will be determined here later by the force-free and torque-free conditions for the particle.
C.2. Stokes’ equation
Inserting the expansion (3.9) into Stokes’ equation (2.32) and using (B17)–(B19) together with the orthogonality of the vectorial spherical harmonics, we obtain the following equations for the corresponding radial functions:
The equations for
$f_{lm}$
and
$g_{lm}$
are coupled together, but they can be decoupled introducing the following functions:
They obey linear second-order inhomogeneous ordinary differential equations,
the solutions of which are given by
with some coefficients
$c_{lm}^{(\xi )}$
,
$d_{lm}^{(\xi )}$
,
$c_{lm}^{(\zeta )}$
and
$d_{lm}^{(\zeta )}$
to be determined.
Therefore, the corresponding radial functions have the following forms:
\begin{align} f_{lm}(r) &= \left [ \frac {l \, a_{lm}}{2(2l+1)} + (l+1) \, c_{lm}^{(\zeta )}\right ] r^{l+1} +\left [ \frac {(l+1)\, b_{lm}}{2(2l+1)} + l \, d_{lm}^{(\xi )}\right ] r^{-l} \nonumber \\ &\quad + \; l \, c_{lm}^{(\xi )} \, r^{l-1} + (l+1) \, d_{lm}^{(\zeta )} \, r^{-l-2} , \\[-9pt] \nonumber \end{align}
\begin{align} g_{lm}(r) &= \left [ \frac {a_{lm}}{2(2l+1)} - c_{lm}^{(\zeta )}\right ] r^{l+1} + \left [-\frac {b_{lm}}{2(2l+1)} + d_{lm}^{(\xi )}\right ] r^{-l} \nonumber \\ &\quad + \; c_{lm}^{(\xi )} \, r^{l-1} - d_{lm}^{(\zeta )} \, r^{-l-2} . \\[9pt] \nonumber \end{align}
Furthermore, (C10) implies that
with some coefficients
$x_{lm}$
and
$y_{lm}$
to be determined.
We note that the functions
$g_{00}(r)$
and
$h_{00}(r)$
are not needed since
$\boldsymbol{\varPsi }_{00}=0$
and
$\boldsymbol{\varPhi }_{00}=0$
.
C.3. Incompressibility condition
Inserting the expansion (3.9) into the incompressibility condition (2.3) and using (B11)–(B13) together with the orthogonality of the vectorial spherical harmonics, we get
so that the radial functions
$f_{lm}$
and
$g_{lm}$
should be related by
Inserting therein the solutions (C17) and (C18), we find that the two following series of coefficients are determined:
For
$l=0$
, the incompressibility condition implies that
$c_{00}^{(\zeta )}=0$
and
$b_{00}=0$
. As a consequence of incompressibility, the functions (C17) and (C18) thus become
If
$l=0$
, (C24) gives
$f_{00}(r) = \, d_{00}^{(\zeta )} r^{-2}$
, because
$b_{00}=0$
. Therefore, the coefficients
$a_{00}$
and
$c_{00}^{(\xi )}$
do not appear in the problem and they can be discarded.
C.4. Velocity field far from the particle
At large distances from the particle, the boundary conditions (2.34) should be satisfied. According to (2.27), the velocity field should thus behave as
$\boldsymbol{v} = \boldsymbol{v}^{\circ } + {\boldsymbol{\unicode{x1D642}}}\boldsymbol{\cdot }\boldsymbol{r}$
for
$r\to \infty$
, i.e. it may not grow faster than the first power of the radial distance
$r$
. According to (C24) and (C25), all the coefficients
$a_{lm}$
with
$l\geqslant 1$
and all the coefficients
$c_{lm}^{(\xi )}$
with
$l\geqslant 3$
should thus be equal to zero. Furthermore, according to (C19), all the coefficients
$x_{lm}$
with
$l\geqslant 2$
should also be equal to zero. Out of the coefficients
$c_{lm}^{(\xi )}$
and
$x_{lm}$
, there remain
$c_{1m}^{(\xi )}$
,
$c_{2m}^{(\xi )}$
and
$x_{1m}$
to be determined by the constant velocity
$\boldsymbol{v}^{\circ }$
and the velocity gradients
$\boldsymbol{\unicode{x1D642}}$
. Using (C1)–(C3), we find that
\begin{align} \boldsymbol{v}^{\circ } + {\boldsymbol{\unicode{x1D642}}}\boldsymbol{\cdot }\boldsymbol{r} = \sum _{m=-1}^{+1} c_{1m}^{(\xi )} \left ( \boldsymbol{\varUpsilon }_{1m} +\boldsymbol{\varPsi }_{1m}\right ) + r \left [\sum _{m=-2}^{+2} c_{2m}^{(\xi )} \left ( 2 \boldsymbol{\varUpsilon }_{2m} + \boldsymbol{\varPsi }_{2m}\right ) + \sum _{m=-1}^{+1} x_{1m}\, \boldsymbol{\varPhi }_{1m} \right ] \end{align}
with
\begin{align} c_{22}^{(\xi )} &= c_{2,-2}^{(\xi )*} = \sqrt {\frac {\pi }{30}} \left [ G_{xx}-G_{yy} - \textrm {i} \left ( G_{xy}+G_{yx} \right )\right ]\! , \nonumber \\ c_{21}^{(\xi )} &= -c_{2,-1}^{(\xi )*} = \sqrt {\frac {\pi }{30}} \left [-G_{xz}-G_{zx} + \textrm {i} \left ( G_{yz}+G_{zy} \right )\right ]\! , \nonumber \\ c_{20}^{(\xi )} &= -\sqrt {\frac {\pi }{5}} \left ( G_{xx} + G_{yy} \right )\! ; \\[-28pt] \nonumber \end{align}
We note that the vorticity at large distances from the particle is given by
\begin{eqnarray} \boldsymbol{\omega }_\infty = (\boldsymbol{\nabla }\times \boldsymbol{v})_\infty &=& (G_{zy}-G_{yz})\, \boldsymbol{1}_x + (G_{xz}-G_{zx})\, \boldsymbol{1}_y + (G_{yx}-G_{xy})\, \boldsymbol{1}_z \nonumber \\ &=& - 2 \sum _{m=-1}^{+1} x_{1m} \left (\boldsymbol{\varUpsilon }_{1m}+\boldsymbol{\varPsi }_{1m}\right ) . \end{eqnarray}
Therefore, the coefficients
$x_{1m}$
determine the vorticity far from the particle according to
$x_{1m}=-(\boldsymbol{\omega }_\infty )_{1m}/2$
.
C.5. Boundary condition normal to the particle surface
Now, we consider the boundary condition (2.6) for a single particle located at the origin
$\boldsymbol{r}=0$
of the reference frame, which reads
because the velocity field of the solid particle is given by (C4) and
$\boldsymbol{n}=\boldsymbol{r}/r$
is orthogonal to
$\boldsymbol{\varOmega }\times \boldsymbol{r}$
. The left-hand side of this boundary condition is the radial component
$ (\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{v} )_R = (\boldsymbol{1}_r\boldsymbol{\cdot }\boldsymbol{v} )_R =(v_r)_R$
of the velocity field (3.9) at
$r=R$
, since
$\boldsymbol{n}=\boldsymbol{1}_r$
. Its right-hand side is the radial component of the particle velocity,
$\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{V} =\boldsymbol{1}_r\boldsymbol{\cdot }\boldsymbol{V} =V_r$
. Using the expansions (3.9) and (C5) with the coefficients (C6), the boundary condition (C31) gives
\begin{align} (v_r)_R=\sum _{lm} f_{lm}(R) \, Y_{lm} = \sum _{m=-1}^{+1} V_{1m} \, Y_{1m} . \end{align}
Using (C24) with the already known coefficients
$b_{00}=0$
,
$a_{lm}=0$
for
$l\geqslant 1$
, and
$c_{lm}^{(\xi )}=0$
for
$l\geqslant 3$
, the boundary condition (C32) is satisfied if
from which we deduce that
Hence, the radial functions of the velocity field become
We note that the components of the particle velocity always appear in the combination
$V_{1m}-c_{1m}^{(\xi )}$
corresponding to the difference
$\boldsymbol{V}-\boldsymbol{v}^{\circ }$
with respect to the constant velocity of the fluid, which is a consequence of Galilean invariance.
The coefficients
$c_{1m}^{(\xi )}$
,
$c_{2m}^{(\xi )}$
and
$x_{1m}$
are given by (C27), (C28) and (C29). Thus, the coefficients
$V_{1m}$
,
$b_{lm}$
and
$y_{lm}$
should still be determined.
C.6. Boundary conditions tangential to the particle surface
In the two directions tangential to the interface between the particle and the fluid, the velocity field should satisfy the boundary conditions (2.7). For a single particle located at
$\boldsymbol{r}=0$
, they read
\begin{align} \left ({\boldsymbol{\unicode{x1D644}}}-\boldsymbol{n}\boldsymbol{n}\right ) \boldsymbol{\cdot }\left [\boldsymbol{v} -\boldsymbol{v}_{ {\textit{solid}}} - b\left (\boldsymbol{\nabla }\boldsymbol{v}+\boldsymbol{\nabla }\boldsymbol{v}^{{T}}\right )\boldsymbol{\cdot }\boldsymbol{n} +\sum _\varkappa b_\varkappa \boldsymbol{\nabla }c_\varkappa \right ]_R = 0 , \end{align}
where the velocity of the solid particle is given by (C4) and the diffusiophoretic terms coupling the velocity with the molecular concentrations by (3.1), the function
$Z(r,\theta ,\varphi )$
being expressed in terms of the coefficients (A16)–(A18) by (A11)–(A15) or, equivalently, by (A19)–(A22).
Since
$\boldsymbol{n}=\boldsymbol{1}_r$
, we have that
${\boldsymbol{\unicode{x1D644}}}-\boldsymbol{n}\boldsymbol{n}=\boldsymbol{1}_\theta \boldsymbol{1}_\theta + \boldsymbol{1}_\varphi \boldsymbol{1}_\varphi$
, which is the projection onto the plane tangential to the particle interface. All the vectors can be decomposed in the orthonormal basis
$\{ \boldsymbol{1}_r, \boldsymbol{1}_\theta , \boldsymbol{1}_\varphi \}$
. In particular, the velocity field can be equivalently expressed as
$\boldsymbol{v} = v_r\boldsymbol{1}_r+ v_\theta \boldsymbol{1}_\theta + v_\varphi \boldsymbol{1}_\varphi$
. In the two tangential directions, the boundary conditions (C50) give
\begin{align} \left (v_\theta \right )_R &= \boldsymbol{1}_\theta \boldsymbol{\cdot }\big ( \boldsymbol{V} + \boldsymbol{\varOmega }\times \boldsymbol{r}\big )_R + b \, \boldsymbol{1}_\theta \boldsymbol{\cdot }\left (\boldsymbol{\nabla }\boldsymbol{v}+\boldsymbol{\nabla }\boldsymbol{v}^{{T}}\right )_R\boldsymbol{\cdot }\boldsymbol{1}_r - \boldsymbol{1}_\theta \boldsymbol{\cdot }\left (\sum _\varkappa b_\varkappa \boldsymbol{\nabla }c_\varkappa \right )_R \!, \\[-9pt] \nonumber \end{align}
\begin{align} \left (v_\varphi \right )_R &= \boldsymbol{1}_\varphi \boldsymbol{\cdot }\big ( \boldsymbol{V} + \boldsymbol{\varOmega }\times \boldsymbol{r}\big )_R + b \, \boldsymbol{1}_\varphi \boldsymbol{\cdot }\left (\boldsymbol{\nabla }\boldsymbol{v}+\boldsymbol{\nabla }\boldsymbol{v}^{{T}}\right )_R\boldsymbol{\cdot }\boldsymbol{1}_r - \boldsymbol{1}_\varphi \boldsymbol{\cdot }\left (\sum _\varkappa b_\varkappa \boldsymbol{\nabla }c_\varkappa \right )_R \!. \\[9pt] \nonumber \end{align}
For the velocity of the solid particle, we have that
\begin{eqnarray} \left ({\boldsymbol{\unicode{x1D644}}}-\boldsymbol{n}\boldsymbol{n}\right ) \boldsymbol{\cdot }\left ( \boldsymbol{V} + \boldsymbol{\varOmega }\times \boldsymbol{r}\right )_R &=& \left (V_\theta + R \, \varOmega _\varphi \right ) \boldsymbol{1}_\theta + \left (V_\varphi - R \, \varOmega _\theta \right ) \boldsymbol{1}_\varphi \nonumber \\ &=& \sum _{lm} \left ( V_{1m} \, \delta _{l1} \, \boldsymbol{\varPsi }_{lm} - R \, \varOmega _{1m} \, \delta _{l1} \, \boldsymbol{\varPhi }_{lm} \right )\! , \end{eqnarray}
because of (C5).
According to (15.20) of Landau & Lifshitz (Reference Landau and Lifshitz1959) and using the expansion (3.9), the components of the symmetrised tensor of velocity gradients are given by
where the coefficients
take the following values:
Consequently, we have that
Taking the scalar product of the boundary conditions (C50) with the vectorial spherical harmonics
$\boldsymbol{\varPsi }_{lm}$
and
$\boldsymbol{\varPhi }_{lm}$
and using the relations (B8) and (B9), we infer that
\begin{align} g_{lm}(R) &= V_{1m} \, \delta _{l1} + b \, G_{lm} - X_{lm}^{(g)} \nonumber \\ &\qquad \mbox{with}\qquad X_{lm}^{(g)} = \frac {1}{l(l+1)} \int \textrm {d}\varOmega \; \boldsymbol{\varPsi }_{lm}^* \boldsymbol{\cdot }\left (\sum _\varkappa b_\varkappa \boldsymbol{\nabla }c_\varkappa \right )_R , \end{align}
\begin{align} h_{lm}(R) &= - R \, \varOmega _{1m} \, \delta _{l1} + b \, H_{lm} - X_{lm}^{(h)} \nonumber \\ &\qquad \mbox{with}\qquad X_{lm}^{(h)} = \frac {1}{l(l+1)} \int \textrm {d}\varOmega \; \boldsymbol{\varPhi }_{lm}^* \boldsymbol{\cdot }\left (\sum _\varkappa b_\varkappa \boldsymbol{\nabla }c_\varkappa \right )_R . \end{align}
As a consequence, we find that
\begin{align} b_{1m} &= - \frac {R}{\displaystyle 1 + \frac {3b}{R}} \left [ \frac {3}{2} \left ( 1 + \frac {2b}{R}\right ) \left ( c_{1m}^{(\xi )} - V_{1m} \right ) + X_{1m}^{(g)} \right ]\! , \\[-12pt] \nonumber \end{align}
\begin{align} b_{2m} &= - \frac {6\, R^2}{\displaystyle 1 + \frac {5b}{R}} \left [ \frac {5}{3} \left ( 1 + \frac {2b}{R}\right ) c_{2m}^{(\xi )} \, R + X_{2m}^{(g)} \right ]\! , \\[-12pt] \nonumber \end{align}
\begin{align} b_{lm} &= - \frac {l(2l-1)\, R^l}{\displaystyle 1 + \frac {(2l+1)b}{R}} \, X_{lm}^{(g)} \quad \mbox{for} \quad l \geqslant 3 , \\[-12pt] \nonumber \end{align}
\begin{align} y_{1m} &= - \frac {R^2}{\displaystyle 1 + \frac {3b}{R}} \left [ R \left ( x_{1m} + \varOmega _{1m} \right ) + X_{1m}^{(h)}\right ]\! , \\[-12pt] \nonumber \end{align}
\begin{align} y_{lm} &= - \frac {R^{l+1}}{\displaystyle 1 + \frac {(l+2)b}{R}} \, X_{lm}^{(h)} \quad \mbox{for} \quad l \geqslant 2 . \\[9pt] \nonumber \end{align}
In addition to the calculation of the coefficients (C63) and (C64), there remains to determine the consequences of the force-free and torque-free conditions.
C.7. Force-free condition
Using the results of Mazur & Bedeaux (Reference Mazur and Bedeaux1974) and Albano et al. (Reference Albano, Bedeaux and Mazur1975), the force (2.13) exerted by the fluid on a single Janus particle was shown by Gaspard & Kapral (Reference Gaspard and Kapral2018, Reference Gaspard and Kapral2019) to be given by
\begin{align} \boldsymbol{F} = - \frac {6\pi \eta _0 R}{\displaystyle 1 + \frac {3b}{R}} \left [\left ( 1 +\frac {2b}{R}\right ) \left (\boldsymbol{V}-\boldsymbol{v}^{\circ }\right ) - \sum _\varkappa \overline {b_\varkappa \left ({\boldsymbol{\unicode{x1D644}}}-\boldsymbol{n}\boldsymbol{n}\right ) \boldsymbol{\cdot }\boldsymbol{\nabla }c_\varkappa }^{ {s}} \right ]\! , \end{align}
where
denotes the average over the spherical surface of the particle. In particular, the surface averages of the vectorial spherical harmonics we here need can be calculated, leading to
Since
\begin{align} \left [\sum _\varkappa b_\varkappa \left ({\boldsymbol{\unicode{x1D644}}}-\boldsymbol{n}\boldsymbol{n}\right ) \boldsymbol{\cdot }\boldsymbol{\nabla }c_\varkappa \right ]_R = \sum _{lm} \left ( X_{lm}^{(g)} \, \boldsymbol{\varPsi }_{lm} + X_{lm}^{(h)} \, \boldsymbol{\varPhi }_{lm} \right ) \quad \mbox{with} \quad l \geqslant 1 , \end{align}
we have that
\begin{align} \sum _\varkappa \overline { b_\varkappa \left ({\boldsymbol{\unicode{x1D644}}}-\boldsymbol{n}\boldsymbol{n}\right ) \boldsymbol{\cdot }\boldsymbol{\nabla }c_\varkappa }^{{s}} = \frac {2}{3} \sum _{m=-1}^{+1} X_{1m}^{(g)} \left ( \boldsymbol{\varUpsilon }_{1m} + \boldsymbol{\varPsi }_{1m} \right ) . \end{align}
Now, the particle velocity can be expanded as
$\boldsymbol{V} = \sum _{m=-1}^{+1} V_{1m} ( \boldsymbol{\varUpsilon }_{1m} + \boldsymbol{\varPsi }_{1m} )$
according to (C4) and (C5) with
$\boldsymbol{\varOmega }=0$
. In addition, the constant velocity has a similar form
$\boldsymbol{v}^{\circ } = \sum _{m=-1}^{+1} c_{1m}^{(\xi )} ( \boldsymbol{\varUpsilon }_{1m} + \boldsymbol{\varPsi }_{1m} )$
because of (C26) with
$\boldsymbol{\unicode{x1D642}} = 0$
. Therefore, the force (C70) can be equivalently expressed as
\begin{align} \boldsymbol{F} = -4\pi \eta _0 \sum _{m=-1}^{+1} b_{1m} \left ( \boldsymbol{\varUpsilon }_{1m} + \boldsymbol{\varPsi }_{1m} \right ) \end{align}
in terms of the coefficients (C65).
The result is that the force-free condition
$\boldsymbol{F}=0$
is equivalent to the condition
$b_{1m}=0$
.
C.8. Torque-free condition
The torque (2.14) exerted by the fluid on the Janus particle was evaluated by Gaspard & Kapral (Reference Gaspard and Kapral2018, Reference Gaspard and Kapral2019) with the methods of Hills (Reference Hills1975) and Felderhof (Reference Felderhof1976a , Reference Felderhof b ), and is given by
\begin{align} \boldsymbol{T} = - \frac {8\pi \eta _0 R^3}{\displaystyle 1 + \frac {3b}{R}} \left [\left (\boldsymbol{\varOmega } - \frac {1}{2}\, \boldsymbol{\omega }_\infty \right ) - \frac {3}{2 R^2} \sum _\varkappa \overline {b_\varkappa \boldsymbol{r}\times \boldsymbol{\nabla }c_\varkappa }^{ {s}} \right ]\! , \end{align}
where
$\boldsymbol{\omega }_\infty =(\boldsymbol{\nabla }\times \boldsymbol{v})_\infty$
is the vorticity at large distances from the particle. Here, we have that
\begin{align} \frac {1}{R}\left (\sum _\varkappa b_\varkappa \boldsymbol{r}\times \boldsymbol{\nabla }c_\varkappa \right )_R = \sum _{lm} \left ( X_{lm}^{(g)} \, \boldsymbol{\varPhi }_{lm} - X_{lm}^{(h)} \, \boldsymbol{\varPsi }_{lm} \right ) \quad \mbox{with} \quad l \geqslant 1 , \end{align}
so that
\begin{align} \frac {1}{R}\sum _\varkappa \overline {b_\varkappa \boldsymbol{r}\times \boldsymbol{\nabla }c_\varkappa }^{ {s}} = - \frac {2}{3}\sum _{m=-1}^{+1} X_{1m}^{(h)} \left ( \boldsymbol{\varUpsilon }_{1m} + \boldsymbol{\varPsi }_{1m} \right )\! , \end{align}
as a consequence of (C72).
Since
\begin{align} \boldsymbol{\varOmega } - \frac {1}{2} \, \boldsymbol{\omega }_\infty &= \sum _{m=-1}^{+1} \left [\varOmega _{1m} - \frac {1}{2} \, (\boldsymbol{\omega }_\infty )_{1m}\right ] \left ( \boldsymbol{\varUpsilon }_{1m} + \boldsymbol{\varPsi }_{1m} \right ) \nonumber\\&= \sum _{m=-1}^{+1} \left (\varOmega _{1m} + x_{1m}\right ) \left ( \boldsymbol{\varUpsilon }_{1m} + \boldsymbol{\varPsi }_{1m} \right ) \end{align}
according to (C30), the torque (C76) has the following equivalent expression:
\begin{align} \boldsymbol{T} = 8\pi \eta _0 \sum _{m=-1}^{+1} y_{1m} \left ( \boldsymbol{\varUpsilon }_{1m} + \boldsymbol{\varPsi }_{1m} \right ) \end{align}
in terms of the coefficients (C68).
Here, the result is that the torque-free condition
$\boldsymbol{T}=0$
is equivalent to the condition
$y_{1m}=0$
.
C.9. Asymptotic velocity and pressure
Since all the coefficients are now determined, we can consider the asymptotic expressions for the velocity (3.9) and pressure (3.6) around a single Janus particle. Indeed, after having taken into account the boundary conditions on the surface of the particle and at infinity, we have obtained the radial functions (C41)–(C49) with the coefficients (C65)–(C69), where
$b_{1m}=0$
and
$y_{1m}=0$
because of the force-free and torque-free conditions. Therefore, the velocity and pressure fields can be expressed as supposed in (2.27) and (2.28), where the disturbance velocity and pressure behave according to (2.36) and (2.37) for
$r\to \infty$
.
On the one hand, the background velocity has the form given by (C26) as confirmed by the terms going as
$c_{1m}^{(\xi )}$
in (C42) and (C45), as
$c_{2m}^{(\xi )}\, r$
in (C43) and (C46), and as
$x_{1m}^{(\xi )}\, r$
in (C48). Next, the leading term of the disturbance velocity
$\Delta \boldsymbol{v}$
goes as
$1/r^2$
. Such a dependence on
$r$
is found in (C43) with the coefficient
$b_{2m}$
given by (C66). There is another such dependence in the radial function (C48), but its coefficient is
$y_{1m}=0$
because of the torque-free condition. Accordingly, the disturbance velocity has the following asymptotic profile:
\begin{align} \Delta \boldsymbol{v} = -\frac {R^2}{r^2} \, \frac {1}{\displaystyle 1 + \frac {5b}{R}} \sum _{m=-2}^{+2} \left [ 5 \left ( 1 + \frac {2b}{R}\right ) c_{2m}^{(\xi )} \, R + 3 \, X_{2m}^{(g)} \right ] \boldsymbol{\varUpsilon }_{2m} + O\left (\frac {1}{r^3}\right )\! , \end{align}
where the coefficients
$X_{2m}^{(g)}$
are given by (C63) with
$l=2$
and the coefficients
$c_{2m}^{(\xi )}$
by (C28).
On the other hand, the pressure is given by the expansion (3.6) into spherical harmonics with the radial functions (3.8), where
$a_{lm}=0$
for
$l\geqslant 1$
to obtain the asymptotic velocity (C26),
$b_{00}=0$
to satisfy the incompressibility condition (C23) for
$l=0$
, and
$b_{1m}=0$
for the force-free condition to hold. In (3.6), the background value of the pressure is given by
$p^{\circ }=p_{00} Y_{00}=\eta _0 a_{00} Y_{00}$
, since
$a_{00}$
is the coefficient of the constant profile going as
$r^0$
in (3.8) with
$l=0$
and because it can take an arbitrary value. Now, (C66) implies that the disturbance pressure has the following asymptotic expression:
\begin{align} \Delta p = -\frac {R^2}{r^3} \, \frac {2\, \eta _0}{\displaystyle 1 + \frac {5b}{R}} \sum _{m=-2}^{+2} \left [ 5 \left ( 1 + \frac {2b}{R}\right ) c_{2m}^{(\xi )} \, R + 3 \, X_{2m}^{(g)} \right ] Y_{2m} + O\left (\frac {1}{r^4}\right ) . \end{align}
The asymptotic profiles (C81) and (C82) are consistent with the expected behaviours (2.36) and (2.37) as a consequence of the force-free and torque-free conditions. These expressions are thus suitable to evaluate the effect due to the presence of Janus particles on the stress tensor using (2.38).
According to (3.1), the coefficients
$X_{2m}^{(g)}$
introduced in (C63) with
$l=2$
can be similarly decomposed as
into a contribution from simple diffusiophoresis and another one from the reaction,
We note that, in (C84), the projector
$ ({\boldsymbol{\unicode{x1D644}}}-\boldsymbol{n}\boldsymbol{n} )$
has been eliminated because the vectorial spherical harmonics
$\boldsymbol{\varPsi }_{2m}$
are oriented tangentially to the particle surface, i.e. inside the subspace of the projector.
In terms of the coefficients
$b_{2m}$
given by (C66), (C81) and (C82) read
\begin{align} \Delta \boldsymbol{v} = \frac {1}{2\, r^2} \sum _{m=-2}^{+2} b_{2m} \, Y_{2m} \, \boldsymbol{1}_r + O\left (\frac {1}{r^3}\right ) \quad \mbox{and}\quad \Delta p = \frac {\eta _0}{r^3} \sum _{m=-2}^{+2} b_{2m} \, Y_{2m} + O\left (\frac {1}{r^4}\right )\! , \end{align}
where the coefficients
$b_{2m}$
in (C86) can be decomposed as
into the three following contributions:
\begin{align} b_{2m}^{({v})} &\equiv -10 \, R^3 \, \frac {\displaystyle 1 + \frac {2b}{R}}{\displaystyle 1 + \frac {5b}{R}} \, c_{2m}^{(\xi )} , \end{align}
\begin{align} b_{2m}^{({d})} &\equiv -\frac {3}{2}\, \frac {R^2}{\displaystyle 1 + \frac {5b}{R}} \int \textrm {d}\varOmega \; \boldsymbol{\varPsi }_{2m}^* \boldsymbol{\cdot }\left ( b_{{A}} \boldsymbol{g}_{{A}} + b_{{B}} \boldsymbol{g}_{{B}}\right )\! , \end{align}
\begin{align} b_{2m}^{({r})} &\equiv -\frac {R^3}{\displaystyle 1 + \frac {5b}{R}} \int \textrm {d}\varOmega \; \boldsymbol{\varPsi }_{2m}^* \boldsymbol{\cdot }\left [\left (\frac {b_{{B}}}{D_{{B}}} - \frac {b_{{A}}}{D_{{A}}}\right ) (\boldsymbol{\nabla } Z)_R\right ]\! . \end{align}
Since the coefficients
$c_{2m}^{(\xi )}$
are related to the tensor of velocity gradients according to (C28), the coefficients
$b_{2m}^{({v})}$
form a contribution that modifies the shear viscosity due to the presence of colloidal particles. The coefficients
$b_{2m}^{({d})}$
are associated with the contribution from simple diffusiophoresis and the coefficients
$b_{2m}^{({r})}$
from reaction. The corresponding velocity and pressure fields will be calculated separately in the following Appendices C.10, C.11 and C.12, respectively.
C.10. Disturbance velocity and pressure due to the velocity gradients
According to (C26) with
$\boldsymbol{v}^{\circ }=0$
and the identities
$\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\varUpsilon }_{2m}=Y_{2m}$
,
$\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\varPsi }_{2m}=0$
and
$\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\varPhi }_{2m}=0$
, we have that
\begin{align} \sum _{m=-2}^{+2} c_{2m}^{(\xi )} \, Y_{2m} = \frac {1}{2} \, G_{\textit{ij}} \, n_i \, n_{\!j} = \frac {1}{2} \, G_{\textit{ij}} \, \frac {r_i \, r_{\!j}}{r^2} , \end{align}
because
$\boldsymbol{n}=\boldsymbol{r}/r$
, i.e.
$n_i=r_i/r$
in Cartesian coordinates. Therefore, the disturbance velocity and pressure (C86) due to the contribution of the coefficients (C88) can be expressed as
\begin{eqnarray} &&\Delta v_i^{({v})} = A_{jk}^{({v})} \, \frac {r_i \, r_{\!j} \, r_k}{r^5} \quad \mbox{and}\quad \Delta p^{({v})} = 2\eta _0 \, A_{jk}^{({v})} \, \frac {r_{\!j} \, r_k}{r^5} \nonumber \\ &&\qquad \mbox{with}\quad A_{jk}^{({v})} \equiv - \frac {5}{2} \, \frac {\displaystyle 1 + \frac {2b}{R}}{\displaystyle 1 + \frac {5b}{R}} \, R^3 \, G_{jk} . \end{eqnarray}
Since the tensor of velocity gradients is traceless because of the incompressibility condition, we have that
$\textrm {tr}\, \boldsymbol{\unicode{x1D63C}}^{({v})}=A_{\textit{jj}}^{({v})}=0$
. This contribution goes as the tensor of velocity gradients and provides the well-known correction to the shear viscosity
$\eta _0$
of the fluid due to the presence of particles in the suspension (Einstein Reference Einstein1906; Happel & Brenner Reference Happel and Brenner1983; Luo & Pozrikidis Reference Luo and Pozrikidis2008; Premlata & Wei Reference Premlata and Wei2022).
C.11. Disturbance velocity and pressure due to simple diffusiophoresis
Since the diffusiophoretic constants of the molecular species may take different values on the two hemispheres according to (2.8), the coefficients (C84) are given by
Since the two hemispheres are mapped onto each other by space inversion
$(r,\theta ,\varphi )\to (r,\pi -\theta ,\varphi +\pi )$
and the spherical harmonics satisfy
$Y_{lm}(\pi -\theta ,\varphi +\pi )=(-1)^l \, Y_{lm}(\theta ,\varphi )$
, the vectorial spherical harmonics (B2) change by
$\boldsymbol{\varPsi }_{lm}(\pi -\theta ,\varphi +\pi )=(-1)^{l+1} \, \boldsymbol{\varPsi }_{lm}(\theta ,\varphi )$
. Therefore,
$\boldsymbol{\varPsi }_{2m}$
is odd under space inversion and (C93) becomes
The calculation of these coefficients with Mathematica (Wolfram Reference Wolfram1988) shows that
\begin{align}& X_{21}^{(g, {d})} = -X_{2,-1}^{(g,{d})*} = \frac {1}{16}\sqrt {\frac {15 \pi }{2}} \sum _\varkappa \Delta b_\varkappa \left ( -g_{\varkappa x} + \textrm {i} \, g_{\varkappa y}\right )\! , \end{align}
Now, we can use the following expressions for the spherical harmonics of order two with
$m=+1,0,-1$
:
and the assumption that the Janus particle is oriented in the direction
$\boldsymbol{u}=\boldsymbol{1}_z$
to obtain
\begin{align} \sum _{m=-2}^{+2} X_{2m}^{(g, {d})} \, Y_{2m} &= \frac {15}{64} \sum _\varkappa \Delta b_\varkappa \big \{ (\boldsymbol{g}_\varkappa \boldsymbol{\cdot }\boldsymbol{u})\big [ 3 (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})^2-1\big ]\big . \nonumber\\&\quad\big . + 2 \big [\boldsymbol{g}_\varkappa \boldsymbol{\cdot }\boldsymbol{n} - (\boldsymbol{g}_\varkappa \boldsymbol{\cdot }\boldsymbol{u})(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})\big ] (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})\big \}\! . \end{align}
Since
$\boldsymbol{n}=\boldsymbol{r}/r=(r_i/r)_{i=x,y,z}$
, the disturbance velocity and pressure due to simple diffusiophoresis are given by
\begin{align} \Delta v_i^{({d})} = A_{jk}^{({d})} \, \frac {r_i \, r_{\!j} \, r_k}{r^5} + B^{({d})} \, \frac {r_i}{r^3} \quad \mbox{and}\quad \Delta p^{({d})} = 2\eta _0 \left ( A_{jk}^{({d})} \, \frac {r_{\!j} \, r_k}{r^5} + \frac {B^{({d})}}{r^3}\right ) \end{align}
with
\begin{align} A_{jk}^{({d})} &\equiv - \frac {3\, R^2}{\displaystyle 1 + \frac {5b}{R}} \, \frac {15}{64} \sum _\varkappa \Delta b_\varkappa \left ( 2\, g_{\varkappa j} + \, \boldsymbol{g}_\varkappa \boldsymbol{\cdot }\boldsymbol{u}\, u_{\!j} \right ) u_k , \end{align}
\begin{align} B^{({d})} &\equiv \frac {3\, R^2}{\displaystyle 1 + \frac {5b}{R}} \, \frac {15}{64} \sum _\varkappa \Delta b_\varkappa \left (\boldsymbol{g}_\varkappa \boldsymbol{\cdot }\boldsymbol{u}\right ) \!. \end{align}
We note that
$\textrm {tr}\, \boldsymbol{\unicode{x1D63C}}^{({d})}=A_{\textit{jj}}^{({d})}=-3 B^{({d})}$
, because
$\boldsymbol{u}$
is a unit vector, so that
which has the form of a so-called stresslet (Batchelor Reference Batchelor1970).
C.12. Disturbance velocity and pressure due to reaction
Again, the diffusiophoretic constants taking different values on the two hemispheres according to (2.8), we can write
\begin{align} \frac {b_{{B}}}{D_{{B}}} - \frac {b_{{A}}}{D_{{A}}} = \sum _{h= {c}, {n}} \varLambda ^h \, H^h(\theta ) \quad \mbox{with}\; \varLambda ^h \equiv \frac {b_{{B}}^h}{D_{{B}}} - \frac {b_{{A}}^h}{D_{{A}}} . \end{align}
Therefore, the coefficients (C85) can be calculated with
Since the function
$Z$
has the expansion (A11), its gradient is given by
We thus have that
because of
$\boldsymbol{\varPsi }_{2m}^* \boldsymbol{\cdot }\boldsymbol{\varUpsilon }_{l'm'}=0$
and (B4). Using the property
$Y_{lm}(\pi -\theta ,\varphi +\pi )=(-1)^l \, Y_{lm}(\theta ,\varphi )$
for the spherical harmonics and the formula
$\int _0^{2\pi } \textrm {e}^{- {i}m\varphi } \textrm {e}^{ {i} m'\varphi } \,\textrm {d}\varphi =\delta _{mm'}$
, we find that the coefficients can be expressed as
with the integrals
These integrals can be evaluated using Mathematica (Wolfram Reference Wolfram1988), leading to
\begin{align} X_{20}^{(g, {r})} &= \frac {\varLambda ^{ {c}}+\varLambda ^{ {n}}}{2} \, z_{20} + \left (\varLambda ^{ {c}}-\varLambda ^{ {n}}\right ) \nonumber \\ & \quad\times \left ( \frac {\sqrt {15}}{16} \, z_{10} + \frac {\sqrt {35}}{16} \, z_{30} - \frac {5\sqrt {55}}{256} \, z_{50} + \frac {7\sqrt {3}}{128} \, z_{70} - \frac {15\sqrt {95}}{2048} \, z_{90} + {\cdots} \right )\! , \\[-12pt] \nonumber \end{align}
\begin{align} X_{21}^{(g, {r})} &= \frac {\varLambda ^{ {c}}+\varLambda ^{{n}}}{2} \, z_{21} + \left (\varLambda ^{ {c}}-\varLambda ^{ {n}}\right ) \nonumber \\ & \quad\times \left ( \frac {\sqrt {5}}{16} \, z_{11} + \frac {1}{8}\sqrt {\frac {35}{2}} \, z_{31} - \frac {25\sqrt {11}}{256} \, z_{51} + \frac {7\sqrt {7}}{64} \, z_{71} - \frac {75\sqrt {57}}{2048} \, z_{91} + {\cdots} \right )\! , \\[-12pt] \nonumber \end{align}
where the last line results from the facts that
$z_{l,-1}=-z_{l1}^*$
and
$z_{l,\pm 2}=0$
according to (A14) and (A15). Now, using the values of the coefficients
$z_{lm}$
given by (A13)–(A15) in terms of the coefficients (A16) with (A17) and (A18), we obtain
where
\begin{align} \alpha _\varLambda &\equiv \frac {\varLambda ^{ {c}}+\varLambda ^{ {n}}}{2} \, a_2 \nonumber \\ &\quad + \left (\varLambda ^{ {c}}-\varLambda ^{ {n}}\right ) \left ( \frac {5}{16} \, a_1 + \frac {5}{16} \, a_3 - \frac {25}{256} \, a_5 + \frac {7}{128} \, a_7 - \frac {75}{2048} \, a_9 + {\cdots} \right )\! , \end{align}
\begin{align} \beta _\varLambda &\equiv \frac {\varLambda ^{{c}}+\varLambda ^{ {n}}}{2} \, b_2 \nonumber \\ &\quad + \left (\varLambda ^{ {c}}-\varLambda ^{ {n}}\right ) \left ( \frac {5}{16} \, b_1 + \frac {5}{16} \, b_3 - \frac {25}{256} \, b_5 + \frac {7}{128} \, b_7 - \frac {75}{2048} \, b_9 + {\cdots} \right )\! , \end{align}
\begin{align} \gamma _\varLambda &\equiv \frac {\varLambda ^{ {c}}+\varLambda ^{ {n}}}{2} \, c_2 \nonumber \\ &\quad + \left (\varLambda ^{ {c}}-\varLambda ^{ {n}}\right ) \left ( \frac {5}{16} \, c_1 + \frac {5}{16} \, c_3 - \frac {25}{256} \, c_5 + \frac {7}{128} \, c_7 - \frac {75}{2048} \, c_9 + {\cdots} \right )\! . \end{align}
Using (C98) and the assumption
$\boldsymbol{u}=\boldsymbol{1}_z$
, we find that
\begin{eqnarray} \sum _{m=-2}^{+2} X_{2m}^{(g, {r})} \, Y_{2m} &=& \frac {1}{2} \left ( \alpha _\varLambda \, \varrho + \frac {3R}{2} \, \beta _\varLambda \, \boldsymbol{g}_{\varrho }\boldsymbol{\cdot }\boldsymbol{u} \right )\big [ 3 (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})^2-1\big ] \nonumber \\&&+ \frac {3R}{2} \, \gamma _\varLambda \left [\boldsymbol{g}_{\varrho }\boldsymbol{\cdot }\boldsymbol{n}-(\boldsymbol{g}_{\varrho }\boldsymbol{\cdot }\boldsymbol{u})(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})\right ] (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}) . \end{eqnarray}
Consequently, the disturbance velocity and pressure due to reaction are given by
\begin{align} \Delta v_i^{({r})} = A_{jk}^{({r})} \, \frac {r_i \, r_{\!j} \, r_k}{r^5} + B^{({r})} \, \frac {r_i}{r^3} \quad \mbox{and}\quad \Delta p^{({r})} = 2\eta _0 \left ( A_{jk}^{({r})} \, \frac {r_{\!j} \, r_k}{r^5} + \frac {B^{( { r})}}{r^3} \right ) \end{align}
with
\begin{align} A_{jk}^{({r})} &\equiv - \frac {3\, R^2}{\displaystyle 1 + \frac {5b}{R}} \, \frac {3}{2} \left [\left ( \alpha _\varLambda \, \varrho + \frac {3R}{2} \, \beta _\varLambda \, \boldsymbol{g}_{\varrho }\boldsymbol{\cdot }\boldsymbol{u}\right ) u_{\!j} + R \, \gamma _\varLambda \left ( g_{\varrho j} - \boldsymbol{g}_{\varrho }\boldsymbol{\cdot }\boldsymbol{u} \, u_{\!j} \right ) \right ] u_k , \end{align}
\begin{align} B^{({r})} &\equiv \frac {3\, R^2}{\displaystyle 1 + \frac {5b}{R}} \, \frac {1}{2} \left ( \alpha _\varLambda \, \varrho + \frac {3R}{2} \, \beta _\varLambda \, \boldsymbol{g}_{\varrho }\boldsymbol{\cdot }\boldsymbol{u} \right )\! . \end{align}
Here also, we have the property that
$\textrm {tr}\, \boldsymbol{\unicode{x1D63C}}^{({r})}=A_{\textit{jj}}^{({ r})}=-3 B^{({r})}$
, because
$\boldsymbol{u}$
is a unit vector, so that this contribution is given by the following stresslet:
Appendix D. Calculation of the change in the stress tensor
Here, we calculate the change (4.2) of the stress tensor due to the disturbance velocity and pressure (3.10). The latter can be equivalently written as
up to terms having higher powers of
$r^{-1}$
and leading to negligible corrections in the limit
$r\to \infty$
, i.e. in the dilute-system limit. Inserting the expressions (D1) into (4.2), we get
The derivatives
$\boldsymbol{\nabla} _i\varXi =\partial \varXi /\partial r_i$
of the function
$\varXi$
can be calculated. Afterwards, the function
$\varXi$
and its derivatives can be expressed in terms of the unit vector
$n_i=r_i/r$
, giving
and a similar expression for
$\boldsymbol{\nabla} _k\varXi$
. Moreover, we have that
$r_i = r \, n_i$
and
$\textrm {d}\varSigma _i = r^2 \, n_i \, \textrm {d}\varOmega$
. Therefore, (D2) becomes
with the notation
$\overline {(\boldsymbol{\cdot })} \equiv (4\pi )^{-1} \int (\boldsymbol{\cdot }) \, \textrm {d}\varOmega$
. Since
we finally obtain the result given by (4.3).
Appendix E. Particle velocities in the case of uniform diffusiophoretic constants
The force-free and torque-free conditions fully determine the rectilinear and angular velocities of the Janus particles in the active suspension. Indeed, using the force-free condition
$\boldsymbol{F}=0$
with (C70) and the torque-free condition
$\boldsymbol{T}=0$
with (C76), we obtain the following expressions for the particle velocities:
\begin{align} \boldsymbol{V} &= \boldsymbol{v}^{\circ } + \frac {1}{\displaystyle 1 +\frac {2b}{R}} \sum _\varkappa \overline {b_\varkappa \left ({\boldsymbol{\unicode{x1D644}}}-\boldsymbol{n}\boldsymbol{n}\right ) \boldsymbol{\cdot }\boldsymbol{\nabla }c_\varkappa }^{ {s}} , \end{align}
These velocities have been calculated by Gaspard & Kapral (Reference Gaspard and Kapral2020) for diffusiophoretic constants taking different values on the two hemispheres according to (2.8).
In the case
$b_\varkappa ^{{c}}=b_\varkappa ^{ {n}}$
, where the diffusiophoretic constants take uniform values, the general expressions obtained by Gaspard & Kapral (Reference Gaspard and Kapral2020) can be simplified using the following relations between the quantities defined therein and those used in the present article:
\begin{eqnarray} &&\alpha ^{ {c}} + \alpha ^{ {n}} = -\frac {4}{3}\, a_1 , \quad \beta ^{ {c}} + \beta ^{ {n}} = -2 \, b_1 , \quad \gamma ^{ {c}} + \gamma ^{ {n}} = -2\, c_1 , \quad \delta ^{ {c}} + \delta ^{{n}} = 0 , \nonumber \\&& \mbox{and}\quad Y^{ {c}} = Y^{ {n}} = - \varLambda , \end{eqnarray}
where the coefficients
$a_1$
,
$b_1$
and
$c_1$
are given by (A16), (A17) and (A18), and
$\varLambda$
is introduced in (6.1). Accordingly, the rectilinear velocity of a Janus particle can be expressed as
\begin{align} \boldsymbol{V} = \boldsymbol{v}^{\circ } + \frac {1}{\displaystyle 1 +\frac {2b}{R}} \left [ b_{{A}} \, \boldsymbol{g}_{{A}} + b_{{B}} \, \boldsymbol{g}_{{B}} + \frac {2}{3} \, a_1 \varLambda \, \varrho \, \boldsymbol{u} + R \, c_1 \, \varLambda \, \boldsymbol{g}_\varrho + R \, (b_1 - c_1) \, \varLambda \, \boldsymbol{g}_\varrho \boldsymbol{\cdot }\boldsymbol{u} \, \boldsymbol{u} \right ] \end{align}
with
$\boldsymbol{v}^{\circ }=\overline {\boldsymbol{v}}$
(Gaspard & Kapral Reference Gaspard and Kapral2020). In the absence of molecular concentration gradients, i.e. if
${\boldsymbol{g}}_{{A}}={\boldsymbol{g}}_{ {B}}=0$
, the vectorial particle velocity (E4) becomes
$\boldsymbol{V} = \boldsymbol{v}^{\circ } + V^{({r})} \, \boldsymbol{u}$
, as written in terms of the scalar particle velocity (6.2) (Reigh et al. Reference Reigh, Huang, Schofield and Kapral2017; Gaspard & Kapral Reference Gaspard and Kapral2018, Reference Gaspard and Kapral2019).
Furthermore, in the case of uniform diffusiophoretic constants, the angular velocity (E2) of a Janus particle reduces to
with a contribution from diffusiophoresis equal to zero (Gaspard & Kapral Reference Gaspard and Kapral2018, Reference Gaspard and Kapral2019, Reference Gaspard and Kapral2020).









































