# Trapped-particle microrheology of active suspensions

Zhiwei Peng  John F. Brady<EMAIL_ADDRESS>Division of Chemistry and
Chemical Engineering, California Institute of Technology, Pasadena, California
91125, USA

###### Abstract

In microrheology, the local rheological properties such as viscoelasticity of
a complex fluid are inferred from the free or forced motion of embedded
colloidal probe particles. Theoretical machinery developed for forced-probe
microrheology of colloidal suspensions focused on either constant-force (CF)
or constant-velocity (CV) probes while in experiments neither the force nor
the kinematics of the probe is fixed. More importantly, the constraint of CF
or CV introduces a difficulty in the meaningful quantification of the
fluctuations of the probe due to a thermodynamic uncertainty relation. It is
known that for a Brownian particle trapped in a harmonic potential well, the
product of the standard deviations of the trap force and the particle position
is $dk_{B}T$ in $d$ dimensions with $k_{B}T$ being the thermal energy. As a
result, if the force (position) is not allowed to fluctuate, the position
(force) fluctuation becomes infinite. To allow the measurement of
fluctuations, in this work we consider a microrheology model in which the
embedded probe is dragged along by a moving harmonic potential so that both
its position and the trap force are allowed to fluctuate. Starting from the
full Smoluchowski equation governing the dynamics of $N$ hard active Brownian
particles, we derive a pair Smoluchowski equation describing the dynamics of
the probe as it interacts with one bath particle by neglecting hydrodynamic
interactions among particles in the dilute limit. From this, we determine the
mean and the variance (i.e., fluctuation) of the probe position in terms of
the pair probability distribution. We then characterize the behavior of the
system in the limits of both weak and strong trap. By taking appropriate
limits, we show that our generalized model can be reduced to the well-studied
CF or CV microrheology models.

##  I Introduction

Rheology is the study of flow and deformation of complex materials in response
to an applied force. Traditional (bulk) rheological measurements are performed
by shearing a macroscopic sample of the material confined between two solid
surfaces, such as in the cone-and-plate rheometer. Bulk rheological studies
such as shear rheometry provide a measurement of the macroscopic rheological
behavior of complex materials.

Recently, particle-tracking microrheology has become a standard tool for
studying the mechanical properties of materials on a much smaller scale
(Weihs, Mason, and Teitell, 2006; Cicuta and Donald, 2007; Wirtz  _et al._ ,
2009; Furst and Squires, 2017). In contrast to bulk rheology, microrheology
only requires a small sample volume and can be used to quantify spatial
heterogeneity. As a result, microrheology is particularly useful for examining
soft biological materials. For example, classical bulk rheometry cannot be
used to probe the microenvironment inside living cells without disrupting
their mechanical structure while particle-tracking microrheology can be
performed (Wilhelm, Gazeau, and Bacri, 2003; Nawaz  _et al._ , 2012; Berret,
2016; Ayala  _et al._ , 2016; Hu  _et al._ , 2017).

To aid in the understanding of experimental measurements and in the prediction
of colloidal microrheology, Squires and Brady (2005) developed a theoretical
framework in which a colloidal probe is pulled through a suspension of
neutrally buoyant bath colloids. This model has been used and generalized to
study the microrheology of passive colloids (Khair and Brady, 2005, 2006;
Meyer  _et al._ , 2006; Zia and Brady, 2010; Swan and Zia, 2013; Zia, 2018)
and active colloids (Burkholder and Brady, 2019, 2020). When the external
pulling force is absent, the probe “collides” with bath particles as it
undergoes Brownian motion—the so-called tracer diffusion problem. To
characterize the nonlinear response, forced microrheology is considered in
which an external force, often larger than the thermodynamic restoring force,
is applied to the probe. Within forced microrheology, two operating
modes—constant-force (CF) and constant-velocity (CV)—are often considered from
a theoretical perspective. In the CF mode, the probe is driven by a constant
external force $\bm{F}^{\text{ext}}$ and the velocity of the probe is
fluctuating. Conversely, for a CV probe, the probe velocity $\bm{U}_{1}$ is a
constant vector (Therefore, the position of the probe is known at all times.)
and the force required to maintain such a steady motion must fluctuate.

To characterize the micro-viscous response of colloidal suspensions, an
effective microviscosity $\eta^{\text{eff}}$ can be defined using the Stokes
drag law. For a spherical probe of radius $a$ in the CF mode, this is given by
$F^{\text{ext}}=6\pi\eta^{\text{eff}}a\langle U_{1}\rangle$, where $\langle
U_{1}\rangle$ the probe velocity in the direction of $\bm{F}^{\text{ext}}$
averaged over Brownian fluctuations. The ratio between the effective
microviscosity and the solvent viscosity, $\eta^{\text{eff}}/\eta$, is the
main quantity of interest in colloidal microrheology. For the CV mode, the
average external force is used in the definition of the effective
microviscosity: $\langle F^{\text{ext}}\rangle=6\pi\eta^{\text{eff}}aU_{1}$.
In order to measure the microviscoelastic response of suspensions, an
oscillatory driving force is considered (Khair and Brady, 2005).

While the CF (or CV) model is successful in quantifying the mean velocity (or
mean force) of a probe driven through colloidal suspensions. The fluctuation
from this mean value is largely unexplored. Taking the CV mode as an example,
one could calculate the variance of the mean force using the probe-distorted
microstructure. The question is what does this variance physically imply? In
particular, how does this variance relate to the fluctuations in the
suspension? In an experimental setting, neither the force nor the velocity of
the probe is fixed; they are both allowed to fluctuate (Meyer  _et al._ ,
2006; Weihs, Mason, and Teitell, 2006; Cicuta and Donald, 2007; Yao  _et al._
, 2009).

To mimic the experimental realization more closely and motivate later
discussions, consider the simple case of an isolated Brownian particle in a
harmonic trap that is centered at the origin (arbitrary). In this physical
picture, both the position and the velocity of the particle is fluctuating. A
statistical mechanical description can be adopted in which one defines the
probability density, $P(\bm{r},t)$, of finding the particle at position
$\bm{r}$ relative to the fixed trap at time $t$. Conservation of probability
dictates that $P(\bm{r},t)$ is governed by the Smoluchowski equation, which
reads $\partial P/\partial t+\nabla\cdot\bm{j}=0$, where the flux vector
$\bm{j}=P\bm{F}^{\text{trap}}/\zeta-D_{T}\nabla P$. Here,
$\bm{F}^{\text{trap}}$ is the trap force and for a harmonic trap is given by
$\bm{F}^{\text{trap}}=-k\bm{r}$ with $k$ being the spring constant; $\zeta$ is
the drag coefficient and $D_{T}$ is the thermal diffusivity given by the
Stokes-Einstein-Sutherland relation, $\zeta D_{T}=k_{B}T$, where $k_{B}T$ is
the thermal energy. The mean external force exerted on the Brownian particle
is
$\langle\bm{F}^{\text{trap}}\rangle=\int\bm{F}^{\text{trap}}Pd\bm{r}=-k\int\bm{r}Pd\bm{r}=-k\langle\bm{r}\rangle$.
Because the trap is harmonic, the mean force is proportional to the mean
displacement with $-k$ being the constant of proportionality. For a fixed
trap, the mean position (therefore the mean force) is zero,
$\langle\bm{r}\rangle=\bm{0}$. The variance of the force,
$\operatorname{Var}\left(\bm{F}^{\text{trap}}\right)=k^{2}\operatorname{Var}(\bm{r})$.
A straightforward calculation leads to the result

$\operatorname{Var}(\bm{r})=\frac{k_{B}T}{k}\bm{I},$ (1)

where $\bm{I}$ is the identity tensor. Introducing the shorthand
$\Delta\bm{F}^{\text{trap}}=\bm{F}^{\text{trap}}-\langle\bm{F}^{\text{trap}}\rangle$,
we can write the fluctuation relation as

$\left<(\Delta\bm{F}^{\text{trap}})^{2}\right>^{1/2}\left<(\Delta\bm{r}_{1})^{2}\right>^{1/2}=dk_{B}T,$
(2)

where $d$ is the spatial dimensionality.

Equation (2) is a fundamental result and a few comments on its implications
are in order. First, by harmonically trapping a particle immersed in a
solvent, the product of the standard deviations of the trap force and the
particle position gives precisely the thermal fluctuations of the
solvent—$dk_{B}T$. Second, one can decrease the uncertainty in the position by
increasing the stiffness of the trap [see equation (1)]. However, the trade-
off is that the fluctuation in the force must increase due to (2). Said
differently, this constitutes a thermodynamic uncertainty relation in which
one cannot decrease the fluctuations in both the force and the position
simultaneously. If the fluctuation in the position vanishes (infinitely stiff
trap), the fluctuation in the force blows up.

We note that (2) is observed elsewhere. For example, consider an ideal
Gaussian polymer chain with one end localized in a harmonic trap. The
fluctuations of the trap force and the position from the trap center satisfy
an identical relation (Wang, 2017).

We are now in a position to consider the fluctuations in the microrheology
problem. Instead of considering either CF or CV, we must allow both the
position of and the force on the probe to fluctuate in order to have a
meaningful quantification of fluctuations. Equation (2) also implies that we
should consider the position not the velocity of the probe. In the CF mode,
therefore, the quantity of interest for fluctuations is the variance of the
position of the probe, which is just the force-induced tracer diffusion
problem. That is, the tracer diffusivity under the influence of a constant
force should be considered—not the variance of the velocity. For the CV mode,
the position of the probe is also prescribed and the fluctuation in the force
is infinite. As a result, in the CV mode the computed variance of the force
does not have a physical meaning.

In this paper, to closely mimic the setup of microrheological experiments, we
consider a trapped-particle microrheology model in which the colloidal probe
particle is driven by a translating harmonic trap. Because biological
materials examined by microrheology such as the microenvironment inside living
cells often contain active “particles”, we model the suspension as an active
colloidal suspension. Compared to passive suspensions, the study of the
microrheology of active suspensions is more recent (Jepson  _et al._ , 2013;
Miño  _et al._ , 2013; Morozov and Marenduzzo, 2014; Kasyap, Koch, and Wu,
2014; Reichhardt and Reichhardt, 2015; Burkholder and Brady, 2017, 2019, 2020;
Knežević, Avilés Podgurski, and Stark, 2021; Peng and Brady, 2022; Seyforth
_et al._ , 2022). The colloidal particles in an active suspension are able to
self-propel, which can be a model for either biologically active microswimmers
or synthetic phoretic particles. This active colloidal suspension model also
includes passive (not self-propelled) colloidal systems, which can be obtained
by setting the self-propulsive swim speed to zero.

The paper is organized as follows. In section II, we present the general
$N$-particle dynamics from a continuum perspective using the Smoluchowski
equation governing the evolution of the positions and orientations of $N$
active Brownian particles. In section III we first derive the mean and
variance (fluctuation) of the probe position relative to the trap center from
the $N$-particle formulation. Neglecting hydrodynamic interactions in the
dilute limit, we then derive the pair-level Smoluchowski equation governing
the dynamics of the probe and one bath particle. We discuss the asymptotic
behavior of the system in the limits of both weak and strong traps. We then
show in section IV that our generalized theoretical framework includes the
well-studied CF and CV microrheology models when appropriate limits are taken.
Finally, we conclude in section V.

##  II Mechanics of active Brownian suspensions

Consider a colloidal suspension consisting of $N$ particles dispersed in an
incompressible Newtonian fluid (solvent) of dynamic viscosity $\eta$. The
particles could be active and are subject to fluctuating thermal (Brownian)
forces from the solvent. Furthermore, the inertia of the fluid and the
particles are assumed to be negligible. In this low-Reynolds-number regime,
the fluid dynamics is governed by the linear Stokes equations and the
probability distribution of the particles are described by the Smoluchowski
equation. In general, all $N$ particles could be active, and we model them as
active Brownian particles. The probability distribution for finding the $N$
particles in positions $\\{\bm{x}_{\alpha}\\}$ and orientations
$\\{\bm{q}_{\alpha}\\}$ at a given time $t$ is denoted as
$P_{N}(\bm{x}^{N},\bm{q}^{N},t)$ where $\alpha=1,\cdot\cdot\cdot,N$ is the
particle label. In the laboratory frame of reference, the $N$-particle
Smoluchowski equation is given by

$\frac{\partial P_{N}}{\partial
t}+\sum_{\alpha=1}^{N}\nabla_{\alpha}^{T}\cdot\bm{j}_{\alpha}^{T}+\sum_{\alpha=1}^{N}\nabla_{\alpha}^{R}\cdot\bm{j}_{\alpha}^{R}=0,$
(3)

where $\nabla_{\alpha}^{T}=\partial/\partial\bm{x}_{\alpha}$ is the spatial
gradient operator with respect to the position vector ($\bm{x}_{\alpha}$) of
particle $\alpha$ in the laboratory frame and
$\nabla_{\alpha}^{R}=\bm{q}_{\alpha}\times\left(\partial/\partial\bm{q}_{\alpha}\right)$
is the orientational gradient operator of particle $\alpha$. The translational
and rotational fluxes in equation (3) are, respectively, given by
$\bm{j}_{\alpha}^{T}=\bm{U}_{\alpha}P_{N}$ and
$\bm{j}_{\alpha}^{R}=\bm{\varOmega}_{\alpha}P_{N}$, where $\bm{U}_{\alpha}$
($\bm{\varOmega}_{\alpha}$) is the instantaneous linear (angular) velocity of
particle labeled $\alpha$ relative to the laboratory frame. The conservation
of probability is

$\int_{\Gamma_{N}}P_{N}d\Gamma^{N}=1,$ (4)

where $d\Gamma^{N}=\prod_{\alpha=1}^{N}d\Gamma_{\alpha}$ denotes the volume
element of the $N$-particle phase space and
$d\Gamma_{\alpha}=d\bm{x}_{\alpha}d\bm{q}_{\alpha}$ is the volume element of
the phase space of particle $\alpha$.

In the absence of a background flow, the linear and angular velocities of any
active particle $\alpha$ are given by

$\displaystyle\begin{pmatrix}\bm{U}_{\alpha}-\bm{U}_{\alpha}^{0}\\\
\bm{\varOmega}_{\alpha}-\bm{\varOmega}_{\alpha}^{0}\end{pmatrix}=$
$\displaystyle\sum_{\beta=1}^{N}\bm{\mathcal{M}}_{\alpha\beta}\cdot\begin{pmatrix}\bm{F}_{\beta}^{e}+\bm{F}_{\beta}^{P}-k_{B}T\nabla_{\beta}^{T}\ln
P_{N}\\\ \bm{L}_{\beta}^{e}+\bm{L}_{\beta}^{P}-k_{B}T\nabla_{\beta}^{R}\ln
P_{N}\end{pmatrix}$ $\displaystyle+\begin{pmatrix}\bm{0}\\\
-D_{\alpha}^{R}\nabla_{\alpha}^{R}\ln P_{N}\end{pmatrix},$ (5)

where $\bm{\mathcal{M}}_{\alpha\beta}$ is the configuration-dependent grand
hydrodynamic mobility tensor coupling the linear and angular velocity of
particle $\alpha$ to the force and torque exerted on particle $\beta$. Note
that for general particle shapes $\bm{\mathcal{M}}_{\alpha\beta}$ is a
function of the instantaneous $N$-particle configuration—both positions and
orientations. The forces on any particle $\beta$ include the external force
$\bm{F}_{\beta}^{e}$, the interparticle colloidal force $\bm{F}_{\beta}^{P}$
and the thermal or entropic force $-k_{B}T\nabla_{\beta}^{T}\ln P_{N}$.
Similarly, the torques on any particle $\beta$ include the external torque
$\bm{L}_{\beta}^{e}$, the interparticle colloidal torque $\bm{L}_{\beta}^{P}$
and the thermal torque $-k_{B}T\nabla_{\beta}^{R}\ln P_{N}$. The interparticle
colloidal forces and torques are assumed to be conservative. For the case of
hard-sphere interactions, the interparticle forces reduce to no-flux boundary
conditions at any surface of contact between particles.

In equation (II), the activity of any particle $\alpha$ is modeled by its
undisturbed swim linear velocity $\bm{U}_{\alpha}^{0}$ and angular velocity
$\bm{\varOmega}_{\alpha}^{0}$ regardless of the presence of any other
particles. For the case of simple ABPs, the swim angular velocity is often
taken to be zero, $\bm{\varOmega}_{\alpha}^{0}=\bm{0}$. Furthermore, a
biological microswimmer may “decide” to change its orientation
$\bm{q}_{\alpha}$ by, for example, actuating the flagella on a different side
of its body without disturbing the flow. In this process, the body of the
microswimmer does not turn. For non-spherical particles, this process means
that the swim orientation $\bm{q}_{\alpha}$ is usually different from the
orientation of the particle shape, in which case the shape orientation needs
to be included as an additional phase space variable. For spherical particles,
only the swim orientation matters and no such difficulty is introduced. This
reorientation process of any particle $\alpha$ is independent of the motion of
other particles and is modeled by a simple rotary diffusion with a constant
rotary diffusivity $D_{\alpha}^{R}$. The reorientation time is
$\tau_{\alpha}^{R}=1/D_{\alpha}^{R}$, which defines the active run or
persistence length of an ABP: $\ell_{\alpha}=U_{\alpha}^{0}\tau_{\alpha}^{R}$.
Because this reorientation process is biological rather than thermal in
origin, $D_{\alpha}^{R}$ is not constrained by the fluctuation-dissipation
theorem and may be inferred from experimental data.

##  III Moving-trap microrheology

In the context of microrheology, the particle with label $1$ is identified as
the probe particle. This particle could be a new particle placed into the
suspension or one of the suspension particles tagged as the probe. Particles
labeled $2-N$ are referred to as bath particles. In the following, we consider
a suspension of neutrally buoyant, hard and active colloidal spheres with
identical radii. The probe may have a different radius than the bath
particles. Instead of fixing the external force $\bm{F}_{1}^{e}$ or the
velocity $\bm{U}_{1}$, the probe particle is trapped in a translating harmonic
potential well. Denoting the position vector of the center of the potential
well as $\bm{x}_{0}(t)$, we have $d\bm{x}_{0}/dt=\bm{U}^{\mathrm{trap}}(t)$,
where $\bm{U}^{\mathrm{trap}}(t)$ is the prescribed velocity of the moving
trap relative to the laboratory frame. The trap force $\bm{F}_{1}^{e}$ is
assumed to be only a function of the relative position between the probe and
the potential well. All bath particles experience no external forces or
torques. We first consider a general derivation in which all particles are
ABPs and the probe is a tagged ABP in the suspension.

In the constant-force or constant-velocity mode of microrheology, the position
of the probe does not matter, and the system is statistically homogeneous. In
contrast, the introduction of a moving trap defines a specific origin in the
system and the position of the probe relative to the trap needs to be
considered explicitly. To this end, we first change to a coordinate system
moving with the instantaneous trap velocity and measure all particle positions
relative to the trap. This change of variables is written as
$\bm{z}_{\alpha}=\bm{z}_{\alpha}(\\{\bm{x}\\},t)=\bm{x}_{\alpha}-\int_{0}^{t}\bm{U}^{\mathrm{trap}}(s)ds-\bm{x}_{0}(0)$
for any $\alpha$ and $t^{\prime}=t^{\prime}(\\{\bm{x}\\},t)=t$. Using the
chain rule we obtain $\partial/\partial
t=-\sum_{\alpha=1}^{N}\bm{U}^{\mathrm{trap}}\cdot\partial/\partial\bm{z}_{\alpha}+\partial/\partial
t^{\prime}$ and
$\partial/\partial\bm{x}_{\alpha}=\partial/\partial\bm{z}_{\alpha}$. The
Smoluchowski equation (3) in the new coordinate system becomes

$\frac{\partial P_{N}}{\partial
t^{\prime}}+\sum_{\alpha=1}^{N}\frac{\partial}{\partial\bm{z}_{\alpha}}\cdot\left(\bm{j}_{\alpha}^{T}-\bm{U}^{\mathrm{trap}}P_{N}\right)+\sum_{\alpha=1}^{N}\nabla_{\alpha}^{R}\cdot\bm{j}_{\alpha}^{R}=0,$
(6)

where $\bm{j}_{\alpha}^{T}$ and $\bm{j}_{\alpha}^{R}$ remain unchanged. In the
context of microrheology, it is more convenient to measure the positions of
all bath particles relative to that of the probe. We therefore introduce
another change of variables such that for the probe
$\bm{r}_{1}=\bm{r}_{1}(\bm{z}^{N},t^{\prime})=\bm{z}_{1}$, and
$\bm{r}_{\alpha}=\bm{r}_{\alpha}(\bm{z}^{N},t^{\prime})=\bm{z}_{\alpha}-\bm{z}_{1}$
for all bath particles ($\alpha=2,\cdot\cdot\cdot,N$). In this coordinate
system, the probe position is measured relative to the trap and the positions
of all bath particles are measured relative to the probe. The change of
variables allows us to write
$\partial/\partial\bm{z}_{1}=\partial/\partial\bm{r}_{1}-\sum_{\alpha=2}^{N}\partial/\partial\bm{r}_{\alpha}$
and $\partial/\partial\bm{z}_{\alpha}=\partial/\partial\bm{r}_{\alpha}$ for
$\alpha=2,\cdot\cdot\cdot,N$. The Smoluchowski equation (6) transforms to

$\displaystyle\frac{\partial P_{N}}{\partial
t}+\nabla_{1}^{T}\cdot\left(\bm{j}_{1}^{T}-\bm{U}^{\mathrm{trap}}P_{N}\right)$
$\displaystyle+\sum_{\alpha=2}^{N}\nabla_{\alpha}^{T}\cdot\left(\bm{j}_{\alpha}^{T}-\bm{j}_{1}^{T}\right)$
$\displaystyle+\sum_{\alpha=1}^{N}\nabla_{\alpha}^{R}\cdot\bm{j}_{\alpha}^{R}=0.$
(7)

It is understood that in equation (III) we have used $t$ for the time variable
and $\nabla_{\alpha}^{T}=\partial/\partial\bm{r}_{\alpha}$ for any $\alpha$.
Formally, the probability density in equation (III) is the conditional
probability of find all particles at a given configuration provided that the
trap is at $\bm{x}_{0}$ at time $t$, i.e.,
$P_{N}=P_{N}\left(\bm{r}^{N},\bm{q}^{N},t\rvert\bm{x}_{0},t\right)$. The
translational flux of particle $\alpha$ can be written as

$\displaystyle\bm{j}_{\alpha}^{T}$ $\displaystyle=$ $\displaystyle
U_{\alpha}^{0}\bm{q}_{\alpha}P_{N}+\bm{M}_{\alpha
1}^{UF}\cdot\bm{F}_{1}^{e}P_{N}$ (8)
$\displaystyle-\sum_{\beta=1}^{N}\left(\bm{D}_{\alpha\beta}^{UF}-\bm{D}_{\alpha
1}^{UF}\right)\cdot\nabla_{\beta}^{T}P_{N}$ $\displaystyle-\bm{D}_{\alpha
1}^{UF}\cdot\nabla_{1}^{T}P_{N}-\sum_{\beta=1}^{N}\bm{D}_{\alpha\beta}^{UL}\cdot\nabla_{\beta}^{R}P_{N},$

where we have taken $\bm{U}_{\alpha}^{0}=U_{\alpha}^{0}\bm{q}_{\alpha}$ and
used the Stokes-Einstein-Sutherland relations
$\bm{D}_{\alpha\beta}^{UF}=k_{B}T\bm{M}_{\alpha\beta}^{UF},\bm{D}_{\alpha\beta}^{UL}=k_{B}T\bm{M}_{\alpha\beta}^{UL}$.
For all accessible configurations, the inter-particle forces are zero and the
hard-particle interaction between two spheres do not induce torques.
Similarly, the rotary flux of particle $\alpha$ is given by

$\displaystyle\bm{j}_{\alpha}^{R}$ $\displaystyle=$
$\displaystyle\bm{M}_{\alpha 1}^{\Omega
F}\cdot\bm{F}_{1}^{e}P_{N}-\sum_{\beta=1}^{N}\left(\bm{D}_{\alpha\beta}^{\Omega
F}-\bm{D}_{\alpha 1}^{\Omega F}\right)\cdot\nabla_{\beta}^{T}P_{N}$ (9)
$\displaystyle-\bm{D}_{\alpha 1}^{\Omega
F}\cdot\nabla_{1}^{T}P_{N}-\sum_{\beta=1}^{N}\bm{D}_{\alpha\beta}^{\Omega
L}\cdot\nabla_{\beta}^{R}P_{N}-D_{\alpha}^{R}\nabla_{\alpha}^{R}P_{N}.$

There are no external force or torque on the bath particles,
$\alpha=2\mbox{--}N$, nor a torque on the probe, $\bm{L}_{1}^{e}=\bm{0}$.

The Smoluchowski equation (III) together with the flux expressions (8) and (9)
fully specify the $N$-particle phase space dynamics. Some comments regarding
equations (III)-(9) are in order. First, the above derivation is an extension
of the model considered by Squires and Brady (2005) for passive Brownian
suspensions. We have generalized their model to a suspension of ABPs in which
one of the particles is tagged as the probe that is driven by a translating
trap. Realizing that the grand mobility tensor does not depend on the swim
orientation vectors of spherical particles, one can set
$\bm{U}_{\alpha}^{0}=0$ and integrate over the orientations of all particles
to obtain the trapped probe microrheology problem of a passive Brownian
suspension. Note that even for passive suspensions, if the probe or the bath
particles are non-spherical, their shape orientations need to be included in
the above formulation. Second, the hydrodynamic interactions between all
$N$-particles are included in the grand mobility tensor. In particular, this
leads to the fact that a gradient in orientation space of particle $\beta$
induces a translational flux of particle $\alpha$, and vice versa, due to the
hydrodynamic translation-rotation coupling. Third, due to the dependence on
particle orientations, the phase space of $N$ ABPs has a dimension of $5N$:
the physical space has a dimension of $3N$ and the orientation space has a
dimension of $2N$ if the orientation of each particle is parametrized by the
azimuthal and polar angles of a spherical coordinate system.

###  III.1 Mean and fluctuation of the probe position

The average position or mean displacement of the probe relative to the trap is
defined by

$\langle\bm{r}_{1}\rangle(t)=\int\bm{r}_{1}P_{N}d\Gamma^{N},$ (10)

where the angle bracket denotes integration against $P_{N}$ over the
configuration space of all particles. Multiplying equation (III) by
$\bm{r}_{1}$ and integrating over the configuration space $\Gamma^{N}$, we
obtain

$\displaystyle\frac{\partial\langle\bm{r}_{1}\rangle}{\partial
t}+\bm{U}^{\mathrm{trap}}$ $\displaystyle=$ $\displaystyle
U_{1}^{0}\langle\bm{q}_{1}\rangle+\bigl{\langle}\bm{M}_{11}^{UF}\cdot\bm{F}_{1}^{e}\bigr{\rangle}-\bigl{\langle}\bm{D}_{11}^{UF}\cdot\nabla_{1}^{T}\ln
P_{N}\bigr{\rangle}$ (11)
$\displaystyle-\sum_{\beta=1}^{N}\bigl{\langle}\left(\bm{D}_{1\beta}^{UF}-\bm{D}_{11}^{UF}\right)\cdot\nabla_{\beta}^{T}\ln
P_{N}\bigr{\rangle}.$

Similarly, the mean squared displacement, a second order tensor, is governed
by

$\frac{\partial\langle\bm{r}_{1}\bm{r}_{1}\rangle}{\partial
t}+2\left[\bm{U}^{\mathrm{trap}}\langle\bm{r}_{1}\rangle\right]^{\mathrm{sym}}=2\int\left[\bm{j}_{1}^{T}\bm{r}_{1}\right]^{\mathrm{sym}}d\Gamma^{N},$
(12)

where the integral

$\displaystyle\int\bm{j}_{1}^{T}\bm{r}_{1}d\Gamma^{N}$ $\displaystyle=$
$\displaystyle
U_{1}^{0}\langle\bm{q}_{1}\bm{r}_{1}\rangle+\Bigl{\langle}\bm{M}_{11}^{UF}\cdot\bm{F}_{1}^{e}\bm{r}_{1}\Bigr{\rangle}$
(13)
$\displaystyle-\sum_{\beta=1}^{N}\Bigl{\langle}\left(\bm{D}_{1\beta}^{UF}-\bm{D}_{11}^{UF}\right)\cdot\left(\nabla_{\beta}^{T}\ln
P_{N}\right)\bm{r}_{1}\Bigr{\rangle}$
$\displaystyle-\Bigl{\langle}\bm{D}_{11}^{UF}\cdot\left(\nabla_{1}^{T}\ln
P_{N}\right)\bm{r}_{1}\Bigr{\rangle},$

and the superscript “sym” denotes the symmetric part of a tensor (see equation
(34)).

The main quantities of interest in the present problem are the mean
displacement $\langle\bm{r}_{1}\rangle$ and the fluctuation

$\displaystyle\operatorname{Var}(\bm{r}_{1})$
$\displaystyle=\operatorname{Cov}(\bm{r}_{1},\bm{r}_{1})=\bigl{\langle}\Delta\bm{r}_{1}\Delta\bm{r}_{1}\bigr{\rangle}$
$\displaystyle=\bigl{\langle}\bm{r}_{1}\bm{r}_{1}\bigr{\rangle}-\langle\bm{r}_{1}\rangle\langle\bm{r}_{1}\rangle,$
(14)

where we have introduced the shorthand
$\Delta\bm{r}_{1}=\bm{r}_{1}-\langle\bm{r}_{1}\rangle$ and
$\operatorname{Var}(\bm{r}_{1})$ denotes the variance tensor of $\bm{r}_{1}$.
For a harmonic trap, the mean force is related to the mean displacement via

$\langle\bm{F}_{1}^{e}\rangle=-k\langle\bm{r}_{1}\rangle,$ (15)

and similarly the fluctuation in the force is given by

$\operatorname{Var}(\bm{F}_{1}^{e})=\bigl{\langle}\Delta\bm{F}_{1}^{e}\Delta\bm{F}_{1}^{e}\bigr{\rangle}=k^{2}\langle\Delta\bm{r}_{1}\Delta\bm{r}_{1}\rangle.$
(16)

###  III.2 The pair problem

To proceed analytically, we restrict the analysis to the dilute limit in which
only pair interactions between a bath particle and the probe is considered.
Furthermore, we neglect hydrodynamic interactions between the bath particle
and the probe, and only consider hard-sphere interactions. The reduction from
the $N$-particle formulation to the pair problem and the consideration of
hydrodynamic interactions are discussed in appendix A.

Because the bath particles are indistinguishable, it is convenient to define
the two-particle probability density function
$\rho_{2}(\bm{r}_{2},\bm{q}_{2},\bm{r}_{1},\bm{q}_{1},t)$, which denotes the
joint probability density function of finding the probe at
$(\bm{r}_{1},\bm{q}_{1})$ and any bath particle at $(\bm{r}_{2},\bm{q}_{2})$
at time $t$. In terms of
$P_{2}(\bm{r}_{2},\bm{q}_{2},\bm{r}_{1},\bm{q}_{1},t)$, which is the joint
probability density function of finding the probe at $(\bm{r}_{1},\bm{q}_{1})$
and the bath particle labeled $2$ (i.e., the first bath particle) at
$(\bm{r}_{2},\bm{q}_{2})$ at time $t$, we have $\rho_{2}=(N-1)P_{2}$. Here,
the factor of $N-1$ comes from removing the labels from the $N-1$ bath
particles. The joint probability can be written as

$\displaystyle\rho_{2}$
$\displaystyle=\rho_{1/1}(\bm{r}_{2},\bm{q}_{2},t\rvert\bm{r}_{1},\bm{q}_{1},t)P_{1}(\bm{r}_{1},\bm{q}_{1},t)$
$\displaystyle=n_{b}g_{1/1}(\bm{r}_{2},\bm{q}_{2},t\rvert\bm{r}_{1},\bm{q}_{1},t)P_{1}(\bm{r}_{1},\bm{q}_{1},t),$
(17)

where $n_{b}=(N-1)/V$ is the number density of bath particles. For a passive
and CF (or CV) probe, $g_{1/1}$ becomes independent of the configuration
($\bm{r}_{1}$ and $\bm{q}_{1}$) of the probe due to statistical homogeneity;
in this case the probe distribution $P_{1}$ can be integrated over and one
only needs to consider $g_{1/1}$ (Squires and Brady, 2005; Burkholder and
Brady, 2018, 2019).

The joint probability $\rho_{2}$ (see appendix A) is governed by

$\displaystyle\frac{\partial\rho_{2}}{\partial t}$
$\displaystyle+\nabla_{1}^{T}\cdot\left(\bm{j}_{1}^{T}-\bm{U}^{\mathrm{trap}}\rho_{2}\right)+\nabla_{2}^{T}\cdot\left(\bm{j}_{2}^{T}-\bm{j}_{1}^{T}\right)$
$\displaystyle+\sum_{\alpha=1}^{2}\nabla_{\alpha}^{R}\cdot\bm{j}_{\alpha}^{R}=0,$
(18)

where

$\displaystyle\bm{j}_{1}^{T}=U_{1}^{0}\bm{q}_{1}\rho_{2}+\frac{1}{\zeta_{1}}\bm{F}_{1}^{e}\rho_{2}+D_{1}^{T}\nabla_{2}^{T}\rho_{2}-D_{1}^{T}\nabla_{1}^{T}\rho_{2},$
(19)
$\displaystyle\bm{j}_{2}^{T}=U_{2}^{0}\bm{q}_{2}\rho_{2}-D_{2}^{T}\nabla_{2}^{T}\rho_{2},$
(20)
$\displaystyle\bm{j}_{\alpha}^{R}=-D_{\alpha}^{R}\nabla_{\alpha}^{R}\rho_{2}.$
(21)

Figure 1: Schematic of the pair problem of a spherical probe particle in a
moving harmonic trap interacting with a spherical bath particle. Both the
probe and the bath particles can be active.

At contact, $r_{2}=R_{c}$, no relative flux is allowed:

$\bm{n}_{2}\cdot\left(\bm{j}_{2}^{T}-\bm{j}_{1}^{T}\right)=0.$ (22)

Far away from the probe, the bath distribution is undisturbed by the probe and
the probe distribution is that in the absence of the bath particles,

$\rho_{2}(\bm{r}_{2},\bm{q}_{2},\bm{r}_{1},\bm{q}_{1},t)\to\frac{n_{b}}{\Omega_{b}}P_{1}(\bm{r}_{1},\bm{q}_{1},t)\quad\mathrm{as}\quad\lvert\bm{r}_{2}\rvert\to\infty,$
(23)

where $\Omega_{b}$ is the total solid angle of the orientation space of the
bath particle. In 3D, $\Omega_{b}=4\pi$. Far away from the trap, the
probability vanishes

$\rho_{2}\to 0\quad\mathrm{as}\quad\lvert\bm{r}_{1}\rvert\to\infty.$ (24)

Equation (11) governing the mean displacement becomes

$\displaystyle\frac{\partial\langle\bm{r}_{1}\rangle}{\partial
t}+\frac{1}{\tau_{k}}\langle\bm{r}_{1}\rangle=$
$\displaystyle-\bm{U}^{\mathrm{trap}}+U_{1}^{0}\langle\bm{q}_{1}\rangle$
$\displaystyle+D_{1}^{T}\int\nabla_{2}^{T}\rho_{2}d\Gamma^{2},$ (25)

where $d\Gamma^{2}=d\Gamma_{1}d\Gamma_{2}$, and we have defined the
viscoelastic timescale

$\tau_{k}=\frac{\zeta_{1}}{k},$ (26)

which is set by the balance between the viscous force
$\zeta_{1}\partial\langle\bm{r}_{1}\rangle/\partial t$ and the elastic force
$k\langle\bm{r}_{1}\rangle$. Using the divergence theorem and the far-field
condition (23), the last term on the rhs of (III.2) can be written as

$D_{1}^{T}\int\nabla_{2}^{T}\rho_{2}d\Gamma^{2}=D_{1}^{T}\int
d\bm{q}_{2}d\Gamma_{1}\oint_{S_{c}}\bm{n}_{2}\rho_{2}dS_{2},$ (27)

where $S_{c}=\\{\bm{r}_{2}:\lvert\bm{r}_{2}\rvert=R_{c}\\}$ is the contact
surface and $\bm{n}_{2}$ is the unit normal vector of $S_{c}$ that points out
of particle $2$.

As shown in appendix B, the position fluctuation of the probe is governed by

$\frac{1}{2}\frac{\partial\operatorname{Var}(\bm{r}_{1})}{\partial
t}+\frac{1}{\tau_{k}}\operatorname{Var}(\bm{r}_{1})=D_{1}^{T}\bm{I}+\left[U_{1}^{0}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})+D_{1}^{T}\int\Delta\bm{r}_{1}\nabla_{2}^{T}\rho_{2}d\Gamma^{2}\right]^{\mathrm{sym}},$
(28)

where the covariance of $\bm{q}_{1}$ and $\bm{r}_{1}$ satisfies

$\displaystyle\frac{\partial\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})}{\partial
t}+\frac{1}{\tau}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})=U_{1}^{0}\operatorname{Var}(\bm{q}_{1})+D_{1}^{T}\int\Delta\bm{q}_{1}\nabla_{2}^{T}\rho_{2}d\Gamma^{2}.$
(29)

In equation (29), we have defined the relaxation time $\tau$ using

$\displaystyle\frac{1}{\tau}=\frac{1}{\tau_{k}}+\frac{d-1}{\tau_{1}^{R}}.$
(30)

Regardless of the presence of the trap or the bath particles, at long times
($t\to\infty$) the net polar and nematic orders of the probe are given by
$\langle\bm{q}_{1}\rangle=\bm{0}$ and
$\bigl{\langle}\bm{q}_{1}\bm{q}_{1}\bigr{\rangle}=\bm{I}/d$, respectively (see
appendix B). As a result, $\operatorname{Var}(\bm{q}_{1})=\bm{I}/d$ at long
times.

It is convenient to consider the rank $m$ polyadic spatial moment tensor

$\bm{M}_{m}(\bm{r}_{2},\bm{q}_{2},\bm{q}_{1},t)=\int\underbrace{\bm{r}_{1}\cdot\cdot\cdot\bm{r}_{1}}_{m}\rho_{2}d\bm{r}_{1},\quad(m=0,1,2,...).$
(31)

Multiplying equation (III.2) by the $m$-adic product of $\bm{r}_{1}$ and
integrating over the physical space of the probe, we obtain

$\displaystyle\frac{\partial\bm{M}_{m}}{\partial
t}-m\left[U_{1}^{0}\bm{q}_{1}\bm{M}_{m-1}-\frac{k}{\zeta_{1}}\bm{M}_{m}+(m-1)D_{1}^{T}\bm{M}_{m-2}\bm{I}+D_{1}^{T}\nabla_{2}^{T}\bm{M}_{m-1}-\bm{U}^{\mathrm{trap}}\bm{M}_{m-1}\right]^{\mathrm{sym}}$
$\displaystyle+\nabla_{2}^{T}\cdot\left(\bm{U}_{r}^{0}\bm{M}_{m}-D_{r}^{T}\nabla_{2}^{T}\bm{M}_{m}+\frac{k}{\zeta_{1}}\bm{M}_{m+1}\right)-mD_{1}^{T}\left[\nabla_{2}^{T}\bm{M}_{m-1}\right]^{\mathrm{sym}}-\sum_{\alpha=1}^{2}D_{\alpha}^{R}\nabla_{\alpha}^{R}\cdot\nabla_{\alpha}^{R}\bm{M}_{m}=\bm{0},$
(32)

where we have defined the relative swim velocity and the relative diffusivity
as, respectively,

$\bm{U}_{r}^{0}=U_{2}^{0}\bm{q}_{2}-U_{1}^{0}\bm{q}_{1},\quad
D_{r}^{T}=D_{1}^{T}+D_{2}^{T},$ (33)

whereas $\left[\bm{A}\right]^{\mathrm{sym}}$ denotes the symmetric part of any
rank $m$ Cartesian tensor $\bm{A}$ such that

$\left[\bm{A}\right]^{\mathrm{sym}}_{i_{1}i_{2}\cdot\cdot\cdot
i_{m}}=\frac{1}{m!}\sum_{\sigma\in\mathfrak{S}_{m}}A_{i_{\sigma 1}i_{\sigma
2}\cdot\cdot\cdot i_{\sigma m}},$ (34)

in which $\mathfrak{S}_{m}$ is the set containing the $m!$ permutations of
indices. For $m=2$, this reduces to the familiar definition of the symmetric
part of a rank $2$ tensor,
$\bm{A}^{\mathrm{sym}}=\left(\bm{A}+\bm{A}^{\intercal}\right)/2$. For any rank
$m$ tensor $\bm{A}$, its symmetric part $\left[\bm{A}\right]^{\mathrm{sym}}$
is invariant under a permutation of all indices. In equation (32),
$\bm{M}_{m}$ for $m<0$ is understood to be zero.

At contact, $r_{2}=R_{c}$, the no-flux boundary condition is satisfied:

$\displaystyle\bm{n}_{2}\cdot\left(\bm{U}_{r}^{0}\bm{M}_{m}-D_{r}^{T}\nabla_{2}^{T}\bm{M}_{m}+\frac{k}{\zeta_{1}}\bm{M}_{m+1}\right)$
$\displaystyle-
mD_{1}^{T}\left[\bm{n}_{2}\bm{M}_{m-1}\right]^{\mathrm{sym}}=\bm{0}.$ (35)

The far-field condition for the spatial moment of rank $m$ is

$\bm{M}_{m}\to\frac{n_{b}}{\Omega_{b}}\bm{\Phi}_{m}(\bm{q}_{1},t)\quad\mathrm{as}\quad
r_{2}\to\infty,$ (36)

where

$\bm{\Phi}_{m}(\bm{q}_{1},t)=\int\underbrace{\bm{r}_{1}\cdot\cdot\cdot\bm{r}_{1}}_{m}P_{1}(\bm{r}_{1},\bm{q}_{1},t)d\bm{r}_{1}$
(37)

is the rank $m$ spatial moment of the single particle probability $P_{1}$ of
the probe. Discussion of the single particle behavior and the method to obtain
$\bm{\Phi}_{m}$ is deferred to section III.3.

From (III.2) and (B), to obtain the mean and mean-squared displacements, one
only needs to calculate the zeroth and first spatial moments, respectively. On
the other hand, the definitions of mean and mean-squared displacements allow
us to write
$\langle\bm{r}_{1}\rangle=\int\bm{M}_{1}/(N-1)d\bm{q}_{1}d\Gamma_{2}$ and
$\langle\bm{r}_{1}\bm{r}_{1}\rangle=\int\bm{M}_{2}/(N-1)d\bm{q}_{1}d\Gamma_{2}$.
Because in obtaining $\langle\bm{r}_{1}\bm{r}_{1}\rangle$ only the integral of
$\bm{M}_{2}$ is required, it’s not necessary to first calculate the
distribution of $\bm{M}_{2}$ explicitly before carrying out the integration.
Instead, one can show that integrating equation (32) for $m=2$ leads to the
same equation as (B). Due to the presence of the harmonic trap force, the
equation for $\bm{M}_{m}$ is coupled to $\bm{M}_{m+1}$. To truncate this
infinite set of equations and obtain a finite set of closed equations, a
closure model may be used.

To see the structure of the spatial moments more clearly, we write out the
first few moment equations explicitly using (32). The zeroth moment,
$M_{0}=\int\rho_{2}d\bm{r}_{1}$, satisfies the equation

$\displaystyle\frac{\partial M_{0}}{\partial t}$
$\displaystyle+\nabla_{2}^{T}\cdot\left[\bm{U}_{r}^{0}M_{0}-D_{r}^{T}\nabla_{2}^{T}M_{0}+\frac{k}{\zeta_{1}}\bm{M}_{1}\right]$
$\displaystyle-\sum_{\alpha=1}^{2}D_{\alpha}^{R}\nabla_{\alpha}^{R}\cdot\nabla_{\alpha}^{R}M_{0}=0,$
(38)

and the normalization $\int M_{0}d\bm{q}_{1}d\Gamma_{2}=N-1$. In addition to
being advected by the relative velocity $\bm{U}_{r}^{0}$ in the physical space
of the bath particle, $M_{0}$ is forced by the trap via the divergence of the
first moment.

The equation governing the evolution of the first spatial moment is

$\displaystyle\frac{\partial\bm{M}_{1}}{\partial
t}-\left(U_{1}^{0}\bm{q}_{1}M_{0}-\frac{k}{\zeta_{1}}\bm{M}_{1}+D_{1}^{T}\nabla_{2}^{T}M_{0}\right)+\bm{U}^{\mathrm{trap}}M_{0}$
$\displaystyle+\nabla_{2}^{T}\cdot\left[\bm{U}_{r}^{0}\bm{M}_{1}-D^{T}_{r}\nabla_{2}^{T}\bm{M}_{1}+\frac{k}{\zeta_{1}}\bm{M}_{2}-D_{1}^{T}\bm{I}M_{0}\right]$
$\displaystyle-\sum_{\alpha=1}^{2}D_{\alpha}^{R}\nabla_{\alpha}^{R}\cdot\nabla_{\alpha}^{R}\bm{M}_{1}=\mathbf{0}.$
(39)

Similarly, the second moment is governed by

$\displaystyle\frac{\partial\bm{M}_{2}}{\partial t}$
$\displaystyle-2\left[U_{1}^{0}\bm{q}_{1}\bm{M}_{1}-\frac{k}{\zeta_{1}}\bm{M}_{2}+D_{1}^{T}M_{0}\bm{I}+D_{1}^{T}\nabla_{2}^{T}\bm{M}_{1}-\bm{U}^{\mathrm{trap}}\bm{M}_{1}\right]^{\mathrm{sym}}$
$\displaystyle+\nabla_{2}^{T}\cdot\left[\bm{U}_{r}^{0}\bm{M}_{2}-D_{r}^{T}\nabla_{2}^{T}\bm{M}_{2}+\frac{k}{\zeta_{1}}\bm{M}_{3}\right]-2D_{1}^{T}\left[\nabla_{2}^{T}\bm{M}_{1}\right]^{\mathrm{sym}}-\sum_{\alpha=1}^{2}D_{\alpha}^{R}\nabla_{\alpha}^{R}\cdot\nabla_{\alpha}^{R}\bm{M}_{2}=\mathbf{0}.$
(40)

###  III.3 The probe distribution in the absence of bath particles

The simplest problem in the above formulation is that of a single particle
(the probe) interacting with the trap. One can formulate this single-particle
problem by neglecting all bath particles or taking the limit $\phi_{b}=4\pi
b^{3}n_{b}/3\to 0$ in the above $N$-particle formulation. The single-particle
probability $P_{1}(\bm{r}_{1},\bm{q}_{1},t\rvert\bm{x}_{0},t)$ of the active
probe satisfies

$\displaystyle\frac{\partial P_{1}}{\partial t}$
$\displaystyle+\nabla_{1}^{T}\cdot\left(\frac{1}{\zeta_{1}}\bm{F}_{1}^{e}P_{1}-D_{1}^{T}\nabla
P_{1}-\bm{U}^{\mathrm{trap}}P_{1}+U_{1}^{0}\bm{q}_{1}P_{1}\right)$
$\displaystyle-D_{1}^{R}\nabla_{1}^{R}\cdot\nabla_{1}^{R}P_{1}=0,$ (41)

where the conservation of probability dictates that $\int P_{1}d\Gamma_{1}=1$
and the harmonic trap force $\bm{F}_{1}^{e}=-k\bm{r}_{1}$. We emphasize that
in equation (III.3) the probe is also considered as an ABP.

The rank $m$ ($m=0,1,...$) spatial moment of $P_{1}$ defined by (37) satisfies

$\displaystyle\frac{\partial\bm{\Phi}_{m}}{\partial t}$
$\displaystyle-m\left[U_{1}^{0}\bm{q}_{1}\bm{\Phi}_{m-1}-\frac{k}{\zeta_{1}}\bm{\Phi}_{m}+(m-1)D_{1}^{T}\bm{\Phi}_{m-2}\bm{I}-\bm{U}^{\mathrm{trap}}\bm{\Phi}_{m-1}\right]^{\mathrm{sym}}$
$\displaystyle-
D_{1}^{R}\nabla_{1}^{R}\cdot\nabla_{1}^{R}\bm{\Phi}_{m}=\bm{0},$ (42)

where $\bm{\Phi}_{m}$ for $m<0$ is defined to be zero. Different from equation
(32) in which the moment $\bm{M}_{m}$ is coupled to $\bm{M}_{m+1}$, the rank
$m$ spatial moment of $P_{1}$ only depends on lower order moments, which leads
to a set of closed equations. The solution to the preceding equation provides
the far-field condition for $\bm{M}_{m}$ as given by equation (36).

The zeroth-order spatial moment $\Phi_{0}$ is the net orientational
distribution, which is unaffected by the trap and is governed by the
orientational diffusion equation:

$\frac{\partial\Phi_{0}}{\partial
t}-D_{1}^{R}\nabla_{1}^{R}\cdot\nabla_{1}^{R}\Phi_{0}=0,$ (43)

where the conservation of $P_{1}$ gives $\int\Phi_{0}d\bm{q}_{1}=1$. At long
times, the solution is simply the uniform distribution,
$\Phi_{0}(\bm{q}_{1},t\to\infty)=1/(4\pi)$ in 3D.

The above formulation also allows us to consider the mean and fluctuation of
the probe displacement in the absence of bath particles. Equation (11) or
(III.2) in the absence of bath particles reduces to

$\frac{\partial\langle\bm{r}_{1}\rangle}{\partial
t}+\frac{1}{\tau_{k}}\langle\bm{r}_{1}\rangle=-\bm{U}^{\mathrm{trap}}+U_{1}^{0}\langle\bm{q}_{1}\rangle,$
(44)

where for the single particle
$\langle\bm{r}_{1}\rangle=\int\bm{r}_{1}P_{1}d\Gamma_{1}=\int\bm{\Phi}_{1}d\bm{q}_{1}$.
Similarly, equation (12) or (B) for the single particle becomes

$\frac{1}{2}\frac{\partial\langle\bm{r}_{1}\bm{r}_{1}\rangle}{\partial
t}+\frac{1}{\tau_{k}}\bigl{\langle}\bm{r}_{1}\bm{r}_{1}\bigr{\rangle}=D_{1}^{T}\bm{I}+\left[U_{1}^{0}\bigl{\langle}\bm{q}_{1}\bm{r}_{1}\bigr{\rangle}-\bm{U}^{\mathrm{trap}}\langle\bm{r}_{1}\rangle\right]^{\mathrm{sym}}.$
(45)

It can be seen from equations (44) and (45) that in order to calculate the
mean and mean-squared displacements, one needs to obtain the net polar order
$\langle\bm{q}_{1}\rangle$ and the covariance of the position and orientation
$\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})$. The governing equation for
$\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})$ follows from (29) and is given by

$\frac{\partial\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})}{\partial
t}+\frac{1}{\tau}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})=U_{1}^{0}\operatorname{Var}(\bm{q}_{1})$
(46)

which depends on the net nematic order $\langle\bm{q}_{1}\bm{q}_{1}\rangle$.

At steady state, it is shown that $\langle\bm{q}_{1}\rangle=\bm{0}$ and
$\langle\bm{q}_{1}\bm{q}_{1}\rangle=\bm{I}/d$, where $d=2,3$ is the
dimensionality of the physical space. This allows us to obtain

$\displaystyle\langle\bm{r}_{1}\rangle$
$\displaystyle=-\frac{\zeta_{1}\bm{U}^{\mathrm{trap}}}{k},$ (47)
$\displaystyle\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})$
$\displaystyle=\frac{U_{1}^{0}}{dk/\zeta_{1}+d(d-1)D_{1}^{R}}\bm{I},$ (48)
$\displaystyle\langle\bm{r}_{1}\bm{r}_{1}\rangle=$
$\displaystyle\frac{\zeta_{1}^{2}}{k^{2}}\bm{U}^{\mathrm{trap}}\bm{U}^{\mathrm{trap}}+\frac{\zeta_{1}D_{1}^{T}}{k}\bm{I}$
$\displaystyle+\frac{\zeta_{1}D_{1}^{\mathrm{swim}}}{k}\frac{1}{1+\frac{k\tau_{1}^{R}}{\zeta_{1}}\frac{1}{d-1}}\bm{I},$
(49)

where $D_{1}^{\mathrm{swim}}=\left(U_{1}^{0}\right)^{2}\tau_{1}^{R}/[d(d-1)]$
is the swim diffusivity of a freely swimming ABP. The average position of the
ABP relative to the trap is given by the balance between the average trap
force $k\langle\bm{r}_{1}\rangle$ and the viscous drag
$\zeta_{1}\bm{U}^{\mathrm{trap}}$. If the trap is strong, $k\to\infty$, the
ABP is tightly confined and pushing against the trap ‘boundary’, which has
been observed in experiments (Takatori  _et al._ , 2016). On the other hand,
for $k\to 0$, the average position of the ABP becomes unbounded. Solving the
steady state first and then taking the limit $k\to 0$ in (47) is singular
because in the absence of the trap ($k=0$) the average position is unbounded
and at long times the particle motion is diffusive. For $k\equiv 0$, we are
simply measuring the motion of an ABP in a frame of reference moving with
velocity $\bm{U}^{\mathrm{trap}}$ relative to the laboratory frame, which
gives $d\langle\bm{r}_{1}\rangle/dt=-\bm{U}^{\mathrm{trap}}$.

Takatori  _et al._ (2016) studied the transient and long-time dynamics of
self-propelled Janus particles in a fixed acoustic trap. They showed that the
experimentally measured density distribution of Janus particles follow closely
the theoretical predictions using a harmonic trap. Equation (III.3) in the
absence of $\bm{U}^{\mathrm{trap}}$ agrees with that obtained in Takatori  _et
al._ (2016).

The fluctuation relation is given by

$\bigl{\langle}\left(\Delta\bm{F}_{1}^{e}\right)^{2}\bigr{\rangle}^{1/2}\bigl{\langle}\left(\Delta\bm{r}_{1}\right)^{2}\bigr{\rangle}^{1/2}=d\left[k_{B}T+\frac{k_{s}T_{s}}{1+\tau_{1}^{R}/[(d-1)\tau_{k}]}\right],$
(50)

where the thermal energy $k_{B}T=\zeta_{1}D_{1}^{T}$ and analogously an
_active_ energy scale $k_{s}T_{s}$ has been defined such that
$k_{s}T_{s}=\zeta_{1}D_{1}^{\mathrm{swim}}$ (Takatori, Yan, and Brady, 2014).
In equation (50), the fluctuation consists of the thermal (passive) energy
$dk_{B}T$ and an active energy contribution. This active energy is different
from $k_{s}T_{s}$ due to the presence of the harmonic trap, which introduces
an orientational decorrelation timescale $\tau_{k}$ in addition to the
reorientation time $\tau_{1}^{R}$ of the ABP. For a weak trap,
$\tau_{1}^{R}/\tau_{k}\ll 1$, the decorrelation occurs on the timescale of
$\tau_{1}^{R}$, and the active contribution scales as
$\zeta_{1}\left(U_{1}^{0}\right)^{2}\tau_{1}^{R}$. As a result, the
fluctuation
$\bigl{\langle}\left(\Delta\bm{F}_{1}^{e}\right)^{2}\bigr{\rangle}^{1/2}\bigl{\langle}\left(\Delta\bm{r}_{1}\right)^{2}\bigr{\rangle}^{1/2}\to
d(k_{B}T+k_{s}T_{s})$ as $\tau_{1}^{R}/\tau_{k}\to 0$. This is often referred
to as Rule #1 of active matter—when all length scales are large compared to
the run length $\ell_{1}$, one can replace $k_{B}T$ with $k_{B}T+k_{s}T_{s}$.
As another example, consider the sedimentation of active colloids under
gravity. At steady state, the number density follows Boltzmann distribution
but with $k_{B}T+k_{s}T_{s}$ in place of $k_{B}T$ (Palacci  _et al._ , 2010).

When $\tau_{1}^{R}/\tau_{k}\gg 1$, the relevant timescale is $\tau_{k}$, and
the active contribution scales as
$\zeta_{1}\left(U_{1}^{0}\right)^{2}\tau_{k}$. In this limit, the ABP is
pushing against the edge of the potential well and the fluctuation comes from
passive Brownian motion alone,
$\bigl{\langle}\left(\Delta\bm{F}_{1}^{e}\right)^{2}\bigr{\rangle}^{1/2}\bigl{\langle}\left(\Delta\bm{r}_{1}\right)^{2}\bigr{\rangle}^{1/2}\to
dk_{B}T$ as $\tau_{1}^{R}/\tau_{k}\to\infty$.

Regardless of the trap strength, the product of the square root of the
fluctuations in the force and the position is always bounded. For a strong
trap, the position fluctuation vanishes,
$\bigl{\langle}\Delta\bm{r}_{1}\Delta\bm{r}_{1}\bigr{\rangle}=O(1/k)\to 0$,
but the force fluctuation blows up linearly since
$\bigl{\langle}\Delta\bm{F}_{1}^{e}\Delta\bm{F}_{1}^{e}\bigr{\rangle}=O(k)\to\infty$
as $k\to\infty$. Conversely, the position fluctuation grows unboundedly while
the force fluctuation vanishes as $k\to 0$.

In the weak trap limit, equation (50) can be equivalently written as

$\frac{k}{\zeta_{1}}\bigl{\langle}\Delta\bm{r}_{1}\Delta\bm{r}_{1}\bigr{\rangle}=\frac{\bigl{\langle}\Delta\bm{r}_{1}\Delta\bm{r}_{1}\bigr{\rangle}}{\tau_{k}}\to\bm{D}_{1}^{\mathrm{eff}}\quad\mathrm{as}\quad\frac{\tau_{1}^{R}}{\tau_{k}}\to
0,$ (51)

where $\bm{D}_{1}^{\mathrm{eff}}=D_{1}^{T}\bm{I}+D_{1}^{\mathrm{swim}}\bm{I}$
is the long-time effective diffusivity of the ABP in the absence of the trap
(see appendix C for the asymptotic analysis). This relation implies the
equivalence of the position fluctuation divided by $\tau_{k}$ in the limit of
vanishing harmonic trapping force and the effective diffusion of a free ABP.
In other words, one could calculate the position fluctuation in a trap and
then take the limit of
$\bigl{\langle}\Delta\bm{r}_{1}\Delta\bm{r}_{1}\bigr{\rangle}/\tau_{k}$ as
$k\to 0$ to obtain the long-time diffusivity that the particle would have in
the absence of the trap, or vice versa. Because the trap is weak, the ABP is
able to explore space via both thermal fluctuation and its undisturbed active
run-and-reorientation, both processes contribute to the position fluctuation.
In the presence of bath particles, this equivalence still holds in which
$\bm{D}_{1}^{\mathrm{eff}}$ is the diffusivity of the probe affected by
collisions with bath particles (i.e., tracer diffusion).

###  III.4 A weak trap

For a weak trap, $\epsilon=\tau_{1}^{R}/\tau_{k}=k\tau_{1}^{R}/\zeta_{1}\ll
1$, the probe is allowed to explore and reorient freely before reaching the
“boundary” of the potential well. The viscoelastic timescale $\tau_{k}$ is
well separated from the reorientation timescale $\tau_{1}^{R}$. In the
intermediate timescale characterized by $t/\tau_{1}^{R}\gg 1$ and
$t/\tau_{k}\ll 1$, the probe has explored the suspension but has not reached
the boundary of the potential; we expect a diffusive behavior of the probe. At
times much longer than the viscoelastic timescale ($t/\tau_{k}\gg 1$), the
variance of the probe position becomes bounded due to the trapping force.
Therefore, the motion of the probe exhibits a transition from diffusive to
bounded behavior.

The separation of the two timescales allows us to consider a multiple-scale
analysis. By defining the fast variable $t_{1}=t$ and the slow variable
$t_{2}=\epsilon t$, we have $\partial/\partial t=\partial/\partial
t_{1}+\epsilon\partial/\partial t_{2}$. Regular perturbation expansions of the
pair probability distribution and its spatial moments in terms of $\epsilon$
are written as

$\displaystyle\rho_{2}$ $\displaystyle=$
$\displaystyle\rho_{2}^{(0)}+\epsilon\rho_{2}^{(1)}+\cdot\cdot\cdot,$ (52)
$\displaystyle\bm{M}_{m}$ $\displaystyle=$
$\displaystyle\bm{M}_{m}^{(0)}+\epsilon\bm{M}_{m}^{(0)}+\cdot\cdot\cdot,$ (53)

where $\bm{M}_{m}^{(k)}$ is the rank $m$ spatial moment of $\rho_{2}^{(k)}$.

At $O(1)$, the zeroth moment satisfies

$\displaystyle\frac{\partial M_{0}^{(0)}}{\partial
t_{1}}+\nabla_{2}^{T}\cdot\left(\bm{U}_{r}^{0}M_{0}^{(0)}-D_{r}^{T}\nabla_{2}^{T}M_{0}^{(0)}\right)$
$\displaystyle-\sum_{\alpha=1}^{2}D_{\alpha}^{R}\nabla_{\alpha}^{R}\cdot\nabla_{\alpha}^{R}M_{0}^{(0)}=0,$
(54)
$\displaystyle\bm{n}_{2}\cdot\left(\bm{U}_{r}^{0}M_{0}^{(0)}-D_{r}^{T}\nabla_{2}^{T}M_{0}^{(0)}\right)=0,\quad\bm{r}_{2}\in
S_{c}.$ (55)

Similarly, the first moment at this order is given by

$\displaystyle\frac{\partial\bm{M}_{1}^{(0)}}{\partial
t_{1}}-\left(U_{1}^{0}\bm{q}_{1}M_{0}^{(0)}+D_{1}^{T}\nabla_{2}^{T}M_{0}^{(0)}\right)+\bm{U}^{\mathrm{trap}}M_{0}^{(0)}$
$\displaystyle+\nabla_{2}^{T}\cdot\left(\bm{U}_{r}^{0}\bm{M}_{1}^{(0)}-D^{T}_{r}\nabla_{2}^{T}\bm{M}_{1}^{(0)}-D_{1}^{T}\bm{I}M_{0}^{(0)}\right)$
$\displaystyle-\sum_{\alpha=1}^{2}D_{\alpha}^{R}\nabla_{\alpha}^{R}\cdot\nabla_{\alpha}^{R}\bm{M}_{1}^{(0)}=\mathbf{0},$
(56)
$\displaystyle\bm{n}_{2}\cdot\left(\bm{U}_{r}^{0}\bm{M}_{1}^{(0)}-D^{T}_{r}\nabla_{2}^{T}\bm{M}_{1}^{(0)}-D_{1}^{T}\bm{I}M_{0}^{(0)}\right)=0,\quad\bm{r}_{2}\in
S_{c}.$ (57)

Expanding the covariances similarly, e.g.,

$\operatorname{Var}(\bm{r}_{1})=\operatorname{Var}^{(0)}(\bm{r}_{1})+\epsilon\operatorname{Var}^{(1)}(\bm{r}_{1})+\cdot\cdot\cdot,$
(58)

we obtain at $O(1)$

$\displaystyle\frac{1}{2}\frac{\partial\operatorname{Var}^{(0)}(\bm{r}_{1})}{\partial
t_{1}}=D_{1}^{T}\bm{I}+\Big{[}U_{1}^{0}\operatorname{Cov}^{(0)}(\bm{q}_{1},\bm{r}_{1})+D_{1}^{T}\int\nabla_{2}^{T}\rho_{2}^{(0)}\Delta\bm{r}_{1}d\Gamma^{2}\Big{]}^{\mathrm{sym}},$
(59)
$\displaystyle\frac{\partial\operatorname{Cov}^{(0)}(\bm{q}_{1},\bm{r}_{1})}{\partial
t_{1}}+\frac{d-1}{\tau_{1}^{R}}\operatorname{Cov}^{(0)}(\bm{q}_{1},\bm{r}_{1})=U_{1}^{0}\operatorname{Var}(\bm{q}_{1})+D_{1}^{T}\int\Delta\bm{q}_{1}\nabla_{2}^{T}\rho_{2}^{(0)}d\Gamma^{2}.$
(60)

Note that $\operatorname{Var}(\bm{q}_{1})$ is not affected by the presence of
the trap (see appendix B) and therefore only has the $O(1)$ term in the
small-$\epsilon$ expansion.

Equations (III.4)–(60) govern the dynamics of a probe in a bath of active
particles in the absence of the trapping force (The presence of
$\bm{U}^{\mathrm{trap}}$ in (III.4) is only due to the fact that we are in a
frame of reference moving with $\bm{U}^{\mathrm{trap}}$ relative to the
laboratory frame). This problem is the so-called tracer—an active
one—diffusion in an active Brownian suspension. Even in the absence of the
trap, the correlation between $\bm{q}_{1}$ and $\bm{r}_{1}$ has a steady-state
(time-independent) solution due to the presence of the decorrelation time
$\tau_{1}^{R}$ in equation (60). Dropping the time derivative in (60) at
steady state, we obtain

$\displaystyle\operatorname{Cov}^{(0)}(\bm{q}_{1},\bm{r}_{1})=$
$\displaystyle\frac{\ell_{1}}{d(d-1)}\bm{I}$
$\displaystyle+\frac{\tau_{1}^{R}}{d-1}D_{1}^{T}\int\Delta\bm{q}_{1}\nabla_{2}^{T}\rho_{2}^{(0)}d\Gamma^{2},$
(61)

where it is understood that the steady-state distribution of $\rho_{2}^{(0)}$
is used, and $\ell_{1}=U_{1}^{0}\tau_{1}^{R}$ is the run length of the active
probe. Therefore, equation (59) is written as

$\displaystyle\bm{D}_{1}^{\mathrm{eff}}=$
$\displaystyle\left(D_{1}^{T}+D_{1}^{\mathrm{swim}}\right)\bm{I}$
$\displaystyle+D_{1}^{T}\left[\int\left(\frac{\ell_{1}}{d-1}\Delta\bm{q}_{1}+\Delta\bm{r}_{1}\right)\nabla_{2}^{T}\rho_{2}^{(0)}d\Gamma^{2}\right]^{\mathrm{sym}},$
(62)

where
$\bm{D}_{1}^{\mathrm{eff}}=\partial\operatorname{Cov}^{(0)}(\bm{r}_{1},\bm{r}_{1})/(2\partial
t_{1})$ is the long-time diffusivity of the probe in the absence of the
trapping force. As expected, one could obtain the same result by setting
$\bm{F}_{1}^{e}=\bm{0}$ from the outset (see section IV.2). This is done in
Burkholder and Brady (2017) but with the free tracer particle being passive.

So long as the trapping force is not identically zero, the probe will
eventually reach the boundary of the trap. This confinement happens at very
large distances from the trap (or at long times if the probe is started near
the trap center).

###  III.5 A strong trap

For a strong trap, the viscoelastic time scale $\tau_{k}=\zeta_{1}/k$ is much
smaller than other timescales (e.g., the reorientation time ) of the problem.
Due to the strong trapping force, both the mean and the variance of the probe
have a steady-state solution that is time independent.

The position fluctuation, governed by equation (28), becomes at steady state

$\displaystyle\frac{k}{\zeta_{1}}\operatorname{Var}(\bm{r}_{1})=$
$\displaystyle
D_{1}^{T}\bm{I}+\Big{[}U_{1}^{0}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})$
$\displaystyle+D_{1}^{T}\int\Delta\bm{r}_{1}\nabla_{2}^{T}\rho_{2}d\Gamma^{2}\Big{]}^{\mathrm{sym}}.$
(63)

Similarly, $\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})$ defined in (29) is
given by

$\frac{k}{\zeta_{1}}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})=\frac{1}{d}U_{1}^{0}\bm{I}+D_{1}^{T}\int\Delta\bm{q}_{1}\nabla_{2}^{T}\rho_{2}d\Gamma^{2}.$
(64)

Because the last term in the preceding equation is finite as $k\to\infty$,
$\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})$ is small and on the order of
$1/k$.

On the other hand, for a strong trap, the relative deviation of the probe
position from the average position is small,
$\Delta|\bm{r}_{1}|/\langle|\bm{r}_{1}|\rangle\ll 1$, which leads to

$\displaystyle\bm{M}_{1}(\bm{r}_{2},\bm{q}_{2},t)$
$\displaystyle=\int\bm{r}_{1}\rho_{2}d\bm{r}_{1}=\bigl{\langle}\bm{r}_{1}\bigr{\rangle}M_{0}+\int\Delta\bm{r}_{1}\rho_{2}d\bm{r}_{1}$
$\displaystyle=\bigl{\langle}\bm{r}_{1}\bigr{\rangle}M_{0}+O(\Delta|\bm{r}_{1}|/\langle|\bm{r}_{1}|\rangle),$
(65)

where the decomposition $\bm{r}_{1}=\langle\bm{r}_{1}\rangle+\Delta\bm{r}_{1}$
is used. Using the first line of (III.5), we have

$\int\nabla_{2}^{T}\rho_{2}\Delta\bm{r}_{1}d\Gamma^{2}=\int\nabla_{2}^{T}\big{[}\bm{M}_{1}-\langle\bm{r}_{1}\rangle
M_{0}\big{]}d\bm{q}_{1}d\Gamma_{2},$ (66)

which is negligible due to the second line of (III.5). Taken together, we
conclude that the last two terms on the rhs of (III.5) are subdominant. To
leading-order, the fluctuation in the strong-trap limit is given by

$k\operatorname{Var}(\bm{r}_{1})=\zeta_{1}D_{1}^{T}\bm{I}=k_{B}T\bm{I},$ (67)

regardless of the presence of the bath particles. Therefore, in this limit we
have

$\displaystyle\bigl{\langle}\left(\Delta\bm{F}_{1}^{e}\right)^{2}\bigr{\rangle}^{1/2}\bigl{\langle}\left(\Delta\bm{r}_{1}\right)^{2}\bigr{\rangle}^{1/2}=dk_{B}T.$
(68)

##  IV Constant-force and constant-velocity microrheology

In this section, we show that the trapped-particle microrheology problem can
be reduced to either the CV or CF problem when appropriate limits are taken.

###  IV.1 Constant-force microrheology

To recover the constant-force microrheology problem, instead of a harmonic
trapping force, we apply a constant force to the probe particle,
$\bm{F}_{1}^{e}=const$, and set the trap velocity
$\bm{U}^{\mathrm{trap}}=\bm{0}$. In this mode of operation, the main quantity
of interest is the average velocity $\langle\bm{U}_{1}\rangle$ of the probe in
response to the constant external driving force. By definition,
$\langle\bm{U}_{1}\rangle=\partial\langle\bm{r}_{1}\rangle/\partial t$, which
can be obtained by considering the rhs of equation (11).

Because the trap is absent, the position $\bm{r}_{1}$ defines an arbitrary
origin in the laboratory frame of reference and the system is statistically
homogeneous (Squires and Brady, 2005). As a result, the conditional
probability $P_{N-1/1}$ defined by

$P_{N}=P_{N-1/1}(\bm{r}^{N-1},\bm{q}^{N-1},t\rvert\bm{r}_{1},\bm{q}_{1},t)P_{1}(\bm{r}_{1},\bm{q}_{1},t)$
(69)

is not a function of $\bm{r}_{1}$. (Note that in general $P_{N-1/1}$ can be a
function of $\bm{q}_{1}$.) The third term on the rhs of equation (11) becomes

$\displaystyle-\bigl{\langle}\bm{D}_{11}^{UF}\cdot\nabla_{1}^{T}\ln
P_{N}\bigr{\rangle}$ $\displaystyle=-\int
d\Gamma^{N-1}P_{N-1/1}\bm{D}_{11}^{UF}\cdot\int
d\Gamma_{1}\nabla_{1}^{T}P_{1}$ $\displaystyle=\bm{0},$ (70)

where we have used the divergence theorem and the fact that $P_{1}$ vanishes
at infinity.

Further manipulations allow us to write equation (11) as

$\displaystyle\langle\bm{U}_{1}\rangle=$ $\displaystyle
U_{1}^{0}\langle\bm{q}_{1}\rangle+\bigl{\langle}\bm{M}_{11}^{UF}\bigr{\rangle}\cdot\bm{F}_{1}^{e}$
$\displaystyle-\sum_{\beta=1}^{N}\int d\Gamma_{1}P_{1}\int
d\Gamma^{N-1}\left(\bm{D}_{1\beta}^{UF}-\bm{D}_{11}^{UF}\right)\cdot\nabla_{\beta}^{T}P_{N-1/1}.$
(71)

If all $N$ particles (including both the probe and the bath particles) are
passive, equation (IV.1) upon integration over $\Gamma_{1}$ reduces to the
average velocity relation originally obtained by Squires and Brady (2005)
(equation A4) for passive colloids.

Neglecting hydrodynamic interactions in the dilute limit, the average velocity
becomes

$\langle\bm{U}_{1}\rangle=\frac{1}{\zeta_{1}}\bm{F}_{1}^{e}+U_{1}^{0}\langle\bm{q}_{1}\rangle+D_{1}^{T}\int\nabla_{2}^{T}\rho_{2}d\Gamma^{2}.$
(72)

Recalling that $M_{0}=\int\rho_{2}d\bm{r}_{1}$, the last term in (72) can be
calculated so long as $M_{0}$ can be obtained. We note that, at long times,
$\langle\bm{q}_{1}\rangle=\bm{0}$. If the probe is under the influence of
external orienting fields, the net polar order $\langle\bm{q}_{1}\rangle$
becomes non-zero (Takatori and Brady, 2014).

In the CF mode of microrheology, the equation governing the spatial moment
$\bm{M}_{m}$ is similar to (32) and can be shown to be

$\displaystyle\frac{\partial\bm{M}_{m}}{\partial
t}-m\left[U_{1}^{0}\bm{q}_{1}\bm{M}_{m-1}+\frac{1}{\zeta_{1}}\bm{F}_{1}^{e}\bm{M}_{m-1}+(m-1)D_{1}^{T}\bm{M}_{m-2}\bm{I}+D_{1}^{T}\nabla_{2}^{T}\bm{M}_{m-1}\right]^{\mathrm{sym}}$
$\displaystyle+\nabla_{2}^{T}\cdot\left(\bm{U}_{r}^{0}\bm{M}_{m}-D_{r}^{T}\nabla_{2}^{T}\bm{M}_{m}-\frac{1}{\zeta_{1}}\bm{F}_{1}^{e}\bm{M}_{m}\right)-mD_{1}^{T}\left[\nabla_{2}^{T}\bm{M}_{m-1}\right]^{\mathrm{sym}}-\sum_{\alpha=1}^{2}D_{\alpha}^{R}\nabla_{\alpha}^{R}\cdot\nabla_{\alpha}^{R}\bm{M}_{m}=\bm{0}.$
(73)

Here, because the external force is constant, the moment equation at rank $m$
only depends on moments of lower ranks and the system up to any rank is a
closed set of equations.

At contact, $r_{2}=R_{c}$, we have the no-flux boundary condition:

$\displaystyle\bm{n}_{2}\cdot\left(\bm{U}_{r}^{0}\bm{M}_{m}-D_{r}^{T}\nabla_{2}^{T}\bm{M}_{m}-\frac{1}{\zeta_{1}}\bm{F}_{1}^{e}\bm{M}_{m}\right)$
$\displaystyle-
mD_{1}^{T}\left[\bm{n}_{2}\bm{M}_{m-1}\right]^{\mathrm{sym}}=\bm{0}.$ (74)

The far-field condition as $r_{2}\to\infty$ is unchanged and given by equation
(36), where $\bm{\Phi}_{m}$ for constant force satisfies

$\displaystyle\frac{\partial\bm{\Phi}_{m}}{\partial
t}-m\left[U_{1}^{0}\bm{q}_{1}\bm{\Phi}_{m-1}+\frac{\bm{F}_{1}^{e}}{\zeta_{1}}\bm{\Phi}_{m-1}+(m-1)D_{1}^{T}\bm{\Phi}_{m-2}\bm{I}\right]^{\mathrm{sym}}$
$\displaystyle-
D_{1}^{R}\nabla_{1}^{R}\cdot\nabla_{1}^{R}\bm{\Phi}_{m}=\bm{0}.$ (75)

To find the average velocity given in equation (72), one needs to consider
equations (73) and (IV.1) for $m=0$.

We note that in the above general formulation, both the probe particle and the
bath particle are ABPs. By setting $U_{1}^{0},U_{2}^{0}=0$ and integrating out
the orientational degrees of freedom of the probe and the bath particle, we
obtain the CF microrheology problem of a passive Brownian probe in a passive
Brownian suspension, which has been considered by Squires and Brady (2005). On
the other hand, the CF microrheology of a passive Brownian probe in an active
Brownian suspension ($U_{1}^{0}=0,U_{2}^{0}\neq 0$) is studied by Burkholder
and Brady (2019).

Taking $m=0$ in equation (73) in the absence of the external force
($\bm{F}_{1}^{e}=\bm{0}$), we obtain

$\displaystyle\frac{\partial M_{0}}{\partial
t}+\nabla_{2}^{T}\cdot\left(\bm{U}_{r}^{0}M_{0}-D_{r}^{T}\nabla_{2}^{T}M_{0}\right)$
$\displaystyle-\sum_{\alpha=1}^{2}D_{\alpha}^{R}\nabla_{\alpha}^{R}\cdot\nabla_{\alpha}^{R}M_{0}=0.$
(76)

Treating the probe as one of the suspension particles, this zeroth spatial
moment is the pair-correlation function of an active Brownian suspension
(subject to proper normalization) in the dilute limit by neglecting all higher
order correlations. Equation (IV.1) governing the pair-correlation at steady
state in 2D has been studied (Poncet  _et al._ , 2021; Dhont, Park, and
Briels, 2021).

###  IV.2 Force-induced tracer diffusion

In the constant-force mode of microrheology, it is also of importance to
consider the force-induced diffusion of the probe particle. In this context,
the probe is often referred to as the tracer, i.e., force-induced tracer
diffusion. If no external force is applied, $\bm{F}_{1}^{e}=\bm{0}$, the
problem is simply called tracer diffusion. The long-time diffusivity of the
tracer in the presence of bath particles can be written as

$\displaystyle\bm{D}_{1}^{\mathrm{eff}}$
$\displaystyle=\lim_{t\to\infty}\frac{1}{2}\frac{d}{dt}\operatorname{Var}(\bm{r}_{1})$
$\displaystyle=\lim_{t\to\infty}\frac{1}{2}\left[\frac{\partial}{\partial
t}\langle\bm{r}_{1}\bm{r}_{1}\rangle-\langle\bm{U}_{1}\rangle\langle\bm{r}_{1}\rangle-\langle\bm{r}_{1}\rangle\langle\bm{U}_{1}\rangle\right],$
(77)

where the covariance tensor of $\bm{r}_{1}$ is governed by

$\displaystyle\frac{1}{2}\frac{d}{dt}\operatorname{Var}(\bm{r}_{1})=$
$\displaystyle
D_{1}^{T}\bm{I}+U_{1}^{0}\big{[}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})\big{]}^{\mathrm{sym}}$
$\displaystyle+D_{1}^{T}\Big{[}\int\Delta\bm{r}_{1}\nabla_{2}^{T}\rho_{2}d\Gamma^{2}\Big{]}^{\mathrm{sym}},$
(78)

and the covariance of $\bm{q}_{1}$ and $\bm{r}_{1}$ satisfies

$\displaystyle\frac{d}{dt}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})+\frac{d-1}{\tau_{1}^{R}}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})=$
$\displaystyle U_{1}^{0}\operatorname{Var}(\bm{q}_{1})$
$\displaystyle+D_{1}^{T}\int\Delta\bm{q}_{1}\nabla_{2}^{T}\rho_{2}d\Gamma^{2}.$
(79)

At long times, we then obtain the diffusivity as

$\displaystyle\bm{D}_{1}^{\mathrm{eff}}=$
$\displaystyle\left(D_{1}^{T}+D_{1}^{\mathrm{swim}}\right)\bm{I}$
$\displaystyle+D_{1}^{T}\left[\int\left(\frac{\ell_{1}}{d-1}\Delta\bm{q}_{1}+\Delta\bm{r}_{1}\right)\nabla_{2}^{T}\rho_{2}d\Gamma^{2}\right]^{\mathrm{sym}}.$
(80)

In (IV.2), the first bracketed term on the rhs is the diffusivity of a single
ABP in free space and the remaining terms are the additional contributions due
to the excluded-volume interaction with the bath particles. As alluded to
earlier, equation (IV.2) is identical to equation (III.4), which is obtained
in the weak-trap limit. We note that in (IV.2) there is a constant external
force while the diffusivity obtained in (III.4) is for a free tracer, i.e.,
force-induced versus free tracer diffusion. It is clear that if the force is
absent the diffusivities obtained from (III.4) and (IV.2) are identical.

Using the divergence theorem, we can relate the integrals on the rhs of (IV.2)
to the zeroth and first spatial moments,

$\displaystyle\int\Delta\bm{q}_{1}\nabla_{2}^{T}\rho_{2}d\Gamma^{2}=\int\Delta\bm{q}_{1}d\bm{q}_{1}d\bm{q}_{2}\oint_{S_{c}}\bm{n}_{2}M_{0}dS_{2},$
(81)
$\displaystyle\int\Delta\bm{r}_{1}\nabla_{2}^{T}\rho_{2}d\Gamma^{2}=\int\Delta\bm{q}_{1}d\bm{q}_{1}d\bm{q}_{2}\oint_{S_{c}}\big{(}\bm{M}_{1}-\langle\bm{r}_{1}\rangle
M_{0}\big{)}\bm{n}_{2}dS_{2},$ (82)

where $\langle\bm{r}_{1}\rangle=\int\bm{M}_{1}/(N-1)d\bm{q}_{1}d\Gamma_{2}$.
Therefore, one only needs to solve for $M_{0}$ and $\bm{M}_{1}$ in equation
(73) in order to calculate the diffusivity.

The above formulation for the forced-induced diffusion of an active tracer in
an active suspension is a direct extension of the generalized Taylor
dispersion theory (GTDT). In particular, we have used the statistical moment
method of Frankel and Brenner (1989). An equivalent approach is to derive the
mean velocity and the diffusivity by first transforming the unbounded
coordinate $\bm{r}_{1}$ into the Fourier space and consider a small wave-
number expansion (Zia and Brady, 2010; Burkholder and Brady, 2017, 2019).

By setting $U_{1}^{0},U_{2}^{0}=0$ and integrating over the orientational
degrees of freedom of both the probe and the bath particles, we recover the
equations governing the force-induced diffusion of a passive probe in a
passive suspension (Zia and Brady, 2010). To recover the problem of a passive
free tracer in an active suspension studied by Burkholder and Brady (2017),
one can set $\bm{F}_{1}^{e}=\bm{0}$, $U_{1}^{0}=0$ and integrate over the
orientational degrees of freedom of the probe.

###  IV.3 Constant-velocity microrheology

To obtain the equations for the CV microrheology problem, we first consider
the probe to have deterministic dynamics with $U_{1}^{0}=0$, $D_{1}^{T}=0$ and
$D_{1}^{R}=0$. Equation (III.2) at steady-state then leads to
$k\langle\bm{r}_{1}\rangle/\zeta_{1}=-\bm{U}^{\mathrm{trap}}$. Furthermore, we
consider the limit of a strong trap in which case the probe tightly follows
the trap velocity. In this limit, the probe velocity is the trap velocity to
leading-order and we then achieve a CV probe.

To see this, we first decompose the position of the probe via
$\bm{r}_{1}=\langle\bm{r}_{1}\rangle+\Delta\bm{r}_{1}$. In the strong trap
limit, the deviation of the probe from the mean position is small,
$\Delta|\bm{r}_{1}|/\langle|\bm{r}_{1}|\rangle\ll 1$.

To leading-order, (III.5) allows us to obtain the first spatial moment as
$k\bm{M}_{1}/\zeta_{1}=-\bm{U}^{\mathrm{trap}}M_{0}$ (this relation can also
be viewed as a closure for the spatial moments). Substitution of this relation
into (III.2) leads to

$\displaystyle\frac{\partial M_{0}}{\partial t}$
$\displaystyle+\nabla_{2}^{T}\cdot\left(U_{2}^{0}\bm{q}_{2}M_{0}-D_{1}^{T}\nabla_{2}^{T}M_{0}-\bm{U}^{\mathrm{trap}}M_{0}\right)$
$\displaystyle-D_{2}^{R}\nabla_{2}^{R}\cdot\nabla_{2}^{R}M_{0}=0.$ (83)

Similarly, the no-flux condition at contact ($\bm{r}_{2}\in S_{c}$) is

$\bm{n}_{2}\cdot\left(U_{2}^{0}\bm{q}_{2}M_{0}-D_{1}^{T}\nabla_{2}^{T}M_{0}-\bm{U}^{\mathrm{trap}}M_{0}\right)=0.$
(84)

Equation (IV.3) describes the distribution of the bath particle measured in a
frame of reference that is co-moving with $\bm{U}^{\mathrm{trap}}$. Realizing
that the probe velocity is the same as the ‘trap’,
$\bm{U}^{\mathrm{probe}}=\bm{U}^{\mathrm{trap}}$, this is the CV microrheology
of an active Brownian suspension. We note that in (IV.3) [cf. (III.2)] the
relative velocity is $U_{2}^{0}\bm{q}_{2}-\bm{U}^{\mathrm{trap}}$ and the
relative diffusivity is $D_{1}^{T}$ because the probe has prescribed
kinematics.

The CV microrheology of an active Brownian suspension governed by (IV.3) and
(84) has been studied by Burkholder and Brady (2020) and Peng and Brady
(2022). To recover the CV microrheology of a passive Brownian suspension
considered by Squires and Brady (2005), one only needs to set $U_{2}^{0}=0$
and integrate over the orientational degrees of freedom of the bath ABP.

##  V Conclusions

In this paper we have considered the trapped-particle microrheology of an
active colloidal suspension consisting of active Brownian spheres. In the
classical models of colloidal microrheology, the applied external force or the
probe velocity are fixed and not subject to random fluctuations. This
constraint of either CF or CV allows a model simpler than that discussed in
the present paper. For the purpose of quantifying the micro-viscous response
of suspensions, the CF or CV models are often sufficient. The challenge arises
if one wishes to consider the fluctuations of the probe as a result of its
interactions with the bath particles and the solvent. More specifically, we
have demonstrated that in order to provide a meaningful quantification of the
fluctuations in the probe position, one must allow both the position of and
the external force on the probe to fluctuate. To achieve this, we developed a
generalized microrheology model in which the probe is driven by a translating
harmonic trap. We explicitly formulated the equations governing the dynamics
of the probe-bath pair in the dilute limit and showed that both the mean
position and the fluctuation of the probe position can be given in terms of
the joint probability distribution.

In the weak-trap limit, we showed that at an intermediate time the probe
exhibits a diffusive behavior in which the diffusivity is the effective
diffusivity of a free tracer immersed in the suspension. At this timescale,
the probe has explored the suspension but hasn’t reached the boundary of the
trap. In other words, it is equivalent to the free-tracer diffusion problem.
For a strong trap, the fluctuations from the activity of the bath particles or
from the collisions between the probe and the bath particles are suppressed
due to the strong confinement of the trap. In this limit, the fluctuation of
the probe originates from the thermal energy alone regardless of the presence
(or activity) of the bath particles.

To conclude, we note that the derived Smoluchowski equation—even at the pair
level—has a high dimensionality, which presents a challenge for the
computation of the probability density. To circumvent this, one can start from
a micromechanical perspective using the Langevin equations and consider a
dynamic simulation of the suspension and the probe (Foss and Brady, 2000;
Carpen and Brady, 2005; Brady and Bossis, 1988). In a dynamic simulation, the
discrete trajectory of the probe is recorded and the calculation of its mean
and fluctuation is straightforward.

###### Acknowledgements.

This work is supported by the National Science Foundation under Grant No. CBET
1803662.

## Data Availability Statement

The data that support the findings of this study are available within the
article.

##  Appendix A Derivation of the pair problem

We integrate equation (III) over the relative positions and the orientations
of the bath particles labeled from $3$ to $N$ to obtain

$\displaystyle\frac{\partial P_{2}}{\partial t}$
$\displaystyle+\nabla_{1}^{T}\cdot\int\left(\bm{j}_{1}^{T}-\bm{U}^{\mathrm{trap}}P_{N}\right)d\Gamma^{N-2}$
$\displaystyle+\nabla_{2}^{T}\cdot\int\left(\bm{j}_{2}^{T}-\bm{j}_{1}^{T}\right)d\Gamma^{N-2}$
$\displaystyle+\sum_{\alpha=1}^{2}\nabla_{\alpha}^{R}\cdot\int\bm{j}_{\alpha}^{R}d\Gamma^{N-2}=0,$
(85)

where $d\Gamma^{N-2}$ is a shorthand for $\prod_{\beta=3}^{N}d\Gamma_{\beta}$
and $P_{2}=\int P_{N}d\Gamma^{N-2}$. In deriving the preceding equation, the
divergence theorem and the no-flux condition are used to eliminate the terms
$\int\nabla_{\beta}^{T}\cdot\left(\bm{j}_{\beta}^{T}-\bm{j}_{1}^{T}\right)d\Gamma^{N-2}$
for $\beta=3,...,N$. In addition, the relation
$\int\nabla_{\alpha}^{R}\cdot\bm{j}_{\alpha}^{R}d\bm{q}_{\alpha}=0$ is used.

To proceed further, we define the conditional probability of finding the
remaining $N-2$ particles, $P_{(N-2)/2}$, given the configuration of the probe
and the _first_ bath particle:

$P_{N}\left(\bm{r}^{N},\bm{q}^{N},t\right)=P_{(N-2)/2}\left(\bm{r}^{N-2},\bm{q}^{N-2},t\big{\rvert}\bm{r}_{2},\bm{q}_{2},\bm{r}_{1},\bm{q}_{1},t\right)P_{2}(\bm{r}_{2},\bm{q}_{2},\bm{r}_{1},\bm{q}_{1},t).$
(86)

Notice that the conditional probability is conserved, $\int
P_{N-2/2}d\Gamma^{N-2}=1$. In equation (85), for $\alpha=1$ or $2$, we have

$\displaystyle\int\bm{j}_{\alpha}^{T}d\Gamma^{N-2}=$ $\displaystyle
U_{\alpha}^{0}\bm{q}_{\alpha}P_{2}+\bigl{\langle}\bm{M}_{\alpha
1}^{UF}\bigr{\rangle}_{(N-2)/2}\cdot\bm{F}_{1}^{e}P_{2}-\bigl{\langle}\bm{D}_{\alpha
2}^{UF}-\bm{D}_{\alpha
1}^{UF}\bigr{\rangle}_{(N-2)/2}\cdot\nabla_{2}^{T}P_{2}$
$\displaystyle-\sum_{\beta=2}^{N}\Bigl{\langle}\left(\bm{D}_{\alpha\beta}^{UF}-\bm{D}_{\alpha
1}^{UF}\right)\cdot\nabla_{\beta}^{T}\ln
P_{(N-2)/2}\Bigr{\rangle}_{(N-2)/2}P_{2}$
$\displaystyle-\Bigl{\langle}\bm{D}_{\alpha 1}^{UF}\cdot\nabla_{1}^{T}\ln
P_{(N-2)/2}\Bigr{\rangle}_{(N-2)/2}P_{2}-\bigl{\langle}\bm{D}_{\alpha
1}^{UF}\bigr{\rangle}_{(N-2)/2}\cdot\nabla_{1}^{T}P_{2}$
$\displaystyle-\sum_{\beta=1}^{2}\Bigl{\langle}\bm{D}_{\alpha\beta}^{UL}\cdot\nabla_{\beta}^{R}\ln
P_{(N-2)/2}\Bigr{\rangle}_{(N-2)/2}P_{2}-\sum_{\beta=1}^{2}\bigl{\langle}\bm{D}_{\alpha\beta}^{UL}\bigr{\rangle}_{(N-2)/2}\cdot\nabla_{\beta}^{R}P_{2},$
(87)

and

$\displaystyle\int\bm{j}_{\alpha}^{R}d\Gamma^{N-2}=$
$\displaystyle\bigl{\langle}\bm{M}_{\alpha 1}^{\Omega
F}\bigr{\rangle}_{(N-2)/2}\cdot\bm{F}_{1}^{e}P_{2}-\Bigl{\langle}\bm{D}_{\alpha
2}^{\Omega F}-\bm{D}_{\alpha 1}^{\Omega
F}\Bigr{\rangle}_{(N-2)/2}\cdot\nabla_{2}^{T}P_{2}$
$\displaystyle-\sum_{\beta=2}^{N}\Bigl{\langle}\left(\bm{D}_{\alpha\beta}^{\Omega
F}-\bm{D}_{\alpha 1}^{\Omega F}\right)\cdot\nabla_{\beta}^{T}\ln
P_{(N-2)/2}\Bigr{\rangle}_{(N-2)/2}P_{2}$
$\displaystyle-\Bigl{\langle}\bm{D}_{\alpha 1}^{\Omega
F}\cdot\nabla_{1}^{T}\ln
P_{(N-2)/2}\Bigr{\rangle}_{(N-2)/2}P_{2}-\bigl{\langle}\bm{D}_{\alpha
1}^{\Omega F}\bigr{\rangle}_{(N-2)/2}\cdot\nabla_{1}^{T}P_{2}$
$\displaystyle-\sum_{\beta=1}^{2}\Bigl{\langle}\bm{D}_{\alpha\beta}^{\Omega
L}\cdot\nabla_{\beta}^{R}\ln
P_{(N-2)/2}\Bigr{\rangle}_{(N-2)/2}P_{2}-\sum_{\beta=1}^{2}\Bigl{\langle}\bm{D}_{\alpha\beta}^{\Omega
L}\Bigr{\rangle}_{(N-2)/2}\cdot\nabla_{\beta}^{R}P_{2}-D_{\alpha}^{R}\nabla_{\alpha}^{R}P_{2}.$
(88)

In equations (A) and (A), we have defined
$\langle(\cdot)\rangle_{(N-2)/2}=\int(\cdot)P_{(N-2)/2}d\Gamma^{N-2}$, and
used the fact that the mobility tensors are independent of $\bm{q}^{N}$ for
spheres, i.e.,
$\bm{M}_{\alpha\beta}=\bm{M}_{\alpha\beta}(\bm{r}_{2},...,\bm{r}_{N})$.

In the dilute limit, neglecting the terms involving the gradients of $\ln
P_{(N-2)/2}$ and using the pair mobility tensor in the absence of other
particles in place of $\langle\bm{M}\rangle_{(N-2)/2}$, we obtain

$\displaystyle\frac{\partial P_{2}}{\partial t}$
$\displaystyle+\nabla_{1}^{T}\cdot\left(\bm{j}_{1}^{T}-\bm{U}^{\mathrm{trap}}P_{2}\right)$
$\displaystyle+\nabla_{2}^{T}\cdot\left(\bm{j}_{2}^{T}-\bm{j}_{1}^{T}\right)+\sum_{\alpha=1}^{2}\nabla_{\alpha}^{R}\cdot\bm{j}_{\alpha}^{R}=0,$
(89)

where using the same symbols as before

$\displaystyle\bm{j}_{\alpha}^{T}$ $\displaystyle=$ $\displaystyle
U_{\alpha}^{0}\bm{q}_{\alpha}P_{2}+\bm{M}_{\alpha
1}^{UF}\cdot\bm{F}_{1}^{e}P_{2}-\left(\bm{D}_{\alpha 2}^{UF}-\bm{D}_{\alpha
1}^{UF}\right)\cdot\nabla_{2}^{T}P_{2}$ (90) $\displaystyle-\bm{D}_{\alpha
1}^{UF}\cdot\nabla_{1}^{T}P_{2}-\sum_{\beta=1}^{2}\bm{D}_{\alpha\beta}^{UL}\cdot\nabla_{\beta}^{R}P_{2}$
$\displaystyle\bm{j}_{\alpha}^{R}$ $\displaystyle=$
$\displaystyle\bm{M}_{\alpha 1}^{\Omega
F}\cdot\bm{F}_{1}^{e}P_{2}-\left(\bm{D}_{\alpha 2}^{\Omega F}-\bm{D}_{\alpha
1}^{\Omega F}\right)\cdot\nabla_{2}^{T}P_{2}-\bm{D}_{\alpha 1}^{\Omega
F}\cdot\nabla_{1}^{T}P_{2}$ (91)
$\displaystyle-\sum_{\beta=1}^{2}\bm{D}_{\alpha\beta}^{\Omega
L}\cdot\nabla_{\beta}^{R}P_{2}-D_{\alpha}^{R}\nabla_{\alpha}^{R}P_{2}.$

In the absence of hydrodynamic interactions, we have
$\bm{M}_{\alpha\beta}^{UF}=\bm{I}\delta_{\alpha\beta}/\zeta_{\alpha}^{T}$,
$\bm{M}_{\alpha\beta}^{\Omega
L}=\bm{I}\delta_{\alpha\beta}/\zeta_{\alpha}^{R}$, and
$\bm{M}_{\alpha\beta}^{UL},\bm{M}_{\alpha\beta}^{\Omega F}=\bm{0}$, where
$\delta_{\alpha\beta}$ is the Kronecker delta. The conditional probability of
finding a bath particle,
$\rho_{1/1}(\bm{r}_{2},\bm{q}_{2},t\rvert\bm{r}_{1},\bm{q}_{1},t)$, can be
related to $P_{2}$ via the relation $\rho_{1/1}=(N-1)P_{1/1}$, where $P_{1/1}$
is defined by $P_{2}=P_{1/1}P_{1}$. The factor of $N-1$ comes from the process
of removing the “labels” of the $N-1$ bath particles. From this, the joint
probability density of finding a bath particle at $\bm{r}_{2}$, $\bm{q}_{2}$
and the probe at $\bm{r}_{1}$, $\bm{q}_{1}$ is defined as
$\rho_{2}=\rho_{1/1}P_{1}$. Furthermore, we can define a dimensionless
conditional distribution function $g_{1/1}$ such that

$\rho_{2}=\rho_{1/1}P_{1}=n_{b}g_{1/1}P_{1},$ (92)

where $n_{b}=(N-1)/V$ is the number density of bath particles.

In the absence of hydrodynamic interactions, equations (89)-(91) reduce to
equations (III.2)-(21) given in the text.

##  Appendix B Derivation of the variance relations

For the pair problem, equation (12) governing the mean-squared displacement
reduces to

$\displaystyle\frac{1}{2}\frac{\partial\langle\bm{r}_{1}\bm{r}_{1}\rangle}{\partial
t}+\frac{1}{\tau_{k}}\bigl{\langle}\bm{r}_{1}\bm{r}_{1}\bigr{\rangle}=$
$\displaystyle
D_{1}^{T}\bm{I}+\left[U_{1}^{0}\bigl{\langle}\bm{q}_{1}\bm{r}_{1}\bigr{\rangle}-\bm{U}^{\mathrm{trap}}\langle\bm{r}_{1}\rangle\right]^{\mathrm{sym}}$
$\displaystyle+D_{1}^{T}\Big{[}\int\nabla_{2}^{T}\rho_{2}\bm{r}_{1}d\Gamma^{2}\Big{]}^{\mathrm{sym}}.$
(93)

Using equations (III.2) and (B), one can show that the position fluctuation of
the probe is governed by

$\displaystyle\frac{1}{2}\frac{\partial\operatorname{Var}(\bm{r}_{1})}{\partial
t}+\frac{1}{\tau_{k}}\operatorname{Var}(\bm{r}_{1})=$ $\displaystyle
D_{1}^{T}\bm{I}+U_{1}^{0}\left[\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})\right]^{\mathrm{sym}}$
$\displaystyle+D_{1}^{T}\Big{[}\int\nabla_{2}^{T}\rho_{2}\Delta\bm{r}_{1}d\Gamma^{2}\Big{]}^{\mathrm{sym}},$
(94)

where
$\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})=\langle\bm{q}_{1}\bm{r}_{1}\rangle-\langle\bm{q}_{1}\rangle\langle\bm{r}_{1}\rangle$
and recall that $\Delta\bm{r}_{1}=\bm{r}_{1}-\langle\bm{r}_{1}\rangle$. To
calculate the covariance of $\bm{q}_{1}$ and $\bm{r}_{1}$ appearing in
equation (B), we need $\langle\bm{q}_{1}\bm{r}_{1}\rangle$,
$\langle\bm{q}_{1}\rangle$ and $\langle\bm{r}_{1}\rangle$.

The net polar order of the probe satisfies

$\frac{\partial\langle\bm{q}_{1}\rangle}{\partial
t}+\frac{d-1}{\tau_{1}^{R}}\langle\bm{q}_{1}\rangle=\bm{0},$ (95)

where $d(=2,3)$ is the dimensionality of the problem. It can be seen that the
net polar order of the probe is not affected by the trap or the bath
particles. The full solution to (95) is readily obtained as

$\langle\bm{q}_{1}\rangle(t)=\exp\left[-(d-1)t/\tau_{1}^{R}\right]\langle\bm{q}_{1}\rangle(0),$
(96)

where any initial net polar order $\langle\bm{q}_{1}\rangle(0)$ decays away
exponentially due to the rotary diffusion.

The average of $\bm{q}_{1}\bm{r}_{1}$ is governed by

$\displaystyle\frac{\partial\bigl{\langle}\bm{q}_{1}\bm{r}_{1}\bigr{\rangle}}{\partial
t}+\frac{1}{\tau}\bigl{\langle}\bm{q}_{1}\bm{r}_{1}\bigr{\rangle}=$
$\displaystyle-\langle\bm{q}_{1}\rangle\bm{U}^{\mathrm{trap}}+U_{1}^{0}\bigl{\langle}\bm{q}_{1}\bm{q}_{1}\bigr{\rangle}$
$\displaystyle+D_{1}^{T}\int\bm{q}_{1}\nabla_{2}^{T}\rho_{2}d\Gamma^{2},$ (97)

where $\langle\bm{q}_{1}\rangle$ is given by (95) and
$\bigl{\langle}\bm{q}_{1}\bm{q}_{1}\bigr{\rangle}$ satisfies

$\frac{\partial\bigl{\langle}\bm{q}_{1}\bm{q}_{1}\bigr{\rangle}}{\partial
t}+\frac{2d}{\tau_{1}^{R}}\left[\bigl{\langle}\bm{q}_{1}\bm{q}_{1}\bigr{\rangle}-\frac{1}{d}\bm{I}\right]=\bm{0}.$
(98)

Similarly to $\langle\bm{q}_{1}\rangle$, the net nematic order of the probe
regardless of the presence of the trap or the bath particles is given by

$\langle\bm{Q}_{1}\rangle(t)=\exp\left[-2dt/\tau_{1}^{R}\right]\bm{Q}_{1}(0),$
(99)

where we have defined the net trace-free nematic tensor
$\langle\bm{Q}_{1}\rangle=\bigl{\langle}\bm{q}_{1}\bm{q}_{1}\bigr{\rangle}-\bm{I}/d$.

At long times ($t\to\infty$), there is no net polar order of the probe,
$\langle\bm{q}_{1}\rangle=\bm{0}$ and the net nematic order is isotropic,
$\bigl{\langle}\bm{q}_{1}\bm{q}_{1}\bigr{\rangle}=\bm{I}/d$.

Using equations (III.2), (95) and (B), we obtain

$\displaystyle\frac{\partial\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})}{\partial
t}+\frac{1}{\tau}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})=$ $\displaystyle
U_{1}^{0}\operatorname{Cov}(\bm{q}_{1},\bm{q}_{1})$
$\displaystyle+D_{1}^{T}\int\Delta\bm{q}_{1}\nabla_{2}^{T}\rho_{2}d\Gamma^{2},$
(100)

where $\Delta\bm{q}_{1}=\bm{q}_{1}-\langle\bm{q}_{1}\rangle$.

##  Appendix C Asymptotic analysis of the probe in the absence of bath
particles

In equation (46), the timescale of transient decay $\tau$ can be written as

$\frac{1}{\tau}=\frac{d-1}{\tau_{1}^{R}}+\frac{1}{\tau_{k}}=\frac{d-1+\epsilon}{\tau_{1}^{R}},$
(101)

where $\epsilon=\tau_{1}^{R}/\tau_{k}=k\tau_{1}^{R}/\zeta_{1}$. Using this
definition, the solution of (46) is given by

$\displaystyle\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})(t)=$ $\displaystyle
e^{-t/\tau}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})(0)$
$\displaystyle+U_{1}^{0}\int_{0}^{t}\exp\left(-\frac{t-s}{\tau}\right)\operatorname{Var}(\bm{q}_{1})(s)ds.$
(102)

From equation (99), the preceding equation becomes

$\displaystyle\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})(t)=$ $\displaystyle
e^{-t/\tau}\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})(0)+\frac{U_{1}^{0}\tau\bm{I}}{d}\left(1-e^{-t/\tau}\right)$
$\displaystyle-\frac{\ell_{1}}{d+1-\epsilon}\left(e^{-2dt/\tau_{R}}-e^{-t/\tau}\right)\langle\bm{Q}_{1}\rangle(0).$
(103)

In the long-time limit ($t/\tau_{R}\gg 1$ and $t/\tau\gg 1$), we obtain
equation (48) in the text. Using equation (B) in the absence of bath
particles, we obtain

$\displaystyle\operatorname{Var}(\bm{r}_{1})(t)=$ $\displaystyle
e^{-2t/\tau_{k}}\operatorname{Var}(\bm{r}_{1})(0)$
$\displaystyle+\tau_{k}\left(D_{1}^{T}+\frac{\left(U_{1}^{0}\right)^{2}\tau}{d}\right)\bm{I}\left(1-e^{-2t/\tau_{k}}\right)$
$\displaystyle+2U_{1}^{0}\int_{0}^{t}\exp\left[-2\frac{t-s}{\tau_{k}}\right]\left[\operatorname{Cov}^{\prime}(\bm{q}_{1},\bm{r}_{1})(s)\right]^{\mathrm{sym}}ds,$
(104)

where
$\operatorname{Cov}^{\prime}(\bm{q}_{1},\bm{r}_{1})(s)=\operatorname{Cov}(\bm{q}_{1},\bm{r}_{1})(s)-U_{1}^{0}\tau\bm{I}/d$
is the time-dependent (transient) part of the covariance of $\bm{q}_{1}$ and
$\bm{r}_{1}$. The integral in (C) can be carried out explicitly but is not
important for the following discussion.

In the presence of the harmonic trap, the system exhibit two timescales that
are important: the reorientation time $\tau_{1}^{R}$ and the viscoelastic
timescale $\tau_{k}$; their relative importance is characterized by the
parameter $\epsilon$. In the weak-trap limit, $\epsilon\to 0$, the two
timescales are well-separated. It is useful to define the fast time variable
$t_{1}=t$ and the slow time variable $t_{2}=\epsilon t$. We now consider the
limit $\epsilon\to 0$ and the intermediate timescale in which the ABP has
experienced many reorientations due to rotary diffusion but hasn’t reached the
“boundary” of the trap, i.e., $t_{1}/\tau_{1}^{R}\gg 1$ but
$t_{2}/\tau_{1}^{R}=\epsilon t/\tau_{R}\ll 1$.

Differentiating equation (C) leads to

$\displaystyle\frac{d}{dt}\operatorname{Var}(\bm{r}_{1})(t)=$
$\displaystyle\frac{-2}{\tau_{k}}e^{-2t/\tau_{k}}\operatorname{Var}(\bm{r}_{1})(0)+2\left(D_{1}^{T}+\frac{\left(U_{1}^{0}\right)^{2}\tau}{d}\right)\bm{I}e^{-2t/\tau_{k}}$
$\displaystyle+2U_{1}^{0}\left[\operatorname{Cov}^{\prime}(\bm{q}_{1},\bm{r}_{1})(t)\right]^{\mathrm{sym}}+2U_{1}^{0}\frac{-2}{\tau_{k}}\int_{0}^{t}\exp\left[-2\frac{t-s}{\tau_{k}}\right]\left[\operatorname{Cov}^{\prime}(\bm{q}_{1},\bm{r}_{1})(s)\right]^{\mathrm{sym}}ds.$
(105)

Since $1/\tau_{k}=\epsilon/\tau_{1}^{R}$ and

$\frac{t}{\tau_{k}}=\frac{t_{2}}{\epsilon\tau_{k}}=\frac{t_{2}}{\tau_{1}^{R}}\ll
1,\quad\frac{t}{\tau}=\frac{t_{1}}{\tau_{1}^{R}}(d+1-\epsilon)\gg 1,$ (106)

we have

$\displaystyle
e^{-2t/\tau_{k}}=e^{-2t_{2}/\tau_{1}^{R}}=1+O(t_{2}/\tau_{1}^{R}),$ (107)
$\displaystyle\frac{\left(U_{1}^{0}\right)^{2}\tau}{d}=\frac{U_{1}^{0}\ell_{1}}{d(d-1)}\left[1+O(\epsilon)\right].$
(108)

Therefore, equation (C) at leading order is

$\displaystyle\frac{1}{2}\frac{d}{dt}\operatorname{Var}(\bm{r}_{1})(t)$
$\displaystyle=\left(D_{1}^{T}+\frac{U_{1}^{0}\ell_{1}}{d(d-1)}\right)\bm{I}$
$\displaystyle=\left(D_{1}^{T}+D_{1}^{\mathrm{swim}}\right)\bm{I}.$ (109)

It is clear that in the weak-trap limit in this intermediate timescale, the
ABP exhibits a diffusive behavior with the free-space diffusivity
$D_{1}^{T}+D_{1}^{\mathrm{swim}}$.

We now consider the weak-trap limit but at long-times, $t/\tau_{1}^{R}\gg
1,t/\tau_{k}\gg 1$. So long as the trap strength is not identically zero, the
ABP will eventually ($t/\tau_{k}\gg 1$) experience the confinement of the
trap. Using equation (C), we obtain at long times

$\frac{1}{\tau_{k}}\operatorname{Cov}(\bm{r}_{1},\bm{r}_{1})\to\left(D_{1}^{T}+D_{1}^{\mathrm{swim}}\right)\bm{I}.$
(110)

In the strong-trap limit ($\epsilon\to\infty$) and at long times, we have
$\tau/\tau_{1}^{R}=O(1/\epsilon)$ and the position fluctuation of the probe
$\operatorname{Var}(\bm{r}_{1})/\tau_{k}=D_{1}^{T}$.

## References

  * Weihs, Mason, and Teitell (2006) D. Weihs, T. G. Mason, and M. A. Teitell, “Bio-microrheology: A frontier in microrheology,” Biophys. J. 91, 4296–4305 (2006).
  * Cicuta and Donald (2007) P. Cicuta and A. M. Donald, “Microrheology: a review of the method and applications,” Soft Matter 3, 1449–1455 (2007).
  * Wirtz  _et al._ (2009) D. Wirtz _et al._ , “Particle-tracking microrheology of living cells: principles and applications,” Annu. Rev. Biophys. 38, 301–326 (2009).
  * Furst and Squires (2017) E. M. Furst and T. M. Squires,  _Microrheology_ (Oxford University Press, 2017).
  * Wilhelm, Gazeau, and Bacri (2003) C. Wilhelm, F. Gazeau, and J.-C. Bacri, “Rotational magnetic endosome microrheology: Viscoelastic architecture inside living cells,” Phys. Rev. E 67, 061908 (2003).
  * Nawaz  _et al._ (2012) S. Nawaz, P. Sánchez, K. Bodensiek, S. Li, M. Simons, and I. A. T. Schaap, “Cell visco-elasticity measured with afm and optical trapping at sub-micrometer deformations,” PLOS ONE 7, 1–9 (2012).
  * Berret (2016) J.-F. Berret, “Local viscoelasticity of living cells measured by rotational magnetic spectroscopy,” Nat. Commun. 7, 1–9 (2016).
  * Ayala  _et al._ (2016) Y. A. Ayala, B. Pontes, D. S. Ether, L. B. Pires, G. R. Araujo, S. Frases, L. F. Romão, M. Farina, V. Moura-Neto, N. B. Viana, _et al._ , “Rheological properties of cells measured by optical tweezers,” BMC Biophys. 9, 1–11 (2016).
  * Hu  _et al._ (2017) J. Hu, S. Jafari, Y. Han, A. J. Grodzinsky, S. Cai, and M. Guo, “Size- and speed-dependent mechanical behavior in living mammalian cytoplasm,” Proc. Natl. Acad. Sci. USA 114, 9529–9534 (2017).
  * Squires and Brady (2005) T. M. Squires and J. F. Brady, “A simple paradigm for active and nonlinear microrheology,” Phys. Fluids 17, 073101 (2005).
  * Khair and Brady (2005) A. S. Khair and J. F. Brady, ““Microviscoelasticity” of colloidal dispersions,” J. Rheol. 49, 1449–1481 (2005).
  * Khair and Brady (2006) A. S. Khair and J. F. Brady, “Single particle motion in colloidal dispersions: a simple model for active and nonlinear microrheology,” J. Fluid Mech. 557, 73–117 (2006).
  * Meyer  _et al._ (2006) A. Meyer, A. Marshall, B. G. Bush, and E. M. Furst, “Laser tweezer microrheology of a colloidal suspension,” J. Rheol. 50, 77–92 (2006).
  * Zia and Brady (2010) R. N. Zia and J. F. Brady, “Single-particle motion in colloids: force-induced diffusion,” J. Fluid Mech. 658, 188–210 (2010).
  * Swan and Zia (2013) J. W. Swan and R. N. Zia, “Active microrheology: Fixed-velocity versus fixed-force,” Phys. Fluids 25, 083303 (2013).
  * Zia (2018) R. N. Zia, “Active and passive microrheology: Theory and simulation,” Annu. Rev. Fluid Mech. 50, 371–405 (2018).
  * Burkholder and Brady (2019) E. W. Burkholder and J. F. Brady, “Fluctuation-dissipation in active matter,” J. Chem. Phys. 150, 184901 (2019).
  * Burkholder and Brady (2020) E. W. Burkholder and J. F. Brady, “Nonlinear microrheology of active Brownian suspensions,” Soft Matter 16, 1034–1046 (2020).
  * Yao  _et al._ (2009) A. Yao, M. Tassieri, M. Padgett, and J. Cooper, “Microrheology with optical tweezers,” Lab Chip 9, 2568–2575 (2009).
  * Wang (2017) Z.-G. Wang, “50th anniversary perspective: Polymer conformation—a pedagogical review,” Macromolecules 50, 9073–9114 (2017).
  * Jepson  _et al._ (2013) A. Jepson, V. A. Martinez, J. Schwarz-Linek, A. Morozov, and W. C. K. Poon, “Enhanced diffusion of nonswimmers in a three-dimensional bath of motile bacteria,” Phys. Rev. E 88, 041002 (2013).
  * Miño  _et al._ (2013) G. L. Miño, J. Dunstan, A. Rousselet, E. Clément, and R. Soto, “Induced diffusion of tracers in a bacterial suspension: theory and experiments,” J. Fluid Mech. 729, 423–444 (2013).
  * Morozov and Marenduzzo (2014) A. Morozov and D. Marenduzzo, “Enhanced diffusion of tracer particles in dilute bacterial suspensions,” Soft Matter 10, 2748–2758 (2014).
  * Kasyap, Koch, and Wu (2014) T. V. Kasyap, D. L. Koch, and M. Wu, “Hydrodynamic tracer diffusion in suspensions of swimming bacteria,” Phys. Fluids 26, 081901 (2014).
  * Reichhardt and Reichhardt (2015) C. Reichhardt and C. J. O. Reichhardt, “Active microrheology in active matter systems: Mobility, intermittency, and avalanches,” Phys. Rev. E 91, 032313 (2015).
  * Burkholder and Brady (2017) E. W. Burkholder and J. F. Brady, “Tracer diffusion in active suspensions,” Phys. Rev. E 95, 052605 (2017).
  * Knežević, Avilés Podgurski, and Stark (2021) M. Knežević, L. E. Avilés Podgurski, and H. Stark, “Oscillatory active microrheology of active suspensions,” Sci. Rep. 11, 1–10 (2021).
  * Peng and Brady (2022) Z. Peng and J. F. Brady, “Forced microrheology of active colloids,” arXiv preprint arXiv:2205.03715 (2022).
  * Seyforth  _et al._ (2022) H. Seyforth, M. Gomez, W. B. Rogers, J. L. Ross, and W. W. Ahmed, “Nonequilibrium fluctuations and nonlinear response of an active bath,” Phys. Rev. Research 4, 023043 (2022).
  * Burkholder and Brady (2018) E. W. Burkholder and J. F. Brady, “Do hydrodynamic interactions affect the swim pressure?” Soft Matter 14, 3581–3589 (2018).
  * Takatori  _et al._ (2016) S. C. Takatori, R. De Dier, J. Vermant, and J. F. Brady, “Acoustic trapping of active matter,” Nat. Commun. 7, 1–7 (2016).
  * Takatori, Yan, and Brady (2014) S. C. Takatori, W. Yan, and J. F. Brady, “Swim pressure: Stress generation in active matter,” Phys. Rev. Lett. 113, 028103 (2014).
  * Palacci  _et al._ (2010) J. Palacci, C. Cottin-Bizonne, C. Ybert, and L. Bocquet, “Sedimentation and effective temperature of active colloidal suspensions,” Phys. Rev. Lett. 105, 088304 (2010).
  * Takatori and Brady (2014) S. C. Takatori and J. F. Brady, “Swim stress, motion, and deformation of active matter: effect of an external field,” Soft Matter 10, 9433–9445 (2014).
  * Poncet  _et al._ (2021) A. Poncet, O. Bénichou, V. Démery, and D. Nishiguchi, “Pair correlation of dilute active Brownian particles: From low-activity dipolar correction to high-activity algebraic depletion wings,” Phys. Rev. E 103, 012605 (2021).
  * Dhont, Park, and Briels (2021) J. K. G. Dhont, G. W. Park, and W. J. Briels, “Motility-induced inter-particle correlations and dynamics: a microscopic approach for active Brownian particles,” Soft Matter 17, 5613–5632 (2021).
  * Frankel and Brenner (1989) I. Frankel and H. Brenner, “On the foundations of generalized Taylor dispersion theory,” J. Fluid Mech. 204, 97–119 (1989).
  * Foss and Brady (2000) D. R. Foss and J. F. Brady, “Brownian dynamics simulation of hard-sphere colloidal dispersions,” J. Rheol. 44, 629–651 (2000).
  * Carpen and Brady (2005) I. C. Carpen and J. F. Brady, “Microrheology of colloidal dispersions by Brownian dynamics simulations,” J. Rheol. 49, 1483–1502 (2005).
  * Brady and Bossis (1988) J. F. Brady and G. Bossis, “Stokesian dynamics,” Annu. Rev. Fluid Mech. 20, 111–157 (1988).