2. Formulation of the problem

We aim to calculate the mobility matrix for two spherical particles of equal radius $a$ with relative distance $s$ between the centres of both particles, which are suspended in a fluid with position-dependent viscosity $\eta (\boldsymbol {r})$. We assume that the fluid viscosity is given by

(2.1)\begin{equation} \eta(\boldsymbol{r}) = \eta^{(0)} + \epsilon\,\eta^{(1)}(\boldsymbol{r}), \end{equation}

where $\eta ^{(0)}$ is the reference viscosity, $\epsilon < 1$ is the perturbation parameter, and $\boldsymbol {r}$ is the position in the fluid. Upper indices in brackets correspond to orders in $\epsilon$. We assume that the viscosity perturbation is at most of the order of $\eta ^{(0)}$, such that $\epsilon$ corresponds to the ratio of the magnitudes of viscosity perturbation and the constant-viscosity background. The Reynolds number (Purcell Reference Purcell1977) associated with the particles is assumed to be zero, such that the fluid is described by the incompressible Stokes equations (Kim & Karilla Reference Kim and Karilla2005, p. 9)

(2.2a,b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma} (\boldsymbol{r}) = 0, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} (\boldsymbol{r}) = 0. \end{equation}

Here, $\boldsymbol {u} (\boldsymbol {r})$ denotes the fluid velocity, $\boldsymbol {\sigma } (\boldsymbol {r})$ is the stress tensor defined as

(2.3)\begin{equation} \boldsymbol{\sigma} (\boldsymbol{r}) ={-} p(\boldsymbol{r})\,{\boldsymbol{\mathsf{I}}} + 2\,\eta(\boldsymbol{r})\, {\boldsymbol{\mathsf{E}}} (\boldsymbol{r}), \end{equation}

and the bold centred dot denotes contraction of the last index of the preceding quantity with the first index of the following quantity, as scalar products of vectors and matrix–vector products. Also, $p (\boldsymbol {r})$ denotes the pressure in the fluid, ${\boldsymbol{\mathsf{E}}}(\boldsymbol {r}) = (\boldsymbol {\nabla } \boldsymbol {u} (\boldsymbol {r}) + (\boldsymbol {\nabla } \boldsymbol {u}(\boldsymbol {r}))^{\rm T} )$/2 is the strain rate in the fluid, and ${\boldsymbol{\mathsf{I}}}$ is the $3 \times 3$ identity matrix. The symbol ^T denotes matrix transposition.

Due to the linearity of the Stokes equations, there exists a linear relation between the (angular) velocities ($\boldsymbol {U}_b$ and $\boldsymbol {\varOmega }_b$) of each particle and the hydrodynamic forces/torques ($\boldsymbol {F}^{hydro}_b$ and $\boldsymbol {T}^{hydro}_b$) that the fluid exerts on each particle $b = 1, 2$. To shorten notation for the subsequent calculations, we introduce generalised velocity and force vectors:

(2.4a,b)\begin{equation} \boldsymbol{\mathcal{U}} = (\boldsymbol{U}_1, \boldsymbol{U}_2, \boldsymbol{\varOmega}_1, \boldsymbol{\varOmega}_2), \quad \boldsymbol{\mathcal{F}}^{hydro} = (\boldsymbol{F}^{hydro}_1, \boldsymbol{F}^{hydro}_2, \boldsymbol{T}^{hydro}_1, \boldsymbol{T}^{hydro}_2). \end{equation}

The components of both vectors as well as of the resistance and mobility matrix (${\boldsymbol{\mathsf{R}}}$ and $\boldsymbol {\mu }$) defined below will be referred to by upper indices $P$, $Q$, which adopt the values $t$ (translation) or $r$ (rotation), and lower indices $b$, $c$ denoting the particle number (1 or 2). Thereby, significant brevity of the subsequent equations is achieved. We note that this notation is different to a common alternative where $F$, $U$, $T$ and $\varOmega$ are used as indices for the translational and rotational degrees of freedom. The conversion between the two notations is straightforward, e.g. $\boldsymbol {\mathcal {U}}^t_1 = \boldsymbol {U}_1$, $(\boldsymbol {\mathcal {F}}^{hydro})^r_2 = \boldsymbol {T}^{hydro}_2$ and ${\boldsymbol{\mathsf{R}}}^{tr}_{11} = {\boldsymbol{\mathsf{R}}}^{F \varOmega }_{11}$. The resistance matrix with components ${\boldsymbol{\mathsf{R}}}^{PQ}_{bc}$ is then defined by

(2.5)\begin{equation} (\boldsymbol{\mathcal{F}}^{hydro})^P_b =: \sum_{Q, c}{\boldsymbol{\mathsf{R}}}^{PQ}_{bc} \boldsymbol{\cdot} \boldsymbol{\mathcal{U}}^Q_c. \end{equation}

Since in the regime of low Reynolds numbers the particle inertia and thus acceleration forces are negligible compared to friction, the hydrodynamic forces on each particle have to equal the negative external forces applied, $\boldsymbol {F}_b = - \boldsymbol {F}^{hydro}_b$, and similarly for torques, $\boldsymbol {T}_b = - \boldsymbol {T}^{hydro}_b$. This implies

(2.6)\begin{equation} \boldsymbol{\mathcal{F}} ={-} \boldsymbol{\mathcal{F}}^{hydro}, \end{equation}

with $\boldsymbol {\mathcal {F}} = (\boldsymbol {F}_1, \boldsymbol {F}_2, \boldsymbol {T}_1, \boldsymbol {T}_2)$. The mobility matrix with components $\boldsymbol {\mu }^{PQ}_{bc}$ is then defined as

(2.7)\begin{equation} \boldsymbol{\mathcal{U}}^P_b =: \sum_{Q, c} \boldsymbol{\mu}^{PQ}_{bc} \boldsymbol{\cdot} \boldsymbol{\mathcal{F}}^Q_c, \end{equation}

and hence $\boldsymbol {\mu } = - {\boldsymbol{\mathsf{R}}}^{-1}$.

For the case of constant fluid viscosity $\eta ^{(0)}$, the mobility matrix is given up to ${O}(s^{-2})$ by (Dhont Reference Dhont1996, pp. 248–256)

(2.8)\begin{equation} \boldsymbol{\mu}^{(0)} = \frac{1}{6 {\rm \pi}a \eta^{(0)}} \begin{pmatrix} {\boldsymbol{\mathsf{I}}} & {\boldsymbol{\mathsf{T}}} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{V}}} \\ {\boldsymbol{\mathsf{T}}} & {\boldsymbol{\mathsf{I}}} & -{\boldsymbol{\mathsf{V}}} & {\boldsymbol{\mathsf{0}}} \\ {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{V}}} & \dfrac{3}{4 a^2} {\boldsymbol{\mathsf{I}}} & {\boldsymbol{\mathsf{0}}} \\ -{\boldsymbol{\mathsf{V}}} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} & \dfrac{3}{4 a^2} {\boldsymbol{\mathsf{I}}} \end{pmatrix}, \end{equation}

with

(2.9)$$\begin{gather} {\boldsymbol{\mathsf{T}}} := \frac{3 a}{4 s} \left( {\boldsymbol{\mathsf{I}}} + \frac{\boldsymbol{s} \otimes \boldsymbol{s}}{s^2} \right), \end{gather}$$

(2.10)$$\begin{gather}{\boldsymbol{\mathsf{V}}} = \frac{3 a}{4 s^3}\,(\boldsymbol{s} \times), \end{gather}$$

and the tensor product denoted by $\otimes$. Furthermore, $\boldsymbol {s} := \boldsymbol {r}_2 - \boldsymbol {r}_1$ is the separation vector between the positions of the two particles, $s = |\boldsymbol{s}|$ and ${\boldsymbol{\mathsf{0}}}$ is the $3 \times 3$ matrix with all entries zero. The tensor $(\boldsymbol {s} \times )$ is defined as $(\boldsymbol {s} \times ) := \epsilon_{ikj} s_k \boldsymbol {e}_i \otimes \boldsymbol {e}_j$, where $\boldsymbol {e}_i$ is the unit vector in the direction associated with the spatial index $i$, $\epsilon_{ikj}$ is the Levi–Civita symbol and we imply summation over repeated spatial indices.

3. Derivation of the correction term to the mobility matrix

We now make use of the Lorentz reciprocal theorem to derive the first-order correction in $\epsilon$ to the mobility matrix of a two-particle system, with particles at positions $\boldsymbol {r}_1$ and $\boldsymbol {r}_2$. Specifically, we consider two systems, where the first has constant viscosity $\eta ^{(0)}$ everywhere in the fluid, and the second is described by the perturbed viscosity profile $\eta (\boldsymbol {r})$ (Oppenheimer et al. Reference Oppenheimer, Navardi and Stone2016). Quantities in the first system will carry a superscript $^{(0)}$ associating them with constant viscosity, whereas quantities corresponding to the second system will carry no superscript in parentheses. In both systems, the particles translate and rotate with velocities $\boldsymbol {\mathcal {U}}^{(0)}$ and $\boldsymbol {\mathcal {U}}$, respectively, inducing velocity, pressure and stress fields $\{\boldsymbol {u}^{(0)} (\boldsymbol {r}), p^{(0)} (\boldsymbol {r}), \boldsymbol {\sigma }^{(0)} (\boldsymbol {r})\}$ and $\{\boldsymbol {u} (\boldsymbol {r}), p (\boldsymbol {r}), \boldsymbol {\sigma } (\boldsymbol {r})\}$, respectively. To shorten notation, we will omit the position-dependence in the following and use explicit notation only for the sake of clarity when necessary. Following the derivation of the standard Lorentz reciprocal theorem (Appendix A, and Oppenheimer et al. Reference Oppenheimer, Navardi and Stone2016), we arrive at the relation

(3.1)\begin{equation} \sum_{P, b} \boldsymbol{\mathcal{F}}^P_b \boldsymbol{\cdot} \boldsymbol{\mathcal{U}}^{(0) P}_b - \sum_{P, b} \boldsymbol{\mathcal{F}}^{(0) P}_b \boldsymbol{\cdot} \boldsymbol{\mathcal{U}}^P_b = 2 \int_V [ \eta(\boldsymbol{r}) - \eta^{(0)} ] {\boldsymbol{\mathsf{E}}} : {\boldsymbol{\mathsf{E}}}^{(0)} \, {\rm d} V. \end{equation}

Here, $V$ denotes the fluid volume and ${\boldsymbol{\mathsf{E}}} : {\boldsymbol{\mathsf{E}}}^{(0)}$ the contraction of the last two indices of ${\boldsymbol{\mathsf{E}}}$ with the first two indices of ${\boldsymbol{\mathsf{E}}}^{(0)}$. Naturally, while for $\eta (\boldsymbol {r}) = \eta ^{(0)}$ one obtains the standard Lorentz reciprocal theorem, we here obtain a relation between the generalised velocities and forces, and a correction term involving the viscosity perturbation and the strain rate in the fluid.

For given particle positions, the fluid velocity at each point is a linear expression in the forces and torques acting on the particles. Therefore, we can express the strain rates ${\boldsymbol{\mathsf{E}}}$ and ${\boldsymbol{\mathsf{E}}}^{(0)}$ as products of the normalised strain rates ${\boldsymbol{\mathsf{B}}}$ and ${\boldsymbol{\mathsf{B}}}^{(0)}$ with respect to the forces and torques acting on the particles, and the generalised forces

(3.2a,b)\begin{equation} {\boldsymbol{\mathsf{E}}} =: \sum_{P, b} {\boldsymbol{\mathsf{B}}}^P_b \boldsymbol{\cdot} \boldsymbol{\mathcal{F}}^P_b, \quad {\boldsymbol{\mathsf{E}}}^{(0)} =: \sum_{P, b} {\boldsymbol{\mathsf{B}}}^{(0) P}_b \boldsymbol{\cdot} \boldsymbol{\mathcal{F}}^{(0) P}_b. \end{equation}

Inserting this expansion together with (2.7) and the analogous equation for the constant viscosity case into (3.1), we obtain

(3.3)\begin{align} &\sum_{P, Q, b, c} ( \boldsymbol{\mathcal{F}}^{P}_b \boldsymbol{\cdot} \boldsymbol{\mu}^{(0) PQ}_{bc} \boldsymbol{\cdot} \boldsymbol{\mathcal{F}}^{(0) Q}_c - \boldsymbol{\mathcal{F}}^{(0) P}_b \boldsymbol{\cdot} \boldsymbol{\mu}^{PQ}_{bc} \boldsymbol{\cdot} \boldsymbol{\mathcal{F}}^{Q}_c ) \nonumber\\ &\quad = 2 \int_V [ \eta(\boldsymbol{r}) - \eta^{(0)} ] \sum_{P, Q, b, c} ({\boldsymbol{\mathsf{B}}}^P_b \boldsymbol{\cdot} \boldsymbol{\mathcal{F}}^P_b ) : ({\boldsymbol{\mathsf{B}}}^{(0) Q}_c \boldsymbol{\cdot} \boldsymbol{\mathcal{F}}^{(0) Q}_c ) \, {\rm d} V. \end{align}

Both sides are linear in $\boldsymbol {\mathcal {F}}$ and $\boldsymbol {\mathcal {F}}^{(0)}$, which are undetermined variables up to here. Hence by suitably renaming the summation indices, using the symmetry of the mobility matrix (Jeffrey & Onishi Reference Jeffrey and Onishi1984), and factoring out $\boldsymbol {\mathcal {F}}$ and $\boldsymbol {\mathcal {F}}^{(0)}$ (Appendix B), we obtain

(3.4)\begin{equation} \boldsymbol{\mu}^{PQ}_{bc} = \boldsymbol{\mu}^{(0) PQ}_{bc} - 2 \int_V [ \eta(\boldsymbol{r}) - \eta^{(0)} ] {\boldsymbol{\mathsf{B}}}^P_b \mathrel{\substack{\bullet \\ \cdot}} {\boldsymbol{\mathsf{B}}}^{(0) Q}_c \, {\rm d} V. \end{equation}

Here, we follow the notation used in Oppenheimer et al. (Reference Oppenheimer, Navardi and Stone2016), defining

(3.5) \begin{equation} \textsf{B}^P_b \mathrel{\substack{\bullet \\ \cdot}} \textsf{B}^{(0) Q}_c := \sum_{k,l} (\textsf{B}^P_b )_{kli}\,( \textsf{B}^{(0) Q}_c )_{klj}\,\boldsymbol{e}_i \otimes \boldsymbol{e}_j. \end{equation}

Realising that $\eta (\boldsymbol {r}) - \eta ^{(0)} = \epsilon \,\eta ^{(1)} (\boldsymbol {r})$ and that we can expand ${\boldsymbol{\mathsf{B}}} = {\boldsymbol{\mathsf{B}}}^{(0)} + {O}(\epsilon ^1)$, we find that the $\epsilon ^1$ correction to the mobility matrix depends on only the normalised strain rate in a constant-viscosity fluid, ${\boldsymbol{\mathsf{B}}}^{(0) P}_b$, instead of ${\boldsymbol{\mathsf{B}}}^P_b$. Consequently, for the $\epsilon ^1$ correction to the mobility matrix, we obtain

(3.6)\begin{equation} \boldsymbol{\mu}^{(1) PQ}_{bc} ={-} 2 \int_V \eta^{(1)} (\boldsymbol{r})\,{\boldsymbol{\mathsf{B}}}^{(0) P}_b \mathrel{\substack{\bullet \\ \cdot}} {\boldsymbol{\mathsf{B}}}^{(0) Q}_c \, {\rm d} V, \end{equation}

which holds irrespective of the actual form of the gradient and the number of particles in the system. Note that in Oppenheimer et al. (Reference Oppenheimer, Navardi and Stone2016), the strain rates have been expanded with respect to the generalised particle velocities instead of the generalised forces. This yields a similar result for the first-order correction to the resistance matrix rather than for the mobility matrix.

While the normalised strain rate for a single particle in constant viscosity can be obtained directly from the flow field associated with translation and rotation of a single particle, in the case of two or more particles, one has to also account for the disturbance flows that other particles induce in response to the flow field produced by a first particle. Depending on whether the expansion of the strain rates is done with respect to the generalised forces or velocities, these other particles must be assumed to be force- and torque-free, or to have zero (angular) velocity in the calculation of the normalised strain rates, respectively. By using a reflection method, one can account for all those disturbance flows ordered by ascending powers of the inverse particle separation (Kim & Karilla Reference Kim and Karilla2005, Chapter 8).

At this point, expanding the strain rates with respect to the generalised forces simplifies the subsequent calculations drastically in comparison to an expansion with respect to the generalised velocities. The reason for this is that the disturbance flow field induced by a force- and torque-free particle scales to leading order with the local gradient of the flow field in which it is immersed, namely the local strain rate of the flow, whereas a particle with zero velocity induces, in general, a disturbance flow proportional to the local flow field (Kim & Karilla Reference Kim and Karilla2005, p. 189). This is also clear intuitively: a force-free particle translates and rotates according to the surrounding flow field, whereas a particle with prescribed zero velocity must resist any external flow, inducing a stronger disturbance flow.

A spherical particle located at the origin and subject to a force in a constant viscosity fluid produces a flow field given by ${\boldsymbol{\mathsf{T}}}^{trans} (\boldsymbol {r}) /(6 {\rm \pi}\eta ^{(0)} a)$ contracted with the external force, with (Dhont Reference Dhont1996, p. 246) the Rotne–Prager tensor

(3.7)\begin{equation} {\boldsymbol{\mathsf{T}}}^{trans} (\boldsymbol{r}) := \left( \frac{3a}{4r} \left({\boldsymbol{\mathsf{I}}} + \frac{\boldsymbol{r} \otimes \boldsymbol{r}}{r^2} \right) + \frac{a^3}{4r^3} \left( {\boldsymbol{\mathsf{I}}} - 3\,\frac{\boldsymbol{r} \otimes \boldsymbol{r}}{r^2} \right)\right), \end{equation}

and $r = |\boldsymbol {r}|$. Similarly, the flow field induced by a particle subject to a torque is given by ${\boldsymbol{\mathsf{T}}}^{rot} (\boldsymbol {r})/(8 {\rm \pi}\eta ^{(0)} a^3)$ contracted with the applied torque, with (Dhont Reference Dhont1996, p. 248)

(3.8)\begin{equation} {\boldsymbol{\mathsf{T}}}^{rot} (\boldsymbol{r}) :={-} \left(\frac{a}{r} \right)^3 (\boldsymbol{r} \times), \end{equation}

being the rotlet. We then define the normalised single particle strain rates for translation and rotation as the symmetrised gradient of the corresponding flow field, where we use explicit index notation for the sake of clarity:

(3.9)$$\begin{gather} (\textsf{P}^{(0) t})_{kli} (\boldsymbol{r}) := \frac{1}{6 {\rm \pi}\eta^{(0)} a}\,\frac{1}{2}\,[\partial_k (\textsf{T}^{trans})_{li} (\boldsymbol{r}) + \partial_l (\textsf{T}^{trans})_{ki} (\boldsymbol{r})], \end{gather}$$

(3.10)$$\begin{gather}(\textsf{P}^{(0) r})_{kli} (\boldsymbol{r}) := \frac{1}{8 {\rm \pi}\eta^{(0)} a^3}\,\frac{1}{2}\,[\partial_k (\textsf{T}^{rot})_{li} (\boldsymbol{r}) + \partial_l (\textsf{T}^{rot})_{ki} (\boldsymbol{r})]. \end{gather}$$

The normalised strain rates in the two-particle system, ${\boldsymbol{\mathsf{B}}}^{(0) P}_b$, can then be expanded by means of the reflection method (Appendix C) in terms of the single-particle strain rates as

(3.11)\begin{equation} {\boldsymbol{\mathsf{B}}}^{(0) t}_b \equiv {\boldsymbol{\mathsf{B}}}^{(0) t}_b (\boldsymbol{r}) = {\boldsymbol{\mathsf{P}}}^{(0) t} (\boldsymbol{r} - \boldsymbol{r}_b) + {\boldsymbol{\mathsf{P}}}^{(0) s} (\boldsymbol{r} - \boldsymbol{r}_c) : {\boldsymbol{\mathsf{P}}}^{(0) t} (\boldsymbol{r}_c - \boldsymbol{r}_b) + {O}(s^{{-}3}),\end{equation}

and

(3.12)\begin{equation} {\boldsymbol{\mathsf{B}}}^{(0) r}_b \equiv {\boldsymbol{\mathsf{B}}}^{(0) r}_b (\boldsymbol{r}) = {\boldsymbol{\mathsf{P}}}^{(0) r} (\boldsymbol{r} - \boldsymbol{r}_b) + {O}(s^{{-}3}). \end{equation}

Here, ${\boldsymbol{\mathsf{P}}}^{(0) s}$ denotes the normalised strain rate associated with the disturbance flow field that a spherical particle immersed in a local strain rate induces (see Appendix C for details).

We have worked out the reflection scheme up to ${O}(s^{-2})$, which will suffice to calculate the correction to the mobility matrix up to order ${O}(s^{-1})$ in a linear viscosity gradient and up to order ${O}(s^{-2})$ for particles with a temperature contrast to the fluid. Inserting the expansions (3.11) and (3.12) into (3.6), we see that the task of calculating the integrals (3.6) reduces to calculating integrals of the form

(3.13)\begin{equation} {\boldsymbol{\mathsf{Q}}}^{(1) PQ}_{bc} :={-} 2 \int_V \eta^{(1)}(\boldsymbol{r})\,{\boldsymbol{\mathsf{P}}}^{(0) P} (\boldsymbol{r} - \boldsymbol{r}_b) \mathrel{\substack{\bullet \\ \cdot}} {\boldsymbol{\mathsf{P}}}^{(0) Q} (\boldsymbol{r} - \boldsymbol{r}_c) \, {\rm d} V, \end{equation}

for $P, Q = t, r, s$, as discussed in the following sections.

4. Particles in a large linear viscosity gradient

We first consider the mobility matrix of two particles within a large linear viscosity gradient, i.e. with an extent $d$ larger by at least an order of magnitude than the distance $s$ between the two particles, $d \gg s$ (figure 1). For simplicity, we assume here that the particles do not perturb the viscosity gradient. In particular, we assume that the viscosity profile is given by

(4.1)\begin{equation} \eta (\boldsymbol{r}) = \eta_0 + \epsilon \boldsymbol{H}^{(1)} \boldsymbol{\cdot} \boldsymbol{r} \end{equation}

for all $|\boldsymbol {r}| < d$, where $\boldsymbol {H}^{(1)}$ is the vectorial gradient, and $\eta _0$ is a constant viscosity background. We introduce a second dimensionless parameter $\kappa := a/d$ as the ratio of bead radius and size of the viscosity gradient. Assuming that the particle distance is sufficiently large as compared to the bead radius $a \ll s$, we restrict our calculation of the correction to the mobility matrix to the leading orders in $a/s$, and include here all terms up to ${O}(s^{-1})$. This implies that for the remainder of the section on linear gradients, we impose $d \gg s \gg a$. However, since in constant viscosity the leading-order approximations for the mobility terms apply with errors below 5 % if the beads are separated by only 5 bead radii, it is expected that the leading-order results for a viscosity gradient are of similar accuracy at such distances. While these seem to be quite some restrictions, we emphasise that viscosity gradients often occur on length scales of millimetres, in liquids filled with micron-sized colloids or active swimmers (Stehnach et al. Reference Stehnach, Waisbord, Walkama and Guasto2021).

Figure 1. Sketch of the prototypical interface-like viscosity gradient. (a) Schematic of the position-dependent viscosity. (b) Sketch of the linear fluid viscosity gradient. Darker blue denotes regions of higher viscosity, and brighter blue denotes areas of lower viscosity. Particles can be placed close to the symmetry plane of the gradient, as shown in green and discussed in § 4.1, or can be placed asymmetrically, shown in yellow and discussed in § 4.2. Enlarged sections are used to illustrate the parameters $a$ and $\boldsymbol {s}$.

A prototypical interface-like linear viscosity gradient is, without restriction of generality along the $y$-direction, given by

(4.2)\begin{equation} \eta (\boldsymbol{r}) = \eta_0 + \begin{cases} - \epsilon \eta_0, & y <{-}d, \\[3pt] \left(0, \epsilon\,\dfrac{\eta_0}{d}, 0\right) \boldsymbol{\cdot} \boldsymbol{r}, & -d \leq y \leq d,\\ \epsilon \eta_0, & d < y, \end{cases}\end{equation}

where $\boldsymbol {r} = (x, y, z)$. It describes a viscosity perturbation which is linear in an infinite slab of diameter $2d$ and is constant in the two-half spaces separated by the slab with a continuous transition (figure 1), avoiding negative viscosity. Such viscosity gradients can be expected to be found e.g. at the interface of two immiscible fluids of different viscosities.

For linear viscosity gradients, we will generally choose the absolute reference viscosity $\eta ^{(0)}$, which defines the $\epsilon ^1$ viscosity perturbation as $\eta ^{(1)} (\boldsymbol {r}) := (\eta (\boldsymbol {r}) - \eta ^{(0)})/\epsilon$, to be the viscosity at some point close to the interacting particles. In contrast, $\eta _0$ defines the viscosity perturbation via (4.2), and $\epsilon \eta _0$ characterises the maximum viscosity contrast between the actual viscosity profile $\eta (\boldsymbol {r})$ and the background viscosity. While we can assume $\eta ^{(0)} = \eta _0$ for the case of both particles being placed close to the plane of symmetry in the interface-like gradient (§ 4.1), the distinction between $\eta ^{(0)}$ and $\eta _0$ becomes relevant when both particles are placed away from the plane of symmetry. The corrections arising in this more general case due to the broken symmetry are investigated in § 4.2.

4.1. Symmetrically placed particles in an odd symmetric viscosity gradient

In this subsection, we assume that the viscosity perturbation is odd symmetric with respect to the origin, i.e. $\eta ^{(1)} (- \boldsymbol {r}) = - \eta ^{(1)} (\boldsymbol {r})$, and both particles are close to the centre of symmetry, which we assume to be the origin of the coordinate system, i.e. $|\boldsymbol {r}_i| \ll d$. For the interface-like viscosity gradient, which by definition satisfies the symmetry assumption, this is the case if the particles are positioned close to the plane of symmetry characterised by $y = 0$. Thus in this case we can assume for the reference viscosity $\eta ^{(0)} = \eta _0$. It can be shown that an infinite linear viscosity gradient, although associated with negative viscosity in parts of the fluid, yields mathematically the same result for the correction to the mobility matrix (§ D.1 of Appendix D) up to order $\epsilon ^1 \kappa ^1$ as the odd symmetric gradient. This allows us to compare our results with previous works considering infinite gradients (Datt & Elfring Reference Datt and Elfring2019; Laumann & Zimmermann Reference Laumann and Zimmermann2019).

We defer the explicit calculation of the ${\boldsymbol{\mathsf{Q}}}^{(1)}$ terms, which can be done with the help of a computer algebra system (Wolfram Research Inc. 2020), to Appendix D. The results are presented in (D14) and (D18). Nonetheless, we use these results to calculate the mobility matrix $\boldsymbol {\mu }^{(1) PQ}_{bc}$ by inserting (3.11) and (3.12) into (3.6), and obtain, after applying the reflection method,

(4.3) \begin{align} \left. \begin{array}{c} \displaystyle \boldsymbol{\mu}^{(1) tt}_{bc} = {\boldsymbol{\mathsf{Q}}}^{(1)tt}_{bc} + \sum_{e \neq c} {\boldsymbol{\mathsf{Q}}}^{(1)ts}_{be} : {\boldsymbol{\mathsf{P}}}^{(0) t} (\boldsymbol{r}_e - \boldsymbol{r}_c) + \sum_{f \neq b} {\boldsymbol{\mathsf{P}}}^{(0) t} (\boldsymbol{r}_f - \boldsymbol{r}_b) \mathrel{\substack{\bullet \\ \cdot}} {\boldsymbol{\mathsf{Q}}}^{(1)st}_{fc} + {O}(s^{{-}3}), \\ \displaystyle \boldsymbol{\mu}^{(1) tr}_{bc} = {\boldsymbol{\mathsf{Q}}}^{(1)tr}_{bc} + \sum_{f \neq b} {\boldsymbol{\mathsf{P}}}^{(0) t} (\boldsymbol{r}_f - \boldsymbol{r}_b) \mathrel{\substack{\bullet \\ \cdot}} {\boldsymbol{\mathsf{Q}}}^{(1)sr}_{fc} + {O}(s^{{-}3}), \\ \displaystyle \boldsymbol{\mu}^{(1) rr}_{bc} = {\boldsymbol{\mathsf{Q}}}^{(1)rr}_{bc} + {O}(s^{{-}3}), \end{array}\right\} \end{align}

where $b, c, e$ and $f$ are independent particle indices. We use the sum notation in (4.3) to cover both cases, $b = c$ and $b \neq c$. It can, however, be shown (§ D.4 of Appendix D) that the equations in (4.3) reduce to

(4.4)\begin{equation} \boldsymbol{\mu}^{(1) PQ}_{bc} = {\boldsymbol{\mathsf{Q}}}^{(1) PQ}_{bc} + {O}(s^{{-}2}) \end{equation}

for the linear viscosity gradient, i.e. all other summands in (4.3) scale at most as ${O}(s^{-2})$ and can be neglected.

Including all terms up to ${O}(s^{-1})$, the $\epsilon ^1 \kappa ^1$ correction to the mobility matrix then is given by

(4.5) \begin{align} \boldsymbol{\mu}^{(1)} &= \tfrac{1}{6 {\rm \pi}a (\eta^{(0)})^2} \nonumber\\ &\quad \times \scriptsize{\begin{pmatrix} - \eta^{(1)} (\boldsymbol{r}_1)\,{\boldsymbol{\mathsf{I}}} & - \eta^{(1)} \left( \dfrac{\boldsymbol{r}_1 + \boldsymbol{r}_2}{2} \right) {\boldsymbol{\mathsf{T}}} + {\boldsymbol{\mathsf{W}}} & \dfrac{1}{4} (\boldsymbol{\nabla} \eta^{(1)} \times) & - \eta^{(1)} \left( \dfrac{\boldsymbol{r}_1 + \boldsymbol{r}_2}{2} \right) {\boldsymbol{\mathsf{V}}} + {\boldsymbol{\mathsf{Y}}} \\[6pt] - \eta^{(1)} \left( \dfrac{\boldsymbol{r}_1 + \boldsymbol{r}_2}{2} \right) {\boldsymbol{\mathsf{T}}} - {\boldsymbol{\mathsf{W}}} & -\eta^{(1)} (\boldsymbol{r}_2)\,{\boldsymbol{\mathsf{I}}} & \eta^{(1)} \left( \dfrac{\boldsymbol{r}_1 + \boldsymbol{r}_2}{2} \right) {\boldsymbol{\mathsf{V}}} + {\boldsymbol{\mathsf{Y}}} & \dfrac{1}{4} (\boldsymbol{\nabla} \eta^{(1)} \times) \\[6pt] -\dfrac{1}{4} (\boldsymbol{\nabla} \eta^{(1)} \times) & - \eta^{(1)} \left( \dfrac{\boldsymbol{r}_1 + \boldsymbol{r}_2}{2} \right) {\boldsymbol{\mathsf{V}}} + {\boldsymbol{\mathsf{Y}}}^{\rm T} & -\dfrac{3}{4 a^2}\,\eta^{(1)} (\boldsymbol{r}_1)\,{\boldsymbol{\mathsf{I}}} & {\boldsymbol{\mathsf{0}}} \\[6pt] \eta^{(1)} \left( \dfrac{\boldsymbol{r}_1 + \boldsymbol{r}_2}{2} \right) {\boldsymbol{\mathsf{V}}} + {\boldsymbol{\mathsf{Y}}}^{\rm T} & -\dfrac{1}{4} (\boldsymbol{\nabla} \eta^{(1)} \times) & {\boldsymbol{\mathsf{0}}} & -\dfrac{3}{4 a^2}\,\eta^{(1)} (\boldsymbol{r}_2)\, {\boldsymbol{\mathsf{I}}} \end{pmatrix}}, \end{align}

with

(4.6)\begin{equation} {\boldsymbol{\mathsf{W}}} := \frac{3 a}{8 s}\,( \boldsymbol{\nabla} \eta^{(1)} \otimes \boldsymbol{s} - \boldsymbol{s} \otimes \boldsymbol{\nabla} \eta^{(1)} ) \end{equation}

and

(4.7)\begin{equation} {\boldsymbol{\mathsf{Y}}} := \frac{3 a}{8 s^3}\,\boldsymbol{s} \otimes (\boldsymbol{s} \times \boldsymbol{\nabla} \eta^{(1)}). \end{equation}

The corrections to the translational and rotational self-mobilities of both particles, i.e. the diagonal elements of (4.5), are proportional to the negative local value of the viscosity perturbation, associating higher viscosity with lower bead self-mobilities as expected. In fact, the corrections to the self-mobilities correspond to the first-order term in a Taylor expansion of the constant-viscosity expressions with the local viscosity $\eta (\boldsymbol {r}_b) = \eta ^{(0)} + \epsilon \,\eta ^{(1)} (\boldsymbol {r}_b)$ inserted, around $\epsilon = 0$.

In contrast to the case of constant viscosity, the viscosity gradient couples the translational and rotational modes of a single particle, i.e. a particle experiences rotation when it is centrally subject to a force. Similarly, translation will appear if the particle is subject to a torque. Assuming a force $\boldsymbol {F}$ acting on a particle, the resulting angular velocity is given from (4.5) as

(4.8)\begin{equation} \boldsymbol{\varOmega} ={-}\frac{1}{24 {\rm \pi}a (\eta^{(0)})^2}\,\boldsymbol{\nabla} \eta \times \boldsymbol{F}. \end{equation}

This can be understood as an effect of the imbalance of the friction forces (Datt & Elfring Reference Datt and Elfring2019). On the side of higher viscosity, the particle experiences higher drag, while on the side of lower viscosity, the effect is inverse. When a force perpendicular to the viscosity gradient acts on a particle, it will therefore induce a rotation as if the particle adhered more strongly on the side of higher viscosity. Similarly, a particle subject to a torque moves translationally when the torque has a component orthogonal to the viscosity gradient. Our findings here reproduce exactly the results presented in Datt & Elfring (Reference Datt and Elfring2019) for infinite viscosity gradients.

In a more intricate fashion, the viscosity gradient alters the interaction between two different particles. From the $\boldsymbol {\mu }^{(1) tt}_{21}$ term in (4.5), we see that the $\epsilon ^1 \kappa ^1$ component of the interaction between two particles is composed of two terms. The first is the $\eta ^{(1)} ((\boldsymbol {r}_1 + \boldsymbol {r}_2)/2) {\boldsymbol{\mathsf{T}}}$ term, which can be understood as an expansion of the constant viscosity interaction, similarly to the $\boldsymbol {\mu }^{(1) tt}_{bb}$ and $\boldsymbol {\mu }^{(1) rr}_{bb}$ terms. Here, the viscosity is evaluated at the geometric centre between the particles. The second term in the $\epsilon ^1 \kappa ^1$ component of the interaction is associated with ${\boldsymbol{\mathsf{W}}}$. It gives rise to novel effects, in particular because it introduces a velocity component orthogonal to the force, when the force is parallel to the connection of both particles. This component is non-zero as soon as the viscosity gradient has a component orthogonal to the connection of both particles.

For the purposes of illustration, we derive the $\epsilon ^1 \kappa ^1$ correction to the flow field produced by a particle subject to some force $\boldsymbol {F}_1$ within the viscosity gradient from the $\boldsymbol {\mu }^{(1) tt}_{21}$ component in (4.5), by treating the second sphere as a test particle. We assume that the particle producing the flow is positioned at the origin of the coordinate system, $\boldsymbol {r}_1 = 0$, where we also assume the viscosity perturbation to vanish, i.e. $\eta ^{(1)} (\boldsymbol {r}_1) = 0$. The first-order correction to the flow field is then given as the velocity of the second particle located at $\boldsymbol {r}_2 = \boldsymbol {r}$,

(4.9)\begin{align} &\boldsymbol{u}^{(1)}_F (\boldsymbol{r}) = \boldsymbol{U}_2^{(1)} \nonumber\\ &\quad = \frac{1}{6 {\rm \pi}a (\eta^{(0)})^2}\,\frac{3a}{8|\boldsymbol{r}|}\left[ - (\boldsymbol{r} \boldsymbol{\cdot} \boldsymbol{\nabla} \eta^{(1)}) \left( {\boldsymbol{\mathsf{I}}} + \frac{\boldsymbol{r} \otimes \boldsymbol{r}}{|\boldsymbol{r}|^2} \right) - ( \boldsymbol{\nabla} \eta^{(1)} \otimes \boldsymbol{r} - \boldsymbol{r} \otimes \boldsymbol{\nabla} \eta^{(1)} ) \right] \boldsymbol{\cdot} \boldsymbol{F}_1, \end{align}

where the first contribution comes from the $- \eta ^{(1)} ((\boldsymbol {r}_1 + \boldsymbol {r}_2)/2) {\boldsymbol{\mathsf{T}}}$ term, and the second contribution comes from the ${\boldsymbol{\mathsf{W}}}$ term (figure 2). The colour in figure 2 indicates the ratio between the $\epsilon ^1 \kappa ^1$ correction flow and the velocity $U^{(0)} = F_1/(6 {\rm \pi}\eta ^{(0)} a)$ at which the particle would move in constant viscosity $\eta ^{(0)}$ when subject to the force $\boldsymbol {F}_1$. This correction flow field is divergence-free and hence satisfies the incompressibility condition as expected. It is in agreement with the generalised Oseen tensor that was derived in Laumann & Zimmermann (Reference Laumann and Zimmermann2019).

Figure 2. The $\epsilon ^1 \kappa ^1$ contribution to the far field flow produced by a sphere subject to a force within the viscosity gradient for different directions of the force. The colour indicates the ratio between the $\epsilon ^1 \kappa ^1$ correction flow and the velocity $U^{(0)} = F_1/(6 {\rm \pi}\eta ^{(0)} a)$ at which the particle would move in constant viscosity $\eta ^{(0)}$ when subject to the force $\boldsymbol {F}_1$.

When the force $\boldsymbol {F}_1$ is applied parallel to the viscosity gradient, it causes an $\epsilon ^1 \kappa ^1$ flow directed inwards along the direction of the viscosity gradient, and a flow directed outwards in the direction orthogonal to the gradient (figure 2a). When the force is acting in the opposite direction, i.e. down the viscosity gradient, also the correction flow field inverts (figure 2b), as a direct consequence of the linearity of the problem. Qualitatively, this effect can be understood as follows: in the region above the particle, the fluid has higher viscosity and hence will react more weakly to the force than a constant-viscosity fluid. Therefore, the correction flow field is oriented oppositely to the force here. In the region below the particle with decreased viscosity, the fluid reacts more strongly and the correction flow is therefore parallel to the force.

When the force is acting perpendicular to the viscosity gradient (figure 2c), the correction flow again follows the direction of the force in the region of lower viscosity, and opposes it in the region of higher viscosity. On the axis perpendicular to the viscosity gradient, a correction flow perpendicular to the force applied is found, such that again a pattern of inwards flow and outwards flow is observed. When the force on the particle at the origin acts in some arbitrary direction (figure 2d), the resulting flow field will be a linear superposition of the flow fields in figures 2(a,c).

As an important consequence of the correction to the mobility matrix (4.5), the far field correction flow at position $\boldsymbol {r}$ depends only on the direction of $\boldsymbol {r}$ but is constant with respect to the absolute value of the distance, i.e. it does not decay radially. This contrasts with the constant-viscosity Oseen tensor, which decays with the first power of the inverse distance from the particle. However, this correction flow field is valid only as long as the particle separation is, by order of magnitude, smaller than $d$. Therefore, the correction flow, being proportional to $\epsilon \kappa$, always compares small to the Stokeslet flow field, which scales as $a/s$. Consequently, the ratio of the correction to the constant-viscosity term at distances of ${O}(s)$ scales as $\epsilon \kappa /(a/s) = \epsilon s/d$. Inserting values from the experimental system in Stehnach et al. (Reference Stehnach, Waisbord, Walkama and Guasto2021), we obtain $\epsilon \approx 0.5$ for the strongest viscosity gradient considered, which is of total size $2d = 1$ mm. Assuming that the interacting objects (particles or Chlamydomonas cells of size of $\sim 10 \mathrm {\mu }$m) have a distance of $100\ \mathrm {\mu }$m, we find that the first-order correction amounts to the order of 10 % of the constant-viscosity term.

The $\boldsymbol {\mu }^{(1) rt}_{21}$ term in (4.5) corresponds to the angular velocity $\boldsymbol {\varOmega }^{(1)}_F$ of particle 2 when particle 1 is subject to a force. With $\boldsymbol {r}_1 = \boldsymbol {0}$ and $\boldsymbol {r}_2 = \boldsymbol {r}$, we find

(4.10)\begin{equation} \boldsymbol{\varOmega}^{(1)}_F(\boldsymbol{r}) = \boldsymbol{\varOmega}_2^{(1)} = \frac{1}{6 {\rm \pi}a (\eta^{(0)})^2}\,\frac{3a}{8|\boldsymbol{r}|^3}[(\boldsymbol{r} \boldsymbol{\cdot} \boldsymbol{\nabla} \eta^{(1)}) (\boldsymbol{r} \times) + (\boldsymbol{r} \times \boldsymbol{\nabla} \eta^{(1)}) \otimes \boldsymbol{r}] \boldsymbol{\cdot} \boldsymbol{F}_1, \end{equation}

where the first summand results from $\eta ^{(1)} ((\boldsymbol {r}_1 + \boldsymbol {r}_2)/2) {\boldsymbol{\mathsf{V}}}$ and the second summand results from ${\boldsymbol{\mathsf{Y}}}^\textrm {T}$.

From the $\boldsymbol {\mu }^{(1) tr}_{21}$ term in the mobility matrix (4.5), we extract the flow field that a particle subject to a torque induces in the viscosity gradient. Assuming that, similarly to before, the viscosity perturbation vanishes at the position of the first particle, the velocity of the second particle acting as a test particle is then given as

(4.11)\begin{align} \boldsymbol{u}^{(1)}_T (\boldsymbol{r}) = \boldsymbol{U}_2^{(1)} = \frac{1}{6 {\rm \pi}a (\eta^{(0)})^2}\,\frac{3a}{8|\boldsymbol{r}|^3} [(\boldsymbol{\nabla} \eta^{(1)} \boldsymbol{\cdot} \boldsymbol{r}) (\boldsymbol{r} \times \boldsymbol{T}_1) + \boldsymbol{r} ( (\boldsymbol{r} \times \boldsymbol{\nabla} \eta^{(1)}) \boldsymbol{\cdot} \boldsymbol{T}_1)]. \end{align}

If a torque acts orthogonal to the viscosity gradient (figure 3a), then a correction flow predominantly in the direction of $\boldsymbol {\nabla } \eta \times \boldsymbol {T}_1$ is induced. This result holds strictly at $z = 0$, while for $z \neq 0$ an additional component of the flow in the $z$-direction is found. This finding is intuitive since the particle itself moves in the same direction as a result of the correction to its self-mobility. For a torque parallel to the viscosity gradient (figure 3b), an additional rotating flow is observed that has the direction of rotation prescribed by $\boldsymbol {T}_1$ in the region of lower viscosity (here $z < 0$) and opposite direction in the region of increased viscosity (here $z > 0$). In the plane $z = 0$, the correction flow vanishes. Also, this is intuitive, because the fluid opposes the rotation more strongly in the region of higher viscosity, thus resulting in a correction flow with opposite direction of rotation, and the opposite effect in the lower-viscosity region.

Figure 3. Far field $\epsilon ^1 \kappa ^1$ correction flow induced by a particle subject to an externally applied torque along the direction indicated by the blue arrow, with the direction of the viscosity gradient (a) orthogonal to it, and (b) parallel to it. The colours of the arrows indicate the ratio of the $\epsilon ^1 \kappa ^1$ correction flow and $\varOmega ^{(0)} a$, where $\varOmega ^{(0)} = 1/(8 {\rm \pi}\eta ^{(0)} a^3) |\boldsymbol {T}_1|$ is the angular velocity of the particle in constant viscosity $\eta ^{(0)}$ when subject to the torque $\boldsymbol {T}_1$.

In both cases, the correction flow decays radially as $r^{-1}$, i.e. by one order slower than in the constant-viscosity case where the rotlet decays as $r^{-2}$. Similarly to the flow field induced by a particle subject to a force, the correction flow field resulting from the linear viscosity gradient decays by one order slower than the constant-viscosity field does, and the ratio of the correction flow to the flow at constant viscosity scales as $\epsilon s/d$. A linear viscosity gradient thus induces long-range interactions between the particles, as long as the particle separation is much shorter than the size of the gradient.

The results presented here for the hydrodynamic interaction of particles in a linear viscosity gradient are also valid for more than two particles (Appendix F) as three-body effects decay quicker than $1/s$, thus particle interaction is pairwise at the level of approximation chosen here.

4.2. Asymmetric particle placement in the viscosity gradient

In the previous subsection, we assumed that the viscosity perturbation is odd symmetric and both particles are located close to the centre of symmetry or the central plane of the interface-like gradient when compared to the gradient size $d$. Since in a realistic setting the particle positions in the viscosity gradient are arbitrary, here we consider the interacting particles at varying positions within a finite-size gradient, and calculate the effect of the viscosity gradient compared to the interaction in constant viscosity with the local viscosity as reference. But we require that first, the distance between both particles is still smaller than the gradient size, $s \ll d$, and second, the distance of each particle to the boundary of the linear gradient in any spatial direction is of the same order of magnitude as $d$, i.e. we assume that both particles do not come close to the edge of the linear gradient at length scales of $d$.

Considering the scaling properties of the integrand in (3.13), it can be shown that only the term ${\boldsymbol{\mathsf{Q}}}^{(1) tt}_{bc} = \boldsymbol {\mu }^{(1) tt}_{bc} + {O}(s^{-2})$ is affected at order $\epsilon ^1 \kappa ^1$ if the particles are placed asymmetrically within the gradient (yellow particles in figure 1). The other ${\boldsymbol{\mathsf{Q}}}^{(1)}$ terms involving rotation remain unchanged (see § D.1 of Appendix D for details). Since the integral (3.13) is linear in the viscosity perturbation, the correction is now given by

(4.12)\begin{equation} {\boldsymbol{\mathsf{Q}}}^{(1) tt}_{asym} ={-} \frac{2}{(6 {\rm \pi}a \eta^{(0)})^2} \int_V [ \eta^{(1)} (\boldsymbol{r}) - \tilde{\eta}^{(1)} (\boldsymbol{r}) ]\,\frac{27 a^2}{8 |\boldsymbol{r} - \boldsymbol{r}_0|^6}\,(\boldsymbol{r} - \boldsymbol{r}_0) \otimes (\boldsymbol{r} - \boldsymbol{r}_0) \, {\rm d} V. \end{equation}

Here, $\boldsymbol{r}_0$ is a reference point close to the positions of both particles, i.e. $|\boldsymbol{r}_b - \boldsymbol{r}_0| \ll d$ for b = 1, 2. $\tilde{\eta}^{(1)} = \boldsymbol{H}^{(1)} \boldsymbol{\cdot} (\boldsymbol{r} - \boldsymbol{r}_0)$ denotes the linear viscosity gradient extended to infinity, and $\eta^{(0)} =\eta(\boldsymbol{r}_0)$ is the corresponding reference viscosity, such that the first-order viscosity perturbation vanishes at $\boldsymbol{r}_0$. As a consistency check, it can be verified easily that (4.12) vanishes due to symmetry if $\eta ^{(1)} (\boldsymbol {r})$ is odd symmetric about $\boldsymbol {r}_0$. We point out that this correction applies to both the translational self-mobilities ($\boldsymbol {\mu }^{(1) tt}_{bb}$) as well as the translational interaction terms in the mobility matrix ($\boldsymbol {\mu }^{(1) tt}_{bc}$ with $b \neq c$). While here we present the correction for particles placed asymmetrically in the interface-like gradient, a similar correction arises if symmetry is broken by the shape of linear region of the viscosity gradient.

The full first-order mobility matrix in the asymmetric case is then given by

(4.13)\begin{equation} \boldsymbol{\mu}^{(1)}_{tot} = \boldsymbol{\mu}^{(1)} + \boldsymbol{\mu}^{(1)}_{asym} = \boldsymbol{\mu}^{(1)} + \begin{pmatrix} {\boldsymbol{\mathsf{Q}}}^{(1) tt}_{asym} & {\boldsymbol{\mathsf{Q}}}^{(1) tt}_{asym} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} \\[3pt] {\boldsymbol{\mathsf{Q}}}^{(1) tt}_{asym} & {\boldsymbol{\mathsf{Q}}}^{(1) tt}_{asym} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} \\[3pt] {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} \\[3pt] {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} & {\boldsymbol{\mathsf{0}}} \\ \end{pmatrix}, \end{equation}

where $\boldsymbol {\mu }^{(1)}$ is the result from (4.5). The correction due to asymmetry yields an additional velocity term for both particles 1 and 2, if the sum of the forces on both particles does not vanish:

(4.14)\begin{equation} \boldsymbol{U}_1^{asym} = \boldsymbol{U}_2^{asym} = {\boldsymbol{\mathsf{Q}}}^{(1) tt}_{asym} \boldsymbol{\cdot} (\boldsymbol{F}_1 + \boldsymbol{F}_2). \end{equation}

Notably, it is the sum of both forces that induces the additional velocity, such that a force-free particle 1 experiences the same additional velocity from the asymmetry as particle 2, if particle 2 is subject to a force.

To illustrate this effect, we consider the particles to be located near a point with $y_0 = y'_0 d$ within the interface-like viscosity gradient, where $y'_0$ is independent of $d$, and $\eta ^{(0)}$ is the reference viscosity at $y_0$. The correction term ${\boldsymbol{\mathsf{Q}}}^{(1) tt}_{asym}$ in dependence of $y_0$ is then calculated in § D.5 of Appendix D

(4.15)\begin{equation} {\boldsymbol{\mathsf{Q}}}^{(1) tt}_{asym} = \frac{1}{6 {\rm \pi}a (\eta^{(0)})^2}\,\frac{9 a |\boldsymbol{\nabla} \eta^{(1)}|}{16} {\rm arctanh}(y_0/d)\,( {\boldsymbol{\mathsf{I}}} + \boldsymbol{e}_y \otimes \boldsymbol{e}_y).\end{equation}

The dependence of this correction term on the relative position in the viscosity gradient is illustrated in figure 4, with the arrows indicating the direction and the colour indicating the strength of $\boldsymbol {U}_1^{asym} = \boldsymbol {U}_2^{asym}$.

Figure 4. Additional velocity $\boldsymbol {U} = \boldsymbol {U}_1^{asym}$ due to asymmetric particle placement in dependence on the position of both particles in the interface-like viscosity gradient. Since the particles have relative distance $s \ll d$, they can be considered as located at the same position at scales of $d$. (a) The correction from asymmetry for a force $\boldsymbol {F} = \boldsymbol {F}_1 + \boldsymbol {F}_2$ applied along the $y$-axis. (b) The correction for a force applied along the $x$-axis. The colour indicates the ratio between the $\epsilon ^1 \kappa ^1$ correction flow and the velocity $U^{(0)} = F/(6 {\rm \pi}\eta ^{(0)} a)$ at which a particle would move in constant viscosity $\eta ^{(0)}$ when subject to the force $\boldsymbol {F}$.

In the half-slab with higher viscosity ($y_0 > 0$ in figure 4), $\boldsymbol {U}_1^{asym}$ is parallel to the force, while in the other half of the slab it is directed oppositely to the force. Consistently with the results for odd symmetric gradients, $\boldsymbol {U}_1^{asym}$ vanishes in the central plane of the gradient. Intuitively, this result can be understood as follows: particles placed in the half-slab of higher viscosity have smaller distance to the edge of the linear gradient at higher viscosity than to the low-viscosity edge. Then, relative to the particle position, the fraction of the linear gradient with higher viscosity is smaller than the fraction with lower viscosity. In comparison to an odd symmetric gradient, where those two fractions are the same, a correction along the force is expected and observed, since the odd symmetric gradient overestimates the viscosity in parts of the fluid. This is also clear from (4.12): for $y_0 > 0$, the actual viscosity is smaller than that associated with an odd symmetric viscosity profile with respect to the particles’ position, thus the term in square brackets is negative for parts of the fluid and zero otherwise. The remaining integrand is in general positive, thus the overall result of (4.12) is also positive, corresponding qualitatively to a correction along the direction of the applied force. For particles placed in the lower viscosity half-slab, the effect inverses.

For $|y_0| \lesssim d/2$, we have $\textrm {arctanh}(y_0/d) \approx y_0/d = {O}(1)$, and ${\boldsymbol{\mathsf{Q}}}^{(1) tt}_{asym}$ is of the same order of magnitude as the first-order mobility term associated with hydrodynamic interaction between the particles in (4.5). The correction due to asymmetry therefore alters the $\epsilon ^1 \kappa ^1$ correction flow produced by a particle significantly as soon as the particles are placed away from the plane of symmetry of the interface-like gradient by a non-negligible fraction of $d$. In figure 5, the correction flow field produced by a particle located at $y_0 = 0.5 d$ with respect to the constant-viscosity flow field at the local viscosity as a reference is displayed. Comparison to the correction flow in a symmetric viscosity gradient (figure 2) shows that an asymmetry has a significant impact on the particle interactions, which is of the same order of magnitude as the correction in the symmetric case. As a general rule, we find that the effect from asymmetry leads to a stronger correction flow field towards the central plane of the gradient, and to a weaker correction field towards the edge of the gradient.

Figure 5. The $\epsilon ^1 \kappa ^1$ correction far field flow produced by a spherical particle subject to a force and placed asymmetrically in the interface-like viscosity gradient at $y_0 = 0.5 d$, where the local viscosity at the particle is taken as reference. In (a), the force acts along the $y$-direction; in (b), it acts along the $x$-direction. The colour indicates the ratio between the $\epsilon ^1 \kappa ^1$ correction flow and the velocity $U^{(0)} = F_1/(6 {\rm \pi}\eta ^{(0)} a)$ at which the particle would move in constant viscosity $\eta ^{(0)}$ when subject to the force $\boldsymbol {F}_1$.

Interestingly, the leading-order correction to the mobility matrix is affected even if the symmetry of the viscosity gradient is broken at distances of ${O}(d)$ from the particles, while the viscosity profile at distances of ${O}(s)$ is linear and thus odd symmetric near the particles. This results ultimately from the scaling properties of (4.12) with respect to the gradient size $d$: rescaling the coordinate system by $d$, the viscosity profile in terms of the rescaled position is independent of $d$, while the remaining integrand factors out $1/d^4$ and the volume element scales as $d^3$. Consequently, the overall integral scales as $1/d \sim \boldsymbol {\nabla } \eta ^{(1)}$ (see § D.1 of Appendix D for details). As a result, the hydrodynamics in linear viscosity gradients includes highly non-local effects proportional to the magnitude of the viscosity gradient.

5. Interaction of hot particles

As a second application of the presented method, we calculate the interaction of spherical particles with a homogeneous surface temperature that is different from the temperature of the surrounding fluid. In contrast to the previous case of a linear viscosity gradient assumed to be unperturbed by the presence of particles, the viscosity profile here emerges as direct consequence of the particles’ presence. Hence we first calculate the viscosity profile, and obtain from this the correction to the mobility matrix in a second step.

The dependence of fluid viscosity on temperature can often be approximated using Andrade's equation, $\eta = A \exp (b/T)$ (Andrade Reference Andrade1930), where $A$ and $b$ are fluid-specific parameters. However, for small temperature contrasts, Taylor expanding Andrade's equation or one of its generalisations (Gutmann & Simmons Reference Gutmann and Simmons1952) yields a linear relation between temperature and fluid viscosity, which we will use in the following. As a result, a single hot particle induces a viscosity perturbation decaying with inverse distance from the particle centre (figure 6), which in turn alters the self-mobility of the particle as shown in Oppenheimer et al. (Reference Oppenheimer, Navardi and Stone2016). However, to the best of our knowledge, no results are available for the interaction of two hot particles, which we determine in the following.

Figure 6. Exemplary sketch of the viscosity perturbation in the presence of a hot particle and a cold particle. The fluid viscosity is smaller in the neighbourhood of the hot particle, and higher in the neighbourhood of the cold particle.

For a two-particle system with prescribed homogeneous surface temperatures, we first assume that a steady state is established by thermal conduction. Consequently, our model requires that the steady state establishes on faster time scales than the particle temperature changes. In practice, this is the case for beads with large heat capacities. The same would occur for objects subject to constant heat radiation if the susceptibility of the particles is different from that of the surrounding fluid. Second, our approach requires a small Péclet number ($Pe$), which compares advective to thermal heat transport:

(5.1)\begin{equation} Pe = \frac{U a}{\kappa} = Re \, Pr, \end{equation}

where $U$ is the typical velocity in the fluid and $\kappa$ is the thermal diffusivity; $Re$ and $Pr$ denote the Reynolds and Prandtl numbers, and the latter depends only on properties of the fluid (Leal Reference Leal2007, p. 596).

The conduction-induced temperature field around a particle with homogeneous surface temperature in a quiescent fluid, which decays as $r^{-1}$, typically breaks down at a length scale of ${O}(a / Pe)$ (Leal Reference Leal2007, p. 604). For the subsequent calculation, we assume that $s \ll a/ Pe$, i.e. the distance between the two interacting particles is smaller than the length scale at which the temperature profile in the fluid becomes advection-dominated. Similarly to the single particle mobilities (Oppenheimer et al. Reference Oppenheimer, Navardi and Stone2016), the corresponding correction due to finite $Pe$ then scales asymptotically as ${O}(Pe^2)$ and can be neglected in the limit of small Péclet numbers. Indeed, for a micron-sized particle, for example, moving at a speed of microns per second in water at room temperature, the Reynolds number is ${O}(10^{-6})$, and with $Pr \sim {O}(10)$ (Dincer & Zamfirescu Reference Dincer and Zamfirescu2015, Appendix B), the Péclet number is ${O}(10^{-5})$, such that the above assumptions are indeed valid e.g. for colloidal particles (Lu & Weitz Reference Lu and Weitz2013) separated by up to tens of radii.

Under these assumptions, the temperature field in the fluid is described by the heat equation with prescribed homogeneous temperatures at the particle surfaces ($T_1$ and $T_2$ respectively) and at infinity ($T_0$) as boundary conditions:

(5.2a–c)\begin{equation} \nabla^2 T(\boldsymbol{r}) = 0, \quad T(\boldsymbol{r}) = T_b\ \mathrm{for} \ \boldsymbol{r} \in S_b, \quad T(\boldsymbol{r}) \to T_0\ \mathrm{for} \ |\boldsymbol{r}| \to \infty, \end{equation}

where $S_b$ is the surface of particle $b$. For a system consisting of two or more particles, the temperature field can be approximated by a reflection scheme, where monopolar, dipolar, etc. single-particle fields are superposed to match the boundary conditions up to successively increasing order in the inverse particle separation, similarly to the fluid velocity field. We obtain for the temperature field including all terms up to ${O}(s^{-2})$,

(5.3)\begin{align} T(\boldsymbol{r}) = T_0 + C^m_1\,\frac{a}{|\boldsymbol{r} - \boldsymbol{r}_1|} + \frac{a^2 \boldsymbol{C}^d_1 \boldsymbol{\cdot} (\boldsymbol{r} - \boldsymbol{r}_1)}{|\boldsymbol{r} - \boldsymbol{r}_1|^3} + C^m_2\,\frac{a}{|\boldsymbol{r} - \boldsymbol{r}_2|} + \frac{a^2 \boldsymbol{C}^d_2 \boldsymbol{\cdot} (\boldsymbol{r} - \boldsymbol{r}_2)}{|\boldsymbol{r} - \boldsymbol{r}_2|^3} + {O}(s^{{-}3}), \end{align}

with

(5.4a,b)$$\begin{gather} C^m_1 = {\rm \Delta} T_1 \left(1 + \frac{a^2}{s^2} \right) - {\rm \Delta} T_2\, \frac{a}{s}, \quad \boldsymbol{C}^d_1 ={-} {\rm \Delta} T_2\, \frac{a^2 \boldsymbol{s}}{s^3}, \end{gather}$$

(5.5a,b)$$\begin{gather}C^m_2 = {\rm \Delta} T_2 \left( 1 + \frac{a^2}{s^2} \right) - {\rm \Delta} T_1\,\frac{a}{s}, \quad \boldsymbol{C}^d_2 = {\rm \Delta} T_1\,\frac{a^2 \boldsymbol{s}}{s^3}. \end{gather}$$

Approximating the dependence of the viscosity on the temperature as linear, we find

(5.6)\begin{equation} \eta \approx \eta^{(0)} + \left. \frac{\partial \eta}{\partial T} \right|_{T_0} (T - T_0)=: \eta^{(0)} + \alpha (T - T_0), \end{equation}

where $\alpha$ is the proportionality coefficient. For particles with a temperature contrast relative to the fluid, it is natural to assume that $\eta ^{(0)}$ is the viscosity at temperature $T_0$. Assuming that ${\rm \Delta} T_1$ and ${\rm \Delta} T_2$ are typically of the order of some ${\rm \Delta} T$, we define the perturbation parameter $\epsilon := - \alpha \,{\rm \Delta} T/\eta ^{(0)}$ as the relative change in viscosity associated with ${\rm \Delta} T$ compared to $\eta ^{(0)}$ (Oppenheimer et al. Reference Oppenheimer, Navardi and Stone2016). Typically, $\alpha$ is negative, associating hot particles ($T > T_0$) with a decrease in viscosity, and cold particles ($T < T_0$) with an increase, thus $\epsilon$ is positive (Gutmann & Simmons Reference Gutmann and Simmons1952).

Using (5.6) and the definition of $\epsilon$, the $\epsilon ^1$ viscosity perturbation is given as

(5.7)\begin{equation} \eta^{(1)} = -\eta^{(0)}\,\frac{T(\boldsymbol{r}) - T_0}{{\rm \Delta} T}. \end{equation}

By inserting the expression for the temperature field (5.3), we obtain the explicit expansion of the viscosity field:

(5.8) \begin{align} \eta^{(1)}(\boldsymbol{r}) &= \eta^{(1) m}_1\,\frac{a}{|\boldsymbol{r} - \boldsymbol{r}_1|} + \frac{a^2 \boldsymbol{\eta}^{(1) d}_1 \boldsymbol{\cdot} (\boldsymbol{r} - \boldsymbol{r}_1)}{|\boldsymbol{r} - \boldsymbol{r}_1|^3} + \eta^{(1) m}_2\,\frac{a}{|\boldsymbol{r} - \boldsymbol{r}_2|}\nonumber\\ &\quad + \frac{a^2 \boldsymbol{\eta}^{(1) d}_2 \boldsymbol{\cdot} (\boldsymbol{r} - \boldsymbol{r}_2)}{|\boldsymbol{r} - \boldsymbol{r}_2|^3} + {O}(s^{{-}3}), \end{align}

with $\eta ^{(1) m}_b = -(\eta ^{(0)} C^m_b)/{\rm \Delta} T$ and $\boldsymbol {\eta }^{(1) d}_b = -(\eta ^{(0)} \boldsymbol {C}^d_b)/{\rm \Delta} T$, $b \in \{1, 2\}$.

Using (5.8), we first calculate the respective ${\boldsymbol{\mathsf{Q}}}^{(1)}$ terms (Appendix E) and from this compute the $\epsilon ^1$ correction to the mobility matrix. One can show that, up to ${O}(s^{-2})$, the corrections to the mobility matrix are given by the corresponding ${\boldsymbol{\mathsf{Q}}}^{(1)}$ terms (§ E.3 of Appendix E):

(5.9)\begin{equation} \boldsymbol{\mu}^{(1)PQ}_{bc} = {\boldsymbol{\mathsf{Q}}}^{(1)PQ}_{bc} + {O}(s^{{-}3}). \end{equation}

We then find the $\epsilon ^1$ correction to the mobility matrix given by

(5.10) \begin{align} \left. \begin{array}{c@{}} \boldsymbol{\mu}^{(1)tt}_{bb} = \dfrac{1}{6 {\rm \pi}a \eta^{(0)}} \left[ \left( \left(\dfrac{5}{12} - \dfrac{7 a^2}{12 s^2} \right) \dfrac{{\rm \Delta} T_b}{{\rm \Delta} T} + \left( \dfrac{7 a}{12 s} - \dfrac{9 a^2}{8 s^2} \right) \dfrac{{\rm \Delta} T_c}{{\rm \Delta} T} \right) {\boldsymbol{\mathsf{I}}} + \dfrac{9 a^2}{8}\, \dfrac{{\rm \Delta} T_c}{{\rm \Delta} T}\,\dfrac{\boldsymbol{s} \otimes \boldsymbol{s}}{s^4} \right]\\ \hspace{-21pc}+\, {O}(s^{{-}3}), \\ \boldsymbol{\mu}^{(1)rr}_{bb} = \dfrac{1}{6 {\rm \pi}a \eta^{(0)}} \left[ \left(\dfrac{9}{16 a^2} - \dfrac{3}{16 s^2} \right) \dfrac{{\rm \Delta} T_b}{{\rm \Delta} T} + \dfrac{3}{16 as}\,\dfrac{{\rm \Delta} T_c}{{\rm \Delta} T} \right] {\boldsymbol{\mathsf{I}}} + {O}(s^{{-}3}), \\ \boldsymbol{\mu}^{(1) tr}_{11} = \dfrac{1}{6 {\rm \pi}a \eta^{(0)}} \left( -\dfrac{{\rm \Delta} T_2}{{\rm \Delta} T}\,\dfrac{a (\boldsymbol{s} \times)}{8 s^3} \right) + {O}(s^{{-}3}) ={-}\boldsymbol{\mu}^{(1) rt}_{11},\\ \boldsymbol{\mu}^{(1) tr}_{22} = \dfrac{1}{6 {\rm \pi}a \eta^{(0)}} \left( \dfrac{{\rm \Delta} T_1}{{\rm \Delta} T}\, \dfrac{a (\boldsymbol{s} \times)}{8 s^3} \right) + {O}(s^{{-}3}) ={-}\boldsymbol{\mu}^{(1) rt}_{22},\\ \boldsymbol{\mu}^{(1)tt}_{bc} = \dfrac{1}{6 {\rm \pi}a \eta^{(0)}} \left( \dfrac{9 a^2}{8}\,\dfrac{{\rm \Delta} T_b + {\rm \Delta} T_c}{{\rm \Delta} T}\,\dfrac{\boldsymbol{s} \otimes \boldsymbol{s}}{s^4} \right) + {O}(s^{{-}3}), \quad b \neq c, \\ \boldsymbol{\mu}^{(1)tr}_{bc} = \boldsymbol{\mu}^{(1) rt}_{bc} = \boldsymbol{\mu}^{(1) rr}_{bc} = {O}(s^{{-}3}), \quad b \neq c . \end{array}\right\} \end{align}

The corrections to the translational and rotational self-mobilities of, say, particle 1 depend at order $s^0$ only on the particle's own temperature, which alters the self-mobility similarly as in the case when there is only one particle in the fluid (Oppenheimer et al. Reference Oppenheimer, Navardi and Stone2016). If the particle is hot, i.e. $T_1 > T_0$, then the particle's mobility increases, and it decreases for $T_1 < T_0$.

In contrast to the case of a single particle, we find additional terms scaling as $s^{-1}$ in the first particle's translational and rotational self-mobilities, which are proportional to the temperature of the second particle at the distance $s$. This is because the temperature of particle 2 induces a radially decaying temperature and viscosity field around it, which can be approximated to leading order as constant near particle 1. In order to satisfy the boundary condition for the temperature field at the surface of particle 1, a reflected monopolar temperature field around particle 1 is induced. However, the reflected flow field decays as $r^{-1}$, where $r$ is the distance from particle 1. Therefore, the effect due to the incident field is stronger than the effect of the reflected field, and the mobility of particle 1 indeed increases for a hot particle 2, and decreases if particle 2 has a lower temperature than the fluid at infinity. Consequently, the ratio of magnitudes of the correction to self-mobility due to the presence of the second particle and the single-particle self-mobility of a neutral particle scales as $\epsilon a/s$.

While for a single particle with homogeneous surface temperature translation and rotation do not couple (Oppenheimer et al. Reference Oppenheimer, Navardi and Stone2016), such coupling is observed at ${O}(s^{-2})$ in the presence of a second particle. This effect on, say, particle 1 can be understood as the result of the gradient of the monopolar temperature field induced by particle 2. However, also here a dipolar temperature field reflected from particle 1 arises to satisfy the boundary conditions at the surface of particle 1, and weakens the effect of the incident temperature gradient.

Considering the components associated with hydrodynamic interactions between particles in (5.10), only the translational term is non-vanishing at our level of approximation, and it scales as $1/s^2$. Interestingly, the translational interaction between two particles increases proportional to the sum ${\rm \Delta} T_1 + {\rm \Delta} T_2$, i.e. it is irrelevant whether the particle causing or experiencing the interaction has an altered temperature. This observation is in agreement with the fact that the mobility matrix is symmetric as a result of the Lorentz reciprocal theorem (Jeffrey & Onishi Reference Jeffrey and Onishi1984), a result that holds independently of whether or not the viscosity is constant. Moreover, if one particle is hotter and the other is colder than the fluid by the same amount, then the effects cancel and the interaction between the particles is, up to order ${O}(s^{-2})$, similar to the interaction of two neutral particles.

Considering the second particle as a neutral test particle, i.e. ${\rm \Delta} T_2 = 0$, for the far field correction flow produced by particle 1 located at $\boldsymbol {r} = \boldsymbol {0}$ with increased temperature ${\rm \Delta} T_1 > 0$ and subject to a force $\boldsymbol {F}_1$, we find

(5.11)\begin{equation} \boldsymbol{u}^{(1)} (\boldsymbol{r}) = \boldsymbol{U}_2^{(1)} = \frac{1}{6 {\rm \pi}a \eta^{(0)}}\, \frac{9 a^2 {\rm \Delta} T_1}{8 r^4 {\rm \Delta} T}\,\boldsymbol{r} (\boldsymbol{r} \boldsymbol{\cdot} \boldsymbol{F}_1). \end{equation}

This flow field decays with the second power of the distance from particle 1 and is depicted in figure 7. Along the direction of the force applied, the correction flow is parallel to the force, thus enhancing the constant viscosity Stokeslet flow. This is reasonable since with increased temperature, the viscosity around the particle is decreased and the flow field becomes stronger. However, on the axis perpendicular to the force, the correction flow vanishes, concentrating the correction flow along the direction of the force applied. As a consistency check, we find that the divergence of the flow field (5.11) vanishes, thus the flow satisfies the incompressibility condition. In contrast to a linear viscosity gradient, this correction flow field decays more strongly than the constant-viscosity Stokeslet by one order. The ratio of the first-order correction in the mobility term for particle interaction and the constant-viscosity term thus also scales as $\epsilon a/s$, and hence can be important in application.

Figure 7. The $\epsilon ^1$ far field correction flow field produced by a hot particle with a force applied along the $x$-direction. The colour indicates the ratio of the first-order correction flow and the velocity $U^{(0)} = F_1/( 6 {\rm \pi}\eta ^{(0)} a)$, at which the particle would move if it had the same temperature as the surrounding fluid, $T_1 = T_0$, and was subject to the force $\boldsymbol {F}_1$.

In comparison to the case of linear viscosity gradients (§ 4), three-body effects affect the correction to the mobility matrix already at the level of approximation considered, namely at ${O}(s^{-2})$. For instance, a third particle 3 will in general induce a reflected temperature and viscosity field in response to the temperature of particle 2. At the position of particle 1, this reflected field then scales as $1/s^2$, assuming that all inter-particle distances are of ${O}(s)$. Consequently, the self-mobility of particle 1 will be altered at ${O}(s^{-2})$ by a proper three-body effect. Even more importantly, the correction to the mobility term associated with interaction between two particles 1 and 2 will also be altered at ${O}(s^{-2})$ by the temperature of a nearby particle 3. This is because particle 3 heats or cools the fluid in between particles 1 and 2 in a similar way as particles 1 and 2 do.