  We present a versatile perturbative calculation scheme to determine the
mobility matrix for two and more particles in a low Reynolds number fluid with
spatially variant viscosity. Assuming an asymptotic non-constant viscosity
perturbation superimposed on a constant viscosity background, we exploit the
Lorentz reciprocal theorem and a reflection method to obtain the leading order
correction to the mobility matrix. We apply our approach firstly to an
interface-like linear viscosity gradient with an extension much larger than the
length scale of the particle separation. We find that the viscosity gradient
gives rise to additional long-range flow fields, for both particle translation
and rotation, which decay by one order slower than their constant-viscosity
counterparts. Also, we reveal that the relative positioning of two interacting
particles in finite-size gradients, a natural precondition to avoid
pathological negative viscosity, affects the mobilities of and interaction
between both particles. Secondly, we apply our result to two interacting
particles with temperatures different from that of the surrounding fluid, which
alters both the fluid temperature and viscosity field. We find that both, the
self-mobilities of the particles as well as their interactions, increase for
hot particles and decrease for cold particles.
