1. Introduction
Phoretic phenomena, i.e. the movement of particles in response to an external field (Anderson Reference Anderson1989; Velegol et al. Reference Velegol, Garg, Guha, Kar and Kumar2016; Khair Reference Khair2022), is pivotal for a range of applications such as separation of biomacromolecules (Heller Reference Heller2001; Lee et al. Reference Lee, Costumbrado, Hsu and Kim2012), purification of nucleic acids from whole blood (Persat, Marshall & Santiago Reference Persat, Marshall and Santiago2009), measurement of zeta potential (Doane et al. Reference Doane, Chuang, Hill and Burda2012; Shin et al. Reference Shin, Ault, Feng, Warren and Stone2017a), banding of colloidal particles (Abécassis et al. Reference Abécassis, Cottin-Bizonne, Ybert, Ajdari and Bocquet2008; Banerjee & Squires Reference Banerjee and Squires2019; Raj, Shields & Gupta Reference Raj, Shields and Gupta2023b), membraneless water filtration (Shin et al. Reference Shin, Shardt, Warren and Stone2017b) and understanding biological pattern formation (Alessio & Gupta Reference Alessio and Gupta2023), among others. In contrast, self-phoretic particles, also known as microswimmers, respond to a self-generated field gradient (Paxton et al. Reference Paxton, Kistler, Olmeda, Sen, St. Angelo, Cao, Mallouk, Lammert and Crespi2004; Howse et al. Reference Howse, Jones, Ryan, Gough, Vafabakhsh and Golestanian2007; Ebbens & Howse Reference Ebbens and Howse2010; Popescu, Uspal & Dietrich Reference Popescu, Uspal and Dietrich2016; Ganguly, Alessio & Gupta Reference Ganguly, Alessio and Gupta2023). Microswimmers are extensively studied for applications in targeted drug delivery (Xuan et al. Reference Xuan, Shao, Lin, Dai and He2014; Luo et al. Reference Luo, Feng, Wang and Guan2018), environmental remediation (Gao et al. Reference Gao, Feng, Pei, Gu, Li and Wang2013; Wang et al. Reference Wang, Kaeppler, Fischer and Simmchen2019), remote sensing of toxic chemicals (Esteban-Fernández de Ávila et al. Reference Esteban-Fernández de Ávila, Loper-Ramirez, Báez, Jodra, Singh, Kaufman and Wang2016), the autonomous motion of microbots (Zarei & Zarei Reference Zarei and Zarei2018; Hu, Liu & Sun Reference Hu, Liu and Sun2022) and collective behaviour of active colloids (Palacci et al. Reference Palacci, Sacanna, Steinberg, Pine and Chaikin2013; Takatori & Brady Reference Takatori and Brady2016; Illien, Golestanian & Sen Reference Illien, Golestanian and Sen2017). The mobility expressions of phoretic and self-phoretic processes are identical, with the key distinction being that the origin of field gradients in the two processes is different. In phoretic processes, this field gradient is externally imposed on colloidal particles, while in self-phoretic particles they are locally generated by the particles themselves typically through surface reactions or other mechanisms.
Studies on external electrophoretic motion have focused on the dependence of electrophoretic mobility on the effect of particle shape (Yoon & Kim Reference Yoon and Kim1989; Solomentsev & Anderson Reference Solomentsev and Anderson1994), surface heterogeneity (Fair & Anderson Reference Fair and Anderson1992; Velegol, Anderson & Garoff Reference Velegol, Anderson and Garoff1996), finite double-layer thickness (Henry Reference Henry1931; O'Brien & White Reference O'Brien and White1978) and, more recently, strong deformation of double-layers (Khair Reference Khair2018, Reference Khair2022) and charge reversal (Kubíčková et al. Reference Kubíčková, Křížek, Coufal, Pavel, Vazdar, Wernersson, Heyda and Jungwirth2012; Gupta et al. Reference Gupta, Rajan, Carter and Stone2020a). Similarly, researchers have predicted the dependence of diffusiophoretic mobility (Anderson Reference Anderson1989; Brady Reference Brady2011) on finite double-layer thickness (Prieve et al. Reference Prieve, Anderson, Ebel and Lowell1984; Keh & Wei Reference Keh and Wei2000), surface chemistry (Gupta, Shim & Stone Reference Gupta, Shim and Stone2020b) and multiple electrolytes (Gupta, Rallabandi & Stone Reference Gupta, Rallabandi and Stone2019; Alessio et al. Reference Alessio, Shim, Mintah, Gupta and Stone2021). Studies on self-phoretic systems (Ramaswamy Reference Ramaswamy2010; Moran & Posner Reference Moran and Posner2017) focus on the impact of particle shape (Shklyaev, Brady & Córdova-Figueroa Reference Shklyaev, Brady and Córdova-Figueroa2014; Nourhani & Lammert Reference Nourhani and Lammert2016; Poehnl, Popescu & Uspal Reference Poehnl, Popescu and Uspal2020; Daddi-Moussa-Ider et al. Reference Daddi-Moussa-Ider, Nasouri, Vilfan and Golestanian2021; Ganguly & Gupta Reference Ganguly and Gupta2023; Lee et al. Reference Lee, Thome, Cruze, Ganguly, Gupta and Shields2023; Raj et al. Reference Raj, Ganguly, Becker, Shields and Gupta2023a), active patch shape (Lisicki, Reigh & Lauga Reference Lisicki, Reigh and Lauga2018; Lee et al. Reference Lee, Al Harraq, Bishop and Bharti2021), surface interaction (Sharifi-Mood, Koplik & Maldarelli Reference Sharifi-Mood, Koplik and Maldarelli2013) and finite Péclet number (Michelin & Lauga Reference Michelin and Lauga2014).
Broadly speaking, there are two approaches for predicting the mobilities described above. The first approach solves the coupled solute conservation equations and the modified Stokes equation (Henry Reference Henry1931; O'Brien & White Reference O'Brien and White1978; Prieve et al. Reference Prieve, Anderson, Ebel and Lowell1984; Prieve & Roman Reference Prieve and Roman1987; Anderson Reference Anderson1989; Keh & Wei Reference Keh and Wei2000; Sharifi-Mood et al. Reference Sharifi-Mood, Koplik and Maldarelli2013; Khair Reference Khair2018, Reference Khair2022; Gupta et al. Reference Gupta, Rallabandi and Stone2019) and employs a force-free and torque-free condition to arrive at the translation velocity, 
$\boldsymbol {U}$, and rotational velocity, 
$\boldsymbol {\varOmega }$, of the particle. Thus, the above approach requires resolving the interaction potential simultaneously with the hydrodynamic equations. While exact and powerful, the methodology described above is cumbersome for analytical results when the particle–fluid interaction potential decays at much larger length scales than the particle size. Further, the solution strategy needs to be revised whenever the interaction potential changes, making it less convenient to be integrated into other analyses.
The second approach employs the reciprocal theorem in the thin interaction length limit (Stone & Samuel Reference Stone and Samuel1996; Brady Reference Brady2011; Michelin & Lauga Reference Michelin and Lauga2014; Lisicki et al. Reference Lisicki, Reigh and Lauga2018; Masoud & Stone Reference Masoud and Stone2019; Poehnl et al. Reference Poehnl, Popescu and Uspal2020; Ganguly & Gupta Reference Ganguly and Gupta2023; Raj et al. Reference Raj, Ganguly, Becker, Shields and Gupta2023a). In this limit, it is assumed that there exists a slip velocity, 
$\boldsymbol {u}_s$, at the interface of the inner region where the interaction potential is non-zero, and the outer region where the interaction potential is zero. This allows one to treat the outer problem as a classical Stokes flow problem with a slip boundary condition. Consequently, 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ can be represented as surface integrals of appropriate functions of 
$\boldsymbol {u}_s$. This approach based on the reciprocal theorem was first utilized by Stone & Samuel (Reference Stone and Samuel1996) to study the impact of distortions in spherical microswimmers. This methodology is particularly powerful because, unlike the first method, computing 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ is relatively straightforward and agnostic to the mechanistic origin of 
$\boldsymbol {u}_s$. However, this approach is valid only when the interaction potential length is significantly smaller than the particle size, restricting its applicability. Additionally, it requires the knowledge of 
$\boldsymbol {u}_s$ a priori, and most studies have to rely on a lumped mobility parameter to estimate the value of 
$\boldsymbol {u}_s$ (Michelin & Lauga Reference Michelin and Lauga2014; Lisicki et al. Reference Lisicki, Reigh and Lauga2018; Poehnl & Uspal Reference Poehnl and Uspal2021; Ganguly & Gupta Reference Ganguly and Gupta2023; Raj et al. Reference Raj, Ganguly, Becker, Shields and Gupta2023a).
In this work, we seek to unify the benefits realized through the two approaches by employing the results of Brady (Reference Brady2021) and demonstrating how it can reconcile a large volume of mobility results for externally driven and self-phoretic propulsion of particles, and using these results for additional analyses. This approach is similar to the prior literature on the inertial correction to Stokes flow (Brenner & Cox Reference Brenner and Cox1963; Hinch Reference Hinch1991; Leal Reference Leal2007), and swimming through non-Newtonian fluids (Datt et al. Reference Datt, Zhu, Elfring and Pak2015; Elfring & Goyal Reference Elfring and Goyal2016; Datt et al. Reference Datt, Natale, Hatzikiriakos and Elfring2017), where the first-order corrections include a body force term from the leading order and reciprocal theorem is employed to find the resulting motion. In § 2, we obtain a general mobility expression for an arbitrary particle shape, subjected to an osmophoretic (combination of osmotic and phoretic force) body force 
$\boldsymbol {b}$, identical to the results in Brady (Reference Brady2021). Subsequently, we take our expression to the thin interaction length scale limit and retrieve the mobility expressions in Stone & Samuel (Reference Stone and Samuel1996), see § 3. In § 3, we also retrieve the expression for the electrophoretic mobility of translation of a charged spherical particle in an external electric field, obtained by Henry (Reference Henry1931). Additionally, we derive the diffusiophoretic mobility of a charged spherical particle in an externally imposed solute gradient at finite interaction lengths, as first obtained by Keh & Wei (Reference Keh and Wei2000). Finally in § 4, we apply our result to study the autophoretic motion of spherical microparticles with catalytic caps. We study how translation velocity depends on the cap size, the surface interaction potential and the interaction length relative to particle size. Since our methodology works for both externally driven and self-propelling particles, we also study particle propulsion by both modes simultaneously, see § 5. These model problems demonstrate the wide applicability of the expressions derived by Brady (Reference Brady2021) and reproduced in this manuscript. Finally, in § 6 we summarize the key findings of our work and outline future ideas.
2. Derivation of the unified mobility expression
In this section, we derive the translation velocity (
$\boldsymbol {U}$) and rotational velocity (
$\boldsymbol {\varOmega }$) of an arbitrary particle with surface 
$S_p$ immersed in a fluid of volume 
$V$ due to an arbitrary osmophoretic body force 
$\boldsymbol {b}$, see figure 1(a). The particle surface is defined via the vector 
$\boldsymbol {x}_S$ relative to the centre of mass of the particle. We define 
$\boldsymbol {e}_r$ as the outward unit normal to the particle surface, 
$r$ is the distance from the centre of mass of the particle and 
$\boldsymbol {r}$ is the position vector defined from the centre of mass of the particle. The fluid velocity around the particle can be resolved through the modified Stokes equation, defined as
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma} + \boldsymbol{b} = \boldsymbol{0}, \end{equation}
where 
$\boldsymbol {\sigma }$ is the hydrodynamic stress tensor. The velocity field 
$\boldsymbol {u}$ is assumed to decay to zero in the far-field, 
$\boldsymbol{u}|_{S_\infty} \to {\boldsymbol{0}}$. At the particle surface, the fluid obeys a no-slip, rigid body boundary condition, 
$\boldsymbol {u}|_{S_p} = \boldsymbol {U} + \boldsymbol {\varOmega } \times \boldsymbol {x}_{S}$. To obtain 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ using Lorentz reciprocal theorem (Masoud & Stone Reference Masoud and Stone2019), we define an auxiliary Stokes flow (
$\hat {\boldsymbol {U}},\hat {\boldsymbol {\varOmega }}$) while preserving particle geometry with the same no-slip rigid surface, 
$\hat {\boldsymbol {u}}|_{S_p} = \hat {\boldsymbol {U}} + \hat {\boldsymbol {\varOmega }} \times \boldsymbol {x}_{S}$, and far-field decay, 
$\hat{\boldsymbol{u}}|_{S_\infty} \to {\boldsymbol{0}}$, boundary conditions.
Two approaches to finding the velocity of a particle of characteristic length scale a by resolving the fluid velocity at the particle surface. The particle surface is defined by a vector 
$\boldsymbol{x}_s$ relative to the centre of mass (COM) of the particle. (a) Obtain the fluid velocity near the particle surface by resolving the modified Stokes equation with an arbitrary body force, 
$\boldsymbol {b}$. The body force, 
$\boldsymbol {b}$, depends on charge, 
$\rho$, salt, 
$s$, and interaction potential, 
$\phi$. (b) When the interaction length is small 
$\lambda /a \ll 1$, the velocity near the fluid surface, at the outer edge of the interaction layer, 
$\boldsymbol {u}_{s}$, is taken to be the velocity at the particle surface. The slip velocity, 
$\boldsymbol {u}_{s}$, depends on the lumped mobility, 
$\mathcal {M}$, which depends on the interaction between the surface and solute, and the solute concentration at the vicinity of the surface, 
$c$.
Using the Lorentz reciprocal theorem, we can relate the phoretic problem (
$\boldsymbol {U}$,
$\boldsymbol {\varOmega }$,
$\boldsymbol {b}$) and the auxiliary problem (
$\hat {\boldsymbol {U}}$,
$\hat {\boldsymbol {\varOmega }}$,
$\hat {\boldsymbol {b}}$) to be
\begin{equation} \int_{S_p} \boldsymbol{e}_r\boldsymbol{\cdot}\boldsymbol{\sigma}\boldsymbol{\cdot}\hat{\boldsymbol{u}}\,{\rm d}S - \int_{S_p} \boldsymbol{e}_r\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}\boldsymbol{\cdot}\boldsymbol{u}\,{\rm d}S = \int_V \hat{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{b}\,{\rm d}V - \int_V \boldsymbol{u}\boldsymbol{\cdot}\hat{\boldsymbol{b}}\,{\rm d}V. \end{equation}
Substituting in the expressions of the fluid velocities, 
$\boldsymbol {u}|_{S_p}$ and 
$\hat {\boldsymbol {u}}|_{S_p}$, at the particle surface we can simplify (2.2). Additionally, we assume that there is no body force in the auxiliary problem, 
${\boldsymbol {\hat b}} = \boldsymbol {0}$. Thus, we can rewrite (2.2) to be
\begin{equation} \int_{S_p} \boldsymbol{e}_r\boldsymbol{\cdot}\boldsymbol{\sigma}\boldsymbol{\cdot} \left(\hat{\boldsymbol{U}}+\hat{\boldsymbol{\varOmega}}\times \boldsymbol{x}_S\right) {\rm d} S - \int_{S_p} \boldsymbol{e}_r\boldsymbol{\cdot}\hat{\boldsymbol{\sigma}}\boldsymbol{\cdot} \left(\boldsymbol{U}+\boldsymbol{\varOmega}\times\boldsymbol{x}_S\right) {\rm d}S = \int_V \hat{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{b}\,{\rm d}V. \end{equation}Since the inertia of the particle is negligible, for both the phoresis and the auxiliary problem, the particle is force-free and torque-free. For the phoretic propulsion, the hydrodynamic force and torque is balanced by the osmophoretic force and torque, or
$$\begin{gather} \underbrace{\displaystyle \int_{S_p} \boldsymbol{e}_r\boldsymbol{\cdot}\boldsymbol{\sigma}\, {\rm d}S}_\textit{hydrodynamic} \ \underbrace{- \displaystyle \int_V \boldsymbol{b}\, {\rm d}V}_\textit{osmophoretic} = \boldsymbol{0}, \end{gather}$$
$$\begin{gather}\underbrace{\displaystyle \int_{S_p} \boldsymbol{x}_s \times \boldsymbol{e}_r \boldsymbol{\cdot} \boldsymbol{\sigma} \,{\rm d}S}_\textit{hydrodynamic} \ \underbrace{- \displaystyle \int_V \boldsymbol{r}\times\boldsymbol{b}\, {\rm d}V}_\textit{osmophoretic} = \boldsymbol{0}, \end{gather}$$where the negative sign in front of the osmophoretic term appears because the osmophoretic force on the particle is equal in magnitude to the osmophoretic force on the fluid but opposite in sign (Brady Reference Brady2011). For the auxiliary system, we balance the hydrodynamic and external forces and torques, or
$$\begin{gather} \underbrace{\displaystyle \int_{S_p} \hat{\boldsymbol{\sigma}}\boldsymbol{\cdot}\boldsymbol{e}_r\,{\rm d}S }_\textit{hydrodynamic} \ + \underbrace{\hat{\boldsymbol{F}}_{ext}}_\textit{external}= \boldsymbol{0}, \end{gather}$$
$$\begin{gather}\underbrace{\displaystyle \int_{S_p} \boldsymbol{x}_s \times \hat{\boldsymbol{\sigma}}\boldsymbol{\cdot}\boldsymbol{e}_r\,{\rm d}S }_\textit{hydrodynamic} \ + \underbrace{\hat{\boldsymbol{L}}_{ext}}_\textit{external}= \boldsymbol{0}, \end{gather}$$
where 
$\hat {\boldsymbol {F}}_{ext}$ and 
$\hat {\boldsymbol {L}}_{ext}$ are the external force and torque required to move the particle in the auxiliary problem, which is to be determined. We note that (2.4)–(2.5) are different from (2.6)–(2.7) since the phoretic propulsion is induced by 
$\boldsymbol {b}$ whereas in the auxiliary problem motion is caused by 
$\hat {\boldsymbol {F}}_{ext}$ and 
$\hat {\boldsymbol {L}}_{ext}$.
To calculate 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ for a given 
$\boldsymbol {b}$ from (2.3)–(2.7), we need to express, 
$\hat {\boldsymbol {F}}_{ext}$, 
$\hat {\boldsymbol {L}}_{ext}$ and 
$\boldsymbol {\hat {u}}$ in the auxiliary problem as functions of 
$\boldsymbol {\hat {U}}$ and 
$\boldsymbol {\hat {\varOmega }}$. To do so, we use a resistance formulation to write
\begin{equation} \begin{bmatrix} \hat{\boldsymbol{F}}_{ext}\\ \hat{\boldsymbol{L}}_{ext} \end{bmatrix} = \begin{bmatrix} \boldsymbol{\mathsf{R}}_{FU} & \boldsymbol{\mathsf{R}}_{F\varOmega} \\ \boldsymbol{\mathsf{R}}_{LU} & \boldsymbol{\mathsf{R}}_{L\varOmega} \end{bmatrix}\boldsymbol{\cdot}\begin{bmatrix} \hat{\boldsymbol{U}} \\ \hat{\boldsymbol{\varOmega}} \end{bmatrix}, \end{equation}
where the resistance matrices 
$\boldsymbol{\mathsf{R}}_{FU}$, 
$\boldsymbol{\mathsf{R}}_{F\varOmega }$, 
$\boldsymbol{\mathsf{R}}_{LU}$ and 
$\boldsymbol{\mathsf{R}}_{L\varOmega }$ relate the driving force (
$\hat {\boldsymbol {F}}_{ext}$) and torque (
$\hat {\boldsymbol {L}}_{ext}$) to the translational (
$\boldsymbol {\hat {U}}$) and rotational velocity (
$\boldsymbol {\hat {\varOmega }}$). Further, we describe 
$\hat {\boldsymbol {u}}$ as
\begin{equation} \hat{\boldsymbol{u}}=\boldsymbol{\mathsf{D}}_T\boldsymbol{\cdot}\hat{\boldsymbol{U}}+ \boldsymbol{\mathsf{D}}_R\boldsymbol{\cdot}\hat{\boldsymbol{\varOmega}}\times\boldsymbol{r}, \end{equation}
where 
$\boldsymbol{\mathsf{D}}_T$ is the translation disturbance tensor and 
$\boldsymbol{\mathsf{D}}_R$ is the rotation disturbance tensor. Equations (2.8)–(2.9) combined with (2.6)–(2.7) provide necessary information to simplify (2.3) as a function of 
$\boldsymbol {\hat {U}}$ and 
$\boldsymbol {\hat {\varOmega }}$.
Next, we choose convenient values of 
$\boldsymbol {\hat {U}}$ and 
$\boldsymbol {\hat {\varOmega }}$ to simplify (2.3). Specifically, we use six auxiliary flow problems, pure translation (
$\hat {\boldsymbol {\varOmega }} = \boldsymbol {0}$ and 
$\hat {\boldsymbol {U}}$ = 
$U_0 \boldsymbol {e}_1$, 
$U_0 \boldsymbol {e}_2$, 
$U_0 \boldsymbol {e}_3$) and pure rotation (
$\hat {\boldsymbol {U}}=\boldsymbol {0}$ and 
$\hat {\boldsymbol {\varOmega }}=U_0/a \boldsymbol {e}_1$, 
$U_0/a \boldsymbol {e}_2$, 
$U_0/a \boldsymbol {e}_3$) with 
$U_0$ being the characteristic velocity scale and 
$a$ being the characteristic particle length, to obtain
\begin{equation} \begin{bmatrix}
\boldsymbol{\mathsf{R}}_{FU} & \boldsymbol{\mathsf{R}}_{F\varOmega} \\
\boldsymbol{\mathsf{R}}_{LU} & \boldsymbol{\mathsf{R}}_{L\varOmega}
\end{bmatrix}\boldsymbol{\cdot}\begin{bmatrix}
{\boldsymbol{U}} \\ {\boldsymbol{\varOmega}} \end{bmatrix}
= \begin{bmatrix} \displaystyle \int_V
\left(\boldsymbol{\mathsf{D}}_T -
\boldsymbol{I}\right)\boldsymbol{\cdot}\boldsymbol{b}\,{\rm
d}{V} \\ \displaystyle \int_V
\left(\boldsymbol{\mathsf{D}}_R -
\boldsymbol{I}\right)\boldsymbol{\cdot}\boldsymbol{r}\times\boldsymbol{b}\,{\rm
d}{V} \end{bmatrix} .
\end{equation} By inverting the resistance tensor, we obtain a mobility formulation that resolves 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ in terms of volume integrals of 
$\boldsymbol {b}$,
\begin{equation} \begin{bmatrix}
\boldsymbol{U} \\ \boldsymbol{\varOmega} \end{bmatrix} =
\begin{bmatrix} \boldsymbol{\mathsf{M}}_{UF} & \boldsymbol{\mathsf{M}}_{UL}
\\ \boldsymbol{\mathsf{M}}_{\varOmega F} & \boldsymbol{\mathsf{M}}_{\varOmega
L} \end{bmatrix}\boldsymbol{\cdot}\begin{bmatrix}
\displaystyle \int_V \left(\boldsymbol{\mathsf{D}}_T -
\boldsymbol{I}\right)\boldsymbol{\cdot}\boldsymbol{b}\,{\rm
d}{V} \\ \displaystyle \int_V
\left(\boldsymbol{\mathsf{D}}_R -
\boldsymbol{I}\right)\boldsymbol{\cdot}
\boldsymbol{r}\times\boldsymbol{b}\,{\rm d}{V}
\end{bmatrix}
\end{equation}
where the matrices 
$\boldsymbol{\mathsf{M}}_{UF}$, 
$\boldsymbol{\mathsf{M}}_{\varOmega F}$, 
$\boldsymbol{\mathsf{M}}_{UL}$ and 
$\boldsymbol{\mathsf{M}}_{\varOmega L }$ are the corresponding mobility tensors. For an in-depth mathematical analysis and mechanistic discussion regarding the various forms of 
$\boldsymbol {b}$, we redirect the reader to Brady (Reference Brady2021).
Physically, (2.11) is insightful as it helps parse apart the difference between the phoretic problem and the auxiliary problem. The rightmost term is the effective force and torque on the particle due to phoretic interactions and has two contributions: (i) the term associated with the identity tensor (
$\boldsymbol{\mathsf{I}}$) is the osmophoretic force and torque acting on the particle, and (ii) the term associated with the disturbance tensors (
$\boldsymbol{\mathsf{D}}_T$, 
$\boldsymbol{\mathsf{D}}_R$) is the hydrodynamic correction to the distribution of body forces around the particle. This correction arises because the phoretic interactions near the particle surface lead to an additional compensating fluid motion (Brady Reference Brady2011) causing a long-range hydrodynamic disturbance. This effect is not captured in the definition of the hydrodynamic mobility tensor and thus manifests separately. If the terms associated with disturbance tensors were not present, (2.11) is essentially identical to (2.8) with osmophoretic force on the particle as the external force.
We note that (2.11) takes an explicit form only when 
$\boldsymbol {b}$ is independent of 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ and is thus most convenient for systems with small Péclet number (
${Pe} \ll 1$), where 
${Pe} = {U_0 a}/{D}$, and 
$D$ is the diffusivity of the solute. The distinction from prior work, such as Stone & Samuel (Reference Stone and Samuel1996), Michelin & Lauga (Reference Michelin and Lauga2014), Lisicki et al. (Reference Lisicki, Reigh and Lauga2018), Poehnl et al. (Reference Poehnl, Popescu and Uspal2020), Poehnl & Uspal (Reference Poehnl and Uspal2021) and Ganguly & Gupta (Reference Ganguly and Gupta2023), that utilize Lorentz reciprocal theorem is that they invoke the thin interaction length limit and apply the analysis in the outer region where 
$\boldsymbol {b}=\boldsymbol {0}$, see figure 1(b). Consequently, they do not arrive at (2.11) but rather represent 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ in terms of a slip velocity at the particle surface 
$\boldsymbol {u}_s$.
We acknowledge that similar results have been presented in Khair (Reference Khair2018) and Brady (Reference Brady2021). However, in Khair (Reference Khair2018), 
$\boldsymbol {b}$ only focused on the electrophoretic contributions. In contrast, Brady (Reference Brady2021) argued that 
$\boldsymbol {b}$ should include both osmotic and phoretic contributions, and we thus refer to 
$\boldsymbol {b}$ as an osmophoretic body force. Care should be taken that the osmotic contribution only includes an excess osmotic effect since a particle cannot move without a phoretic interaction; interested readers are referred to Brady (Reference Brady2021). The phoretic contribution arises from the interaction of the particle with a macroscopically established potential field. The nature of this field depends on the specific model problem under consideration. We refer the readers to Brady (Reference Brady2021) for an in-depth mathematical analysis and a general discussion on the mechanistic origin of 
$\boldsymbol {b}$. Building on the work by Brady (Reference Brady2021), we systematically illustrate how both phoretic and osmotic contributions to the body force term are required to reconcile a broad range of results in the literature and arrive at universal mobility relationships. Additionally, through this framework, we quantify the impact of interaction length on microswimmer motion in electrolytic solutions, elaborating on the suggestion established in Brady (Reference Brady2021).
For a spherical particle (see Duprat & Stone (Reference Duprat and Stone2016) for derivation), the relevant hydrodynamic parameters are 
$\boldsymbol{\mathsf{D}}_T = {3a}/{4r}(\boldsymbol{\mathsf{I}}+\boldsymbol {e}_r\boldsymbol {e}_r)+{a^3}/{4r^3}(\boldsymbol{\mathsf{I}}-3\boldsymbol {e}_r \boldsymbol {e}_r)$, 
$\boldsymbol{\mathsf{D}}_R = ({a^3}/{r^3})\boldsymbol{\mathsf{I}}$, 
$\boldsymbol{\mathsf{M}}_{UF} = ({1}/{6{\rm \pi} \mu a})\boldsymbol{\mathsf{I}}$, 
$\boldsymbol{\mathsf{M}}_{UL}=0$, 
$\boldsymbol{\mathsf{M}}_{\varOmega F}=0$ and 
$\boldsymbol{\mathsf{M}}_{\varOmega L} = ({1}/{8{\rm \pi} \mu a^3}) \boldsymbol{\mathsf{I}}$; where 
$a$ is the radius of the sphere, 
$\mu$ is the fluid viscosity, 
$r$ is the radial distance from the centre of the sphere and 
$\boldsymbol {e}_r$ is the radial vector pointing away from the centre. Substituting, these definitions of the hydrodynamic disturbance and mobility in (2.11) we obtain
$$\begin{gather} \boldsymbol{U} = \frac{1}{6{\rm \pi}\mu a} \displaystyle \int_V \left[\left(\frac{3a}{2r}-\frac{a^3}{2r^3}-1\right)\boldsymbol{b}_\perp{+} \left(\frac{3a}{4r}+\frac{a^3}{4r^3}-1\right)\boldsymbol{b}_\parallel\right]{\rm d}{V}, \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega} = \frac{1}{8{\rm \pi}\mu a^3} \displaystyle \int_V r\left(\frac{a^3}{r^3}-1\right)\boldsymbol{e}_r \times \boldsymbol{b}_\parallel\,{\rm d}{V}, \end{gather}$$
where the body force is decomposed into 
$\boldsymbol {b}=\boldsymbol {b}_\perp +\boldsymbol {b}_\parallel$. The perpendicular subscript denotes the component normal to the sphere and the parallel subscript denotes the component parallel to the surface. Equations (2.12)–(2.13) were also reported in the prior literature for phoretic systems (Brady Reference Brady2021) as well as for different physical systems (Brenner & Cox Reference Brenner and Cox1963; Hinch Reference Hinch1991; Leal Reference Leal2007; Datt et al. Reference Datt, Zhu, Elfring and Pak2015; Elfring & Goyal Reference Elfring and Goyal2016; Datt et al. Reference Datt, Natale, Hatzikiriakos and Elfring2017). We extensively validate this result in the next section and show it relaxes to the various well-known expressions present in the literature, for both microswimmers and externally driven particles. Further, we employ this expression to study a microswimmer in the arbitrary interaction layer limit and a microswimmer driven by an external gradient in addition to its self-propelling mode of swimming.
3. Validation
3.1. Simplification at the limit of the thin interaction length scale
In this subsection, we aim to simplify (2.12)–(2.13) in the limit of the thin interaction length scale for a spherical microswimmer and recover the equations discussed in Stone & Samuel (Reference Stone and Samuel1996).
We proceed to simplify (2.12)–(2.13) at the thin interaction length limit, 
${\lambda }/{a} \ll 1$, where 
$\lambda$ is the interaction length scale. To this end, we define a stretched radial coordinate 
$\rho = {(r-a)}/{\lambda }$. Next, we expand and rewrite (2.12)–(2.13) in orders of 
${\lambda }/{a}$. Subsequently, the leading-order contribution to the translation and rotation velocities are obtained to be
$$\begin{gather} \boldsymbol{U} ={-}\frac{\lambda}{4{\rm \pi}\mu a^2} \displaystyle \int_V \rho \boldsymbol{b}_\parallel\,{\rm d}{V}, \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega} ={-}\frac{3\lambda}{8{\rm \pi}\mu a^3} \displaystyle \int_V \rho \boldsymbol{e}_r \times \boldsymbol{b}_\parallel\,{\rm d}{V} . \end{gather}$$
Since the volume of interest at the thin interaction limit is a spherical shell of thickness 
$\lambda$ surrounding the particle, we can rewrite the differential volume element to be 
${\rm d}{V} = \lambda \,{\rm d}\rho \,{\rm d}{S}$ and the volume integrals as
$$\begin{gather} \boldsymbol{U} ={-}\frac{\lambda^2}{4{\rm \pi}\mu a^2} \displaystyle \int_{S_p}\left[\displaystyle \int_0^\infty \rho \boldsymbol{b}_\parallel\,{\rm d}\rho\right]{\rm d}{S}, \end{gather}$$
$$\begin{gather}\varOmega ={-}\frac{3\lambda^2}{8{\rm \pi}\mu a^3} \displaystyle \int_{S_P} \boldsymbol{e}_r \times \left[\displaystyle \int_0^\infty \rho \boldsymbol{b}_\parallel\,{\rm d}\rho\right]{\rm d}S. \end{gather}$$In the thin interaction limit, the shear force is balanced by the parallel body force, or
\begin{equation} \frac{\mu}{\lambda^2}\frac{\partial^2 \boldsymbol{u}_\parallel}{\partial \rho^2} + \boldsymbol{b}_\parallel= 0 . \end{equation}
We note that in (3.5), since 
$\boldsymbol {b}$ is the osmophoretic force, 
$\boldsymbol {b}_\parallel$ also includes the excess osmotic term. Multiplying (3.5) by 
$\rho$ and employing integration by parts, we arrive at
\begin{equation} \displaystyle \int_0^\infty \rho \boldsymbol{b}_\parallel\,{\rm d}{\rho} = \frac{\mu}{\lambda^2} \boldsymbol{u}_\parallel |_0^\infty = \frac{\mu}{\lambda^2} \boldsymbol{u}_s, \end{equation}
where 
$\boldsymbol {u}_s = \boldsymbol {u}_{\parallel,\infty }-\boldsymbol {u}_{\parallel,0}$, is the phoretic slip velocity. Substituting (3.6) into (3.3)–(3.4), we get the widely used result derived in Stone & Samuel (Reference Stone and Samuel1996),
$$\begin{gather} \boldsymbol{U} ={-}\frac{1}{4{\rm \pi} a^2} \displaystyle \int_{S_P} \boldsymbol{u}_s\,{\rm d}{S}, \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega} ={-}\frac{3}{8{\rm \pi} a^3} \displaystyle \int_{S_P} \boldsymbol{e}_r \times \boldsymbol{u}_s\,{\rm d}{S}, \end{gather}$$where the integral is over the surface of the sphere.
3.2. Electrophoretic mobility at arbitrary interaction length scales
To further validate (2.12)–(2.13) by the determination of the electrophoretic mobility of a sphere in the Debye–Hückel limit for an arbitrary Debye length (Henry Reference Henry1931; Teubner Reference Teubner1982; Kim & Karrila Reference Kim and Karrila2013), we assume a homogeneous sphere of radius, 
$a$, immersed in a binary monovalent electrolytic solution such that the electrical permittivity of the solution is denoted as 
$\varepsilon$. Our objective is to analyse the motion of the particle with a given surface zeta potential driven by an external electric field. First, we assume that the surface zeta potential, 
$\zeta$, falls in the Debye–Hückel limit, 
${e\zeta }/{k_B T} \ll 1$, where 
$e$ is the charge of an electron, 
$k_B$ is the Boltzmann constant and 
$T$ is the absolute temperature. We also assume an electric field disturbance of 
$\boldsymbol {E}_\infty$ far away from the particle such that 
$\boldsymbol {E}_\infty = \epsilon E_0 \boldsymbol {e}_z$, where 
$\epsilon$ is a small parameter and physically indicates that the length scale of the far-field potential decay is much larger than the particle size, see figure 2(a). Note that 
$\varepsilon$ is the electrical permittivity and should not be confused with 
$\epsilon$, which is a small parameter in our analysis. Finally, the total osmophoretic body force (
$\boldsymbol {b}$) driving the particle arises through a combination of the electrostatic interaction and net excess osmotic pressure in the fluid, or
\begin{equation} \boldsymbol{b} ={-}e\left(c_{+} - c_{-}\right) \boldsymbol{\nabla} \phi - k_B T \boldsymbol{\nabla}\left(c_{+} + c_{-}\right), \end{equation}
where 
$c_{+}$ and 
$c_{-}$ are the concentrations of the positive and negative electrolytic species, respectively, and 
$\phi$ is the electric potential. As a convenient choice, we can represent the solute concentrations in terms of net charge, 
$\rho = e ( c_{+} - c_{-} )$, and salt, 
$s = c_{+} + c_{-}$, and rewrite the body force to be 
$\boldsymbol {b} = - \rho \boldsymbol {\nabla } \phi - k_B T \boldsymbol {\nabla } s$. It should be noted that care should be exercised in choosing the appropriate expression for 
$\boldsymbol {b}$. Specifically, the osmotic contribution 
$-k_B T \boldsymbol {\nabla }s$ refers to the excess osmotic contribution arising out of an interaction that locally drives the solute out of equilibrium. This effectively implies that in the absence of such interactions, an external salt gradient on its own cannot induce net particle motion, as demonstrated in Brady (Reference Brady2021). The equivalence of (2.12) and the results of Brady (Reference Brady2021) can be seen by defining an additional surface stress contribution, 
$\boldsymbol {\sigma }_p = -k_B T s \boldsymbol{\mathsf{I}}$, as per (2.17) in Brady (Reference Brady2021), due to the excess osmotic pressure. The divergence of 
$\boldsymbol {\sigma }_p$ leads to the second term in (3.9), 
$-k_B T \boldsymbol {\nabla } s$. We refer the readers to our discussion on the mechanistic origin of 
$\boldsymbol {b}$ in § 2 and redirect them to Brady (Reference Brady2021) for more details.
Methodology to validate proposed mobility expressions for a charged particle with a zeta potential (
$\zeta$) in the Debye–Hückel limit for (a) electrophoresis with an external field 
$\boldsymbol {E}_\infty =\epsilon E_0 \boldsymbol {e}_z$ and (b) diffusiophoresis with externally imposed solute gradient 
$\boldsymbol {\nabla }s_\infty = 2\epsilon c_0 \boldsymbol {e}_z$. The expressions of dimensionless osmophoretic force 
$\tilde {\boldsymbol {b}}$ are provided. Substituting the appropriate 
$\boldsymbol {b}$ in (2.12) enables us to recover mobility relationships that otherwise require cumbersome calculations.
To appropriately derive the particle motion, we are required to obtain the solutions to 
$\rho$, 
$s$ and 
$\phi$ for a given 
$\boldsymbol {E}_\infty$ and 
$\zeta$. As mentioned earlier, we assume that the species are monovalent, 
$z_\pm =\pm 1$ and have diffusivities 
$D_\pm$, the species balance is given by the steady Nernst–Planck equations,
$$\begin{gather} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}c_+=D_+ \left[\nabla^2 c_{+} + \frac{e}{k_B T} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(c_{+}\boldsymbol{\nabla}\phi\right)\right] , \end{gather}$$
$$\begin{gather}\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}c_-= D_- \left[\nabla^2 c_{-} - \frac{e}{k_B T} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(c_{-}\boldsymbol{\nabla}\phi\right)\right]. \end{gather}$$
For 
${Pe} = aU/D \ll 1$ (
$U$ is the velocity scale for the particle, 
$D = {2 D_+ D_-}/({D_+ + D_-})$ is the ambipolar diffusivity), we ignore the convective effects. Consequently, (3.10) and (3.11) can be rewritten in terms of 
$\rho$ and 
$s$ as
$$\begin{gather} \nabla^2 s + \frac{1}{k_B T} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\rho\boldsymbol{\nabla}\phi\right) = 0, \end{gather}$$
$$\begin{gather}\nabla^2\rho + \frac{e^2}{k_B T} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(s\boldsymbol{\nabla}\phi\right) = 0. \end{gather}$$Finally, the system of equations is closed by using Poisson's equation to resolve the electric potential,
\begin{equation} - \varepsilon \nabla^2 \phi = \rho. \end{equation}
In the far-field, at 
$r\to \infty$, the potential gradient is the externally imposed electric field and the fluid is electroneutral, or
$$\begin{gather} \left. -\boldsymbol{\nabla}\phi \right| _{r \rightarrow \infty}= \epsilon E_0 \boldsymbol{e}_z. \end{gather}$$
$$\begin{gather}\left. \rho \right|_{r \rightarrow \infty} = 0. \end{gather}$$
Moreover, in the far-field 
$c_{+} = c_{-} = c_0$, where 
$c_0$ is a characteristic solute concentration, we can write 
$s$ to follow
\begin{equation} s |_{r \rightarrow \infty} = 2c_0. \end{equation}
At the particle surface, at 
$r=a$, the electrostatic potential is equal to the zeta potential at the surface, or
\begin{equation} \phi |_{r = a} = \zeta. \end{equation}Additionally, there is no salt or charge flux normal to the particle surface,
$$\begin{gather} \boldsymbol{e}_r\boldsymbol{\cdot}\left[\boldsymbol{\nabla}s + \frac{\rho}{k_B T} \boldsymbol{\nabla}\phi\right]_{r = a} = 0, \end{gather}$$
$$\begin{gather}\boldsymbol{e}_r\boldsymbol{\cdot}\left[\boldsymbol{\nabla}\rho + \frac{s e^2}{k_B T} \boldsymbol{\nabla}\phi\right]_{r = a} = 0. \end{gather}$$Equations (3.12)–(3.20) are non-dimensionalized using the following appropriate scales:
\begin{equation} \tilde{\boldsymbol{\nabla}} = a \boldsymbol{\nabla},\quad \tilde{\nabla}^2 = a^2 \nabla^2,\quad \tilde{\phi} = \frac{e\phi}{k_B T},\quad\tilde{\rho} = \frac{\rho}{e c_0},\quad \tilde{s} = \frac{s}{c_0}, \quad\tilde{r} = \frac{r}{a}. \end{equation}Thus, the non-dimensional Poisson–Nernst–Planck equations are given as
$$\begin{gather} \tilde{\nabla}^2\tilde{s} + {\tilde{\boldsymbol{\nabla}}}\boldsymbol{\cdot} \left(\tilde{\rho}\tilde{\boldsymbol{\nabla}}\tilde{\phi}\right) = 0, \end{gather}$$
$$\begin{gather}\tilde{\nabla}^2\tilde{\rho} + \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot} \left(\tilde{s}\tilde{\boldsymbol{\nabla}}\tilde{\phi}\right) = 0, \end{gather}$$
$$\begin{gather}\tilde{\nabla}^2 \tilde{\phi} ={-}\frac{\kappa^2}{2} \tilde{\rho}, \end{gather}$$
where 
$\kappa = ({2 a^2 e^2 c_0 }/{\varepsilon k_B T})^{1/2}$ is the dimensionless inverse Debye length. In the far-field, thus
$$\begin{gather} -\left. \tilde{\boldsymbol{\nabla}}\tilde{\phi} \right|_{\tilde{r} \rightarrow \infty} = \epsilon \tilde{E}_0 \boldsymbol{e}_z, \end{gather}$$
$$\begin{gather}\left. \tilde{\rho} \right|_{\tilde{r} \rightarrow \infty} = 0 , \end{gather}$$
$$\begin{gather}\left. \tilde{s}\right|_{\tilde{r} \rightarrow \infty} = 2 , \end{gather}$$
where 
$\tilde{E}_0 = a e E_0/(k_B T )$ is the non-dimensional electric-field. Similarly, at the particle surface, the non-dimensional boundary conditions read
$$\begin{gather} \boldsymbol{e}_r \boldsymbol{\cdot} \left[\tilde{\boldsymbol{\nabla}}\tilde{s} + \tilde{\rho} \tilde{\boldsymbol{\nabla}}\tilde{\phi}\right]_{\tilde{r} = 1} = 0, \end{gather}$$
$$\begin{gather}\boldsymbol{e}_r \boldsymbol{\cdot} \left[\tilde{\boldsymbol{\nabla}}\tilde{\rho} + \tilde{s} \tilde{\boldsymbol{\nabla}} \tilde{\phi}\right]_{\tilde{r} = 1} = 0, \end{gather}$$
$$\begin{gather}\left. \tilde{\phi} \right|_{\tilde{r} = 1} = \frac{e \zeta}{k_B T} = \tilde{\zeta}. \end{gather}$$As we discuss later, it is more appropriate to write (3.30) as a constant charge boundary condition, which renders the gradient of the potential to be constant instead. However, for the weak-field, these surface boundary conditions are equivalent and hence we retain the constant potential boundary condition for simplicity.
For the remainder of the calculation until (3.67), we will drop the tilde superscript in (3.22)–(3.30) for convenience and restore the dimensions once the non-dimensional calculations are complete. We expand 
$\phi$, 
$\rho$ and 
$s$ in the small parameters 
$\zeta$ and 
$\epsilon$ as
$$\begin{gather} \phi = \phi_{00} + \zeta \phi_{01} + \epsilon \left(\phi_{10} + \zeta\phi_{11}\right), \end{gather}$$
$$\begin{gather}\rho = \rho_{00} + \zeta \rho_{01} + \epsilon\left(\rho_{10} + \zeta \rho_{11}\right), \end{gather}$$
$$\begin{gather}s = s_{00} + \zeta s_{01} + \epsilon \left(s_{10} + \zeta s_{11}\right). \end{gather}$$The asymptotic expansions in (3.31)–(3.33) are substituted into (3.22)–(3.30), and the corresponding equations are solved at each asymptotic order.
Order 
$O(1)$: an uncharged particle without any electric field. The governing equations and boundary conditions are obtained to be
$$\begin{gather} \nabla^2 s_{00} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\rho_{00} \boldsymbol{\nabla}\phi_{00}\right) = 0, \end{gather}$$
$$\begin{gather}\nabla^2 \rho_{00} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(s_{00}\boldsymbol{\nabla}\phi_{00}\right) = 0, \end{gather}$$
$$\begin{gather}\nabla^2 \phi_{00} ={-}\frac{\kappa^2}{2} \rho_{00}, \end{gather}$$
$$\begin{gather}\boldsymbol{e}_r\boldsymbol{\cdot}\left[\boldsymbol{\nabla}s_{00}+\rho_{00}\boldsymbol{\nabla}\phi_{00}\right]=0,\quad\mathrm{at}\ r=1, \end{gather}$$
$$\begin{gather}\boldsymbol{e}_r \boldsymbol{\cdot} \left[\boldsymbol{\nabla} \rho_{00} + s_{00} \boldsymbol{\nabla}\phi_{00}\right]=0,\quad \mathrm{at}\ r=1, \end{gather}$$
$$\begin{gather}\phi_{00} = 0,\quad \mathrm{at}\ r=1, \end{gather}$$
$$\begin{gather}s_{00}=2,\quad \mathrm{at}\ r\to\infty, \end{gather}$$
$$\begin{gather}\rho_{00}=0,\quad \mathrm{at}\ r\to\infty, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\phi_{00}=0,\quad\mathrm{at}\ r\to\infty. \end{gather}$$ The system of (3.34)–(3.40) have a trivial solution, i.e. 
$\phi _{00}=0$, 
$\rho _{00}=0$, 
$s_{00}=2$. Physically, the solution simply implies that the ion concentration is uniform because the particle is uncharged and there is no electric field.
Order 
$O(\epsilon )$: perturbation due to the electric field for an uncharged particle. The governing equations for the salt and charge dynamics, and electrostatic potential after substituting the expressions of 
$\rho _{00}$, 
$s_{00}$, and 
$\phi _{00}$ are given as
$$\begin{gather} \nabla^2 s_{10}=0, \end{gather}$$
$$\begin{gather}\nabla^2 \rho_{10} + 2 \nabla^2 \phi_{10}=0, \end{gather}$$
$$\begin{gather}\nabla^2 \phi_{10} ={-}\frac{\kappa^2}{2} \rho_{10}. \end{gather}$$
The corresponding reduced boundary conditions when 
$r \to \infty$ are
$$\begin{gather} s_{10} = 0, \end{gather}$$
$$\begin{gather}\rho_{10}=0, \end{gather}$$
$$\begin{gather}-\boldsymbol{\nabla}\phi_{10} = E_0 \boldsymbol{e}_z, \end{gather}$$
and when 
$r=1$, they read
$$\begin{gather} \boldsymbol{e}_r\boldsymbol{\cdot}\boldsymbol{\nabla}s_{10}=0, \end{gather}$$
$$\begin{gather}\boldsymbol{e}_r\boldsymbol{\cdot}\left[\boldsymbol{\nabla}\rho_{10}+2\boldsymbol{\nabla}\phi_{10}\right]=0. \end{gather}$$
$$\begin{gather}\boldsymbol{e}_r \boldsymbol{\cdot} \boldsymbol{\nabla} \phi_{10}=0, \end{gather}$$
where the derivative of the potential is set to be zero to ensure that there is no excess charge on the surface, see discussion below (3.30). Mathematically, at 
$O(\epsilon )$, the charge on the particle surface is zero. Through Gauss's law, the no surface charge boundary condition necessitates that 
$\boldsymbol {e}_r \boldsymbol{\cdot}\boldsymbol {\nabla }\phi _{10} = 0$. Hence, (3.50) implies that 
$\boldsymbol {e}_r\boldsymbol{\cdot}\boldsymbol {\nabla }\rho _{10}=0$ at the particle surface. Along with (3.44) and (3.47), we obtain 
$\rho _{10}=0$. Similarly, (3.43), (3.46) and (3.49) reveal 
$s_{10}=0$.
Since 
$\rho _{10}=0$, 
$\phi _{10}$ is governed by the Laplace equation, and the solution with appropriate boundary conditions reads (Griffiths Reference Griffiths2005)
\begin{equation} \phi_{10}\left(r,\theta\right) ={-}E_0 r \left(1+\frac{1}{2r^3}\right) \cos \theta. \end{equation}Physically, this order implies that the perturbed potential and corresponding electric field lines get modified due to the geometry of the particle but there is no charge and salt accumulation.
Order 
$O(\zeta )$: perturbation of a charged particle without an external electric field. The equations governing the charge, salt and potential at 
$O(\zeta )$ are analogous to (3.43)–(3.45) at 
$O(\epsilon )$ but have different boundary conditions. The governing equations are
$$\begin{gather} \nabla^2 s_{01} = 0, \end{gather}$$
$$\begin{gather}\nabla^2 \rho_{01} + 2 \nabla^2 \phi_{01} =0, \end{gather}$$
$$\begin{gather}\nabla^2 \phi_{01} ={-}\frac{\kappa^2}{2}\rho_{01}. \end{gather}$$
As 
$r \to \infty$ the net charge, salt and electric potential gradient all decay to zero, or
$$\begin{gather} s_{01} = 0, \end{gather}$$
$$\begin{gather}\rho_{01} = 0, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \phi_{01} = 0. \end{gather}$$
At the particle surface, 
$r=1$, we obtain
$$\begin{gather} \boldsymbol{e}_r\boldsymbol{\cdot}\boldsymbol{\nabla}s_{01} = 0, \end{gather}$$
$$\begin{gather}\boldsymbol{e}_r\boldsymbol{\cdot}\left[\boldsymbol{\nabla}\rho_{01} + 2\boldsymbol{\nabla}\phi_{01}\right]=0, \end{gather}$$
$$\begin{gather}\phi_{01}=1. \end{gather}$$
Equations (3.53), (3.56) and (3.59) yield 
$s_{01} = 0$. The remainder of the equations reveal
$$\begin{gather} \phi_{01}(r) = \frac{1}{r} \exp({-\kappa(r-1)}) , \end{gather}$$
$$\begin{gather}\rho_{01}(r) ={-}\left(\frac{2}{r}\right)\exp({-\kappa(r-1)}) . \end{gather}$$Physically, the results indicate the distribution of potential and charge arising due to a charged particle.
Order 
$O(\epsilon \zeta )$: perturbation of both imposed electric field and surface zeta potential. To fully resolve the body force to the order of 
$O(\epsilon \zeta )$, we have to obtain 
$s_{11}$, this can be observed by the expansion of (3.9) and collecting the respective orders. The governing equations and boundary conditions for salt, 
$s_{11}$, are given by
$$\begin{gather} \nabla^2 s_{11} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\rho_{01}\boldsymbol{\nabla} \phi_{10}\right) = 0, \end{gather}$$
$$\begin{gather}s_{11} \to 0,\quad \mathrm{at}\ r \to \infty, \end{gather}$$
$$\begin{gather}\boldsymbol{e}_r \boldsymbol{\cdot} \left[\boldsymbol{\nabla}s_{11} + \rho_{01}\boldsymbol{\nabla}\phi_{10}\right]=0\quad \mathrm{at}\ r=1. \end{gather}$$
Substituting in the expression of 
$\rho _{01}$ from (3.63) and 
$\phi _{10}$ from (3.52) we can solve for 
$s_{11}(r,\theta )$. The salt dynamics are in the form of 
$s_{11}=E_0\,f(r)\cos \theta$ where 
$f(r)$ is
\begin{align} f(r) &=
\frac{1}{r^2}\left(\frac{\kappa}{6}-\frac{5}{3}-\frac{10}{\kappa}-\frac{2}{\kappa^2}\right)+
\frac{\kappa^2}{3}\exp\left({\kappa}\right)\,{Ei} \left(-\kappa
r\right)+\frac{2}{\kappa^2 r^2}
\exp({-\kappa(r-1)})\nonumber\\
&\quad-\frac{4}{3\kappa r^2}\exp({-\kappa(r-1)})+\frac{ \kappa}{3}
\exp({-\kappa(r-1)})\nonumber\\
&\quad -\frac{1}{3r}\exp({-\kappa(r-1)})+\frac{2}{\kappa
r} \exp({-\kappa(r-1)}),
\end{align}
where 
${Ei}()$ is the elliptic integral.
Reintroducing dimensions: at this stage, we restore the dimensions and reintroduce the tilde for dimensionless variables. The total osmophoretic body force 
$\boldsymbol {b}$ is made dimensionless by writing 
$\boldsymbol {b} = ({k_B T c_0}/{a} )\tilde { \boldsymbol {b} } = ({ \varepsilon (k_B T )^2 \kappa ^2 }/{2 e^2 a^3} )\tilde {\boldsymbol {b}}$. The first relevant order of 
$\tilde {\boldsymbol {b}}$ for electrophoresis is 
$\epsilon \tilde {\zeta }$ because it is the order at which a charged particle is being driven by an electric field. To this end, we write
\begin{equation} \tilde{\boldsymbol{b}} ={-} \epsilon \tilde{\zeta} \left[ \tilde{\rho}_{00} \tilde{\boldsymbol{\nabla}} \tilde{\phi}_{11} + \tilde{\rho}_{11} \tilde{\boldsymbol{\nabla}} \tilde{\phi}_{00} + \tilde{\rho}_{01} \tilde{\boldsymbol{\nabla}} \tilde{\phi}_{10} + \tilde{\rho}_{10} \tilde{\boldsymbol{\nabla}} \tilde{\phi}_{01} + \tilde{\boldsymbol{\nabla}} \tilde{s}_{11} \right]. \end{equation}
Based on the solutions at different orders, it is straightforward to see that 
$\boldsymbol {\tilde {b}}$ reduces to
\begin{equation} \boldsymbol{\tilde{b}} ={-}\epsilon \tilde{\zeta} \left[ \tilde{\rho}_{01} \tilde{\boldsymbol{\nabla}} \tilde{\phi}_{10} + \tilde{\boldsymbol{\nabla}} \tilde{s}_{11} \right]. \end{equation}
We note that the body force term of 
$\tilde {s}_{11}$ integrates out to zero in the calculation of 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ for electrophoresis and is generally not included in prior analyses. However, we retain this term for consistency as it does become crucial for diffusiophoretic phenomena, as we detail in § 3.3.
After substituting the values of 
$\tilde {\rho }_{01}$, 
$\tilde {\phi }_{10}$ and 
$\tilde {s}_{11}$, the resultant equation in dimensional form is
\begin{align} \boldsymbol{b}&={-}\frac{\epsilon \varepsilon E_0 \zeta \kappa^2} {2 a^2} \left\{\left[\frac{2}{\tilde{r}}\exp({-\kappa(\tilde{r}-1)})\left(1-\frac{1}{\tilde{r}^3} \right)\cos\theta+\frac{{\rm d} f(\tilde{r})}{{\rm d} \tilde{r}}\cos\theta\right]\boldsymbol{e}_r \right. \nonumber\\ &\quad \left.-\left[\frac{2}{\tilde{r}}\exp({-\kappa(\tilde{r}-1)}) \left(1+\frac{1}{2 \tilde{r}^3}\right)\sin\theta+\frac{f(\tilde{r})}{\tilde{r}}\sin \theta\right]\boldsymbol{e}_\theta \right\}. \end{align}Equation (3.70) is substituted into (2.12) to obtain the translational velocity to be
\begin{equation} \boldsymbol{U}=\mathcal{M} \boldsymbol{E}_{\infty}, \end{equation}
where the mobility 
$\mathcal {M}$ after restoring dimensions is
\begin{equation} \mathcal{M}=\frac{\varepsilon\zeta}{6\mu}\left[(1+\kappa)+(12-\kappa^2)\displaystyle \int_1^\infty \frac{e^{\kappa (1-t)}}{t^5} {\rm d}t\right], \end{equation}
which is the seminal result of Henry (Reference Henry1931) for arbitrary double layer thickness, and has also been reported by Teubner (Reference Teubner1982) and Kim & Karrila (Reference Kim and Karrila2013). For a homogeneous sphere, our analysis reveals 
$\boldsymbol {\varOmega }=\boldsymbol {0}$, as expected. We emphasize that it is straightforward to extend the calculations to heterogeneous spheres (Velegol et al. Reference Velegol, Anderson and Garoff1996; Teubner Reference Teubner1982) and obtain results for electrorotation in arbitrary double-layer thicknesses, which otherwise requires considerable effort.
3.3. Electrolytic diffusiophoretic mobility at arbitrary interaction length scales
Next, we focus on the process of electrolytic diffusiophoresis in the Debye–Hückel limit and for arbitrary double-layer thickness. We assume that the external concentration gradient of a binary monovalent electrolyte is given as 
$\boldsymbol {\nabla } s_{\infty } = 2 \epsilon \boldsymbol {\nabla } c_0$, where 
$\epsilon$ is a small parameter, much like in § 3.2, see figure 2(b). Here, 
$\boldsymbol {b}$ is required to be expanded to an additional higher order of 
$\epsilon \tilde {\zeta }^2$. The term of order 
$O ( \epsilon \tilde {\zeta } )$ is identical to electrophoresis and represents the electrophoretic component of the diffusiophoretic mobility. The second term of order 
$O ( \epsilon \tilde {\zeta }^2 )$ denotes the chemiphoretic component. We employ the expression of 
$\boldsymbol {b}$, derived in this subsection, to (2.12)–(2.13). This allows us to retrieve the expression of the translation velocity of a charged spherical particle in an unbounded solution of a symmetrically charged electrolyte for an arbitrary double-layer thickness, which otherwise requires considerable efforts, see Keh & Wei (Reference Keh and Wei2000).
We acknowledge the electrokinetic equations used to describe such diffusiophoretic systems are analogous to our treatment of the motion of electrophoretically propelled particles in § 3.2. However, the key mechanistic difference is the presence of an external gradient of solute 
$\boldsymbol {\nabla } s_{\infty }$ instead of an imposed electric field 
$\boldsymbol {E}_\infty$; this results in a change of boundary conditions and subsequently the solutions at different asymptotic orders. To preserve the pedagogical nature of our manuscript, we will rederive the electrophoretic contribution and subsequently solve for the chemiphoretic contribution to the osmophoretic body force term and attempt to emphasize key physical and mathematical differences between the derivations laid out in §§ 3.2 and 3.3.
Consider a colloidal particle with a surface zeta potential, 
$\zeta$, in an external solute gradient of a symmetric binary electrolyte. We assume that the electrolytes are monovalent such that 
$z_\pm = \pm 1$. The ions are assumed to have different diffusivities 
$D_+ \neq D_-$. The governing equations of the concentration of the ionic species, 
$c_\pm$, are identical to (3.10)–(3.11) and the interaction potential is governed by Poisson's equation, as given in (3.14). However, the far-field boundary conditions are different. Specifically, as 
$r \to \infty$, we assume that the concentration of the ionic species is linear with position 
$z$, or
\begin{equation} c_\pm= c_0\left(1+ \frac{\epsilon z}{a}\right). \end{equation}Further, it is assumed that the electric current in the far-field is zero, which yields (Prieve et al. Reference Prieve, Anderson, Ebel and Lowell1984; Velegol et al. Reference Velegol, Garg, Guha, Kar and Kumar2016; Gupta et al. Reference Gupta, Rallabandi and Stone2019)
\begin{equation} - {\boldsymbol{\nabla}} \phi |_{r\to\infty} = \epsilon \beta \frac{k_B T}{a e} = \epsilon E_0 \boldsymbol{e}_z, \end{equation}
where 
$\beta = ({D_+ - D_-})/({D_+ + D_-})$ and 
$E_0 = \beta ({k_B T}/{a e})$. The electric field is thus induced due to unequal diffusivities and a non-zero salt gradient.
Similar to electrophoresis, modification of the governing equations in terms of charge, 
$\rho = e(c_+ - c_-)$, salt, 
$s=c_+ + c_-$, and potential, 
$\phi$, result in (3.12)–(3.14). In the far-field, the boundary conditions for 
$\phi$ and 
$\rho$ are identical to (3.15) and (3.16) with the aforementioned definition of 
$E_0$. However, the boundary condition of 
$s$ is modified to
\begin{equation} s|_{r \to\infty } = 2c_0 \left( 1+\frac{\epsilon z}{a} \right). \end{equation}
Note that 
$\epsilon /a = \boldsymbol {\nabla } \log s_\infty$. The boundary conditions are identical at the particle surface, i.e. (3.18)–(3.20). The objective is to solve 
$\rho$, 
$s$ and 
$\phi$ with the modified boundary conditions above and subsequently evaluate the total osmophoretic body force, following the procedure used to obtain (3.9). We non-dimensionalize the equations using the same scales as (3.21a–f) and also define 
$\boldsymbol {b} = ({k_B T c_0}/{a} )\tilde {\boldsymbol {b}} = ({\varepsilon (k_B T)^2 \kappa ^2}/{2 e^2 a^3}) \tilde {\boldsymbol {b}}$, where the definition of 
$\kappa$ is also identical.
For simplicity, we drop the tilde from our analysis until (3.92) and reintroduce them afterward. Thus, the non-dimensional osmophoretic body force is given as 
$\boldsymbol {b} = -\rho \boldsymbol {\nabla }\phi -\boldsymbol {\nabla }s$. We expand 
$\rho$, 
$s$ and 
$\phi$ until 
$O(\epsilon \zeta ^2)$, and solve the equations at each order.
Order 
$O(1)$: an uncharged particle without any external salt gradient. The results at this order are identical to electrophoresis and thus yield 
$s_{00}=2$, 
$\rho _{00}=0$, and 
$\phi _{00}=0$, indicating a uniform concentration of ion with no charge and potential.
Order 
$O(\epsilon )$: perturbation of the external salt concentration to an uncharged particle. This order is distinct compared with electrophoresis since the far-field boundary condition for salt is different, while the remainder of the equations and boundary conditions are identical. We note that the boundary condition for the electric field is similar to electrophoresis since we have defined 
$E_0$. The solution simply reduces to zero and uniform charge density 
$\rho _{10}=0$, while both salt and potential follow the Laplace equation. The results read
$$\begin{gather} s_{10}\left(r,\theta\right)=2r\left(1+\frac{1}{2r^3}\right)\cos\theta, \end{gather}$$
$$\begin{gather}\phi_{10} ={-}E_0 r \left(1+\frac{1}{2r^3}\right)\cos\theta. \end{gather}$$
Physically, at this order, a gradient in the salt concentration far away perturbs the salt field and induces a potential field if diffusivity asymmetry is present (
$E_0 \neq 0$ only when 
$\beta \neq 0$). However, since the surface is uncharged, 
$\rho _{10}=0$.
Order 
$O(\zeta )$: perturbation in the surface charge of the particle without an external field. Since there is no external field at this order, the solution is identical to electrophoresis with 
$s_{01}=0$ and
$$\begin{gather} \phi_{01}(r) = \frac{1}{r}\exp({-\kappa(r-1)}), \end{gather}$$
$$\begin{gather}\rho_{01}(r) ={-}\frac{2}{r} \exp({-\kappa(r-1)}). \end{gather}$$This order represents the potential and charge profiles due to the surface charge of the particle. However, there is no salt accumulation at this order since the reduction in the coion concentration is balanced by the increase in the counter-ion concentration.
Order 
$O(\epsilon \zeta )$: perturbation in both the imposed salt concentration and surface charge. The governing equations for charge (
$\rho _{11}$), salt (
$s_{11}$), and potential (
$\phi _{11}$) are
$$\begin{gather} \nabla^2 s_{11} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left[\rho_{01}\boldsymbol{\nabla}\phi_{10} \right]=0, \end{gather}$$
$$\begin{gather}\nabla^2 \rho_{11}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left[s_{10}\boldsymbol{\nabla}\phi_{01}+ 2 \boldsymbol{\nabla}\phi_{11}\right]=0, \end{gather}$$
$$\begin{gather}\nabla^2 \phi_{11} ={-}\frac{\kappa^2}{2}\rho_{11}. \end{gather}$$
The boundary conditions at the particle surface, 
$r=1$, are
$$\begin{gather} \boldsymbol{e}_r\boldsymbol{\cdot}\left[\boldsymbol{\nabla}s_{11}+\rho_{01}\boldsymbol{\nabla}\phi_{10}\right]=0, \end{gather}$$
$$\begin{gather}\boldsymbol{e}_r\boldsymbol{\cdot}\left[\boldsymbol{\nabla}\rho_{11}+s_{10}\boldsymbol{\nabla}\phi_{01}+2\boldsymbol{\nabla}\phi_{11}\right]=0, \end{gather}$$
$$\begin{gather}\boldsymbol{e}_r\boldsymbol{\cdot}\boldsymbol{\nabla} \phi_{11}=0. \end{gather}$$
Again, there is no external field as 
$r\to \infty$. Substituting in the solutions obtained in 
$O(\epsilon )$, and 
$O(\zeta )$ we begin to solve 
$s_{11}$, 
$\rho _{11}$ and 
$\phi _{11}$. Further, we can separate the 
$r$ and 
$\theta$ contributions by redefining 
$s_{11}(r, \theta ) = f_{s_{11}}(r)\cos \theta$, 
$\rho _{11}(r,\theta )=f_{\rho _{11}}(r)\cos \theta$, and 
$\phi _{11}(r, \theta )=f_{\phi _{11}}(r)\cos \theta$ and solve the equations numerically, see Appendix A. It is possible to find analytical solutions to (A1)–(A6), similar in form to (3.67). In the scope of our current paper, we choose to resolve the dynamics at 
$O(\epsilon \zeta )$ and 
$O(\epsilon \zeta ^2)$ numerically. For an analytical derivation of such higher-order effects the reader is directed to Keh & Wei (Reference Keh and Wei2000).
Order 
$O(\epsilon \zeta ^2)$: first-order perturbation in salt field and second-order perturbation in surface charge. We only seek to solve 
$s_{12}$ at this order since it is the only quantity required to resolve the body force up to 
$O(\epsilon \zeta ^2)$, see (3.91). The equation governing the dynamics of 
$s_{12}$ is
\begin{equation} \nabla^2 s_{12} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left[\rho_{01}\boldsymbol{\nabla}\phi_{11}+\rho_{11}\boldsymbol{\nabla}\phi_{01}\right]=0, \end{equation}with
\begin{equation} \boldsymbol{e}_r\boldsymbol{\cdot}\left[\boldsymbol{\nabla}s_{12}+\rho_{01}\boldsymbol{\nabla}\phi_{11}+\rho_{11} \boldsymbol{\nabla}\phi_{01}\right]_{r=1}=0, \end{equation}
and 
$s_{12}(r \rightarrow \infty, \theta ) =0$. We write 
$s_{12}$ as 
$s_{12}(r,\theta )=f_{s_{12}}(r)\cos \theta$ and solve the equations numerically, see Appendix A. As discussed previously, we solve (A7), (A8), and the far-field constraint numerically.
Restoring dimensions: we now restore dimensions and reintroduce tilde to describe dimensionless variables. Therefore, we write the body force 
$\boldsymbol {b} = ({\varepsilon (k_B T)^2 \kappa ^2}/{2 e^2 a^3} )\tilde {\boldsymbol {b}}$ such that
$$\begin{gather} \tilde{\boldsymbol{b}} = \epsilon \tilde{\zeta} \boldsymbol{\tilde{b}}_{11} + \epsilon \tilde{\zeta}^2 \boldsymbol{\tilde{b}}_{12}, \end{gather}$$
$$\begin{gather}\tilde{\boldsymbol{b}}_{11} ={-} \left( \tilde{\rho}_{01} \boldsymbol{\tilde{\nabla}}\tilde{\phi}_{10}+\boldsymbol{\tilde{\nabla}}s_{11} \right), \end{gather}$$
$$\begin{gather}\tilde{\boldsymbol{b}}_{12} ={-} \left( \tilde{\rho}_{01}\boldsymbol{\tilde{\nabla}} \tilde{\phi}_{11}+\tilde{\rho}_{11}\boldsymbol{\tilde{\nabla}}\tilde{\phi}_{01}+\boldsymbol{\tilde{\nabla}}\tilde{s}_{12} \right). \end{gather}$$ At this point, some comments are in order. We note that 
$\tilde { \boldsymbol {b}}$ could also include a term at 
$O(\epsilon )$ since 
$s_{10} \neq 0$. However, 
$\tilde {\boldsymbol {b}}$ only includes excess osmotic pressure and not the osmotic pressure itself. This is because the osmotic pressure contribution due to 
$\boldsymbol {\nabla }s_\infty$ would lead to particle motion even in the absence of phoretic interactions. Only the terms at subsequent orders are included to ignore this effect. Further, we highlight that both 
$\tilde {\phi }_{10}$ and 
$\tilde {s}_{11}$ are proportional to 
$\tilde {E}_0$. Therefore, the 
$O (\epsilon \tilde {\zeta } )$ term only depends on 
$E_0$ and is referred to as the electrophoretic contribution. Since 
$\tilde {\rho }_{01}$, 
$\tilde {\phi }_{10}$ and 
$\tilde {s}_{11}$ are identical to electrophoretic solution, the 
$O(\epsilon \tilde {\zeta })$ is equal to the one described earlier in (3.72).
However, in contrast, for the 
$O(\epsilon \tilde {\zeta }^2)$ contribution, 
$\tilde {\rho }_{01}$, 
$\tilde {\phi }_{11}$, 
$\tilde {\rho }_{11}$, 
$\tilde {\phi }_{01}$ and 
$\tilde {s}_{12}$ are independent of 
$E_0$. Furthermore, 
$\tilde {\phi }_{11}$, 
$\tilde {\rho }_{11}$ and 
$\tilde {s}_{12}$ are all proportional to 
$\tilde {s}_{10}$, which is consistent with the prior literature (Anderson Reference Anderson1989; Keh & Wei Reference Keh and Wei2000; Gupta et al. Reference Gupta, Rallabandi and Stone2019).
Since 
$\tilde {\boldsymbol {b}}_{11}$ is identical to the electrophoretic motion, we focus our attention on 
$\tilde {\boldsymbol {b}}_{12}$, which reads
\begin{equation} \tilde{\boldsymbol{b}}_{12}={-} \left(\rho_{01}\frac{{\rm d}f_{\phi_{11}}}{{\rm d}r}+f_{\rho_{11}}\frac{{\rm d}\phi_{01}}{{\rm d}r}+\frac{{\rm d}f_{s_{12}}}{{\rm d}r} \right) \cos\theta \boldsymbol{e}_r - \left(\frac{\rho_{01} f_{\phi_{11}}}{r}+\frac{f_{s_{12}}}{r}\right)\sin\theta \boldsymbol{e}_\theta . \end{equation} Upon substituting the value of 
$\boldsymbol {b}$ in (2.12), the translation velocity of the particle could be simplified to read
\begin{equation} \boldsymbol{U} = \mathcal{M} \boldsymbol{\nabla} \log s_{\infty}, \end{equation}where
\begin{equation} \mathcal{M} = \frac{\varepsilon}{\mu}\left[\frac{k_B T}{e}\beta \zeta \varTheta_1\left(\kappa\right) + \frac{\zeta^2}{8}\varTheta_2\left(\kappa\right)\right], \end{equation}
where 
$\varTheta _1$ and 
$\varTheta _2$ are evaluated numerically, see Appendix A. Figure 3 demonstrates good quantitative agreement between the values obtained in Keh & Wei (Reference Keh and Wei2000) and our results.
Comparison of proposed mobility expressions of diffusiophoretic mobility in (3.93) with the mobility reported in Keh & Wei (Reference Keh and Wei2000). Quantitative agreement of both (a) the electrophoretic component 
$\varTheta_1 (\kappa)$ and (b) the chemiphoretic component 
$\varTheta_2 (\kappa)$ is observed.
This section highlighted the generality of (2.12) and (2.13). We showed they are able to recover the mobilities for microswimmers in thin interaction layer limit, electrophoresis for arbitrary double-layer thickness and electrolytic diffusiophoresis for arbitrary double-layer thickness.
4. Autophoretic motion of microswimmers
In this section, we use the formula obtained in (2.12)–(2.13) to study the translation of Janus-like particles with a spherical cap, see figure 4(a). The key novelty of our analysis is that (2.12)–(2.13) do not impose the restriction on interaction length scale. As we show later, if the interaction length is comparable to particle size, the particle velocity is significantly impacted.
(a) Self-phoretic Janus particle where propulsion is controlled by the size of the spherical cap 
$\theta _0$, the reactive flux 
$J$ and the interaction length scale 
$\kappa ^{-1}$. (b) Here, 
$U$ versus 
$\theta _0$ for different 
$\kappa$ values demonstrates a maximum velocity for 
$\theta _0={{\rm \pi} }/{2}$ irrespective of 
$\kappa$. (c) Here, 
${U\kappa ^2}/{U_0}$ versus 
$\kappa$ asymptotically approaches the thin limit calculations as 
$\kappa \to 100$. However, considerable dampening is observed even for 
$\kappa = O(10)$. The values are reported for 
$\theta _0={{\rm \pi} }/{2}$. (d) The dampening of 
${U\kappa ^2}/{U_0}$ with 
$\kappa$ is observed for all 
$\theta _0$. The value 
$U_0$ is the asymptotic limit of 
$U \kappa ^2$ from the thin interaction layer calculations.
We define catalytic surface activity through a non-dimensional outward surface flux of strength, 
$J$ (scaled by a characteristic flux 
$D c_0 / a$, where 
$D$ is the diffusivity of the solute, 
$c_0$ is a reference concentration of the solute and 
$a$ is the particle radius). The non-dimensional interaction length is characterized by 
$\kappa ^{-1}$ (scaled by 
$a$). Lastly, the size of the catalytic cap is controlled by the polar angle 
$\theta _0$, such that 
$\theta _0=0$ indicates no catalytic cap on the particle and 
$\theta _0={{\rm \pi} }/{2}$ represents a hemispherical cap.
We assume a Helmholtz-like equation governs the interaction potential (
$\phi$, scaled by 
$k_B T$) with a constant surface potential (
$\phi _0$) and a far-field decay. We take this opportunity to highlight the choice of Helmholtz-like potential. While the potential is not representative of different surface interactions possible, it provides a convenient choice to explore the impact of 
$\kappa$ and thus has been chosen for this analysis. We note that our analysis can be easily extended to other interaction potentials provided that the integrals in (2.11) are convergent. For a detailed analysis of phoretic motion due to a general particle–solute interaction the reader is directed to Brady (Reference Brady2021).
To resolve particle translation for a given surface activity and interaction, we write
$$\begin{gather} \nabla^2 \phi = \kappa^2 \phi, \end{gather}$$
$$\begin{gather}\phi = \phi_0,\quad {\rm at}\ r=1, \end{gather}$$
$$\begin{gather}\phi \to 0,\quad {\rm at}\ r \to \infty. \end{gather}$$Solute transport is governed by diffusion and phoretic interactions with the following boundary conditions:
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\boldsymbol{\nabla} c + c \boldsymbol{\nabla} \phi \right) = 0, \end{gather}$$
$$\begin{gather}-\boldsymbol{n}\boldsymbol{\cdot}\left(\boldsymbol{\nabla} c + c\boldsymbol{\nabla}\phi\right) = J,\quad {\rm at}\ r=1, \end{gather}$$
$$\begin{gather}c\to0,\quad {\rm at}\ r\to\infty. \end{gather}$$
Our model problem is illustrated in figure 4(a). We solve the coupled equations (4.1)–(4.6) numerically. A scaled model geometry was constructed with the particle radius given to be 
$a=1$ and an outer radius of 
$r|_\infty =20$, representing the far-field. The interaction potential (
$\phi$) and solute concentration (
$c$) were defined by (4.1)–(4.6). To obtain the translation velocity through (2.12), we define the body force 
$\boldsymbol {b} = -c\boldsymbol {\nabla }\phi -\boldsymbol {\nabla } c$. Note that the arbitrary body force has both a phoretic and an osmotic contribution. Details of the computational method are provided in Appendix B.
We first analyse the effects of the spherical cap size (
$\theta _0$) and interaction length scale (
$\kappa ^{-1}$) for a fixed surface flux (
$J=1$) and surface potential (
$\phi _0=-1$), see figure 4. First, we note that the particle moves in the direction of the catalytic cap, since when 
$\phi _0=-1$ the particle is attracted towards regions of higher solute concentration. From figure 4(b), we observe that the propulsion velocity is maximum for 
$\theta _0={{\rm \pi} }/{2}$. This is in agreement with prior observations in the literature (Golestanian, Liverpool & Ajdari (Reference Golestanian, Liverpool and Ajdari2007), Michelin & Lauga (Reference Michelin and Lauga2014), Popescu et al. (Reference Popescu, Uspal, Eskandari, Tasinkevych and Dietrich2018) and others). Additionally, the dependence of translation velocity on the cap size is symmetric about 
$\theta _0={{\rm \pi} }/{2}$. As described in Michelin & Lauga (Reference Michelin and Lauga2014), when 
$\theta _0={{\rm \pi} }/{2}$, for small Péclet numbers, the sharpest concentration gradients are located near the equator, consequently leading to a larger slip velocity over an extended inert surface. In contrast, for smaller or larger catalytic caps, the aforementioned solute front is closer to the pole and thus involves a smaller share of the particle surface in generating slip velocities and therefore a smaller swimming speed.
Thin interaction layer: we also study the dependence of propulsion velocity on the interaction length scale (
$\kappa ^{-1}$), see figure 4(c). To compare our calculations with the thin interaction layer limit, we perform analogous calculations following the approach of Anderson (Reference Anderson1989) and Derjaguin et al. (Reference Derjaguin, Sidorenkov, Zubashchenkov and Kiseleva1947). We consider diffusive transport of solute 
$\nabla ^2 c = 0$ through the fluid volume. The surface flux condition is given by 
$-\boldsymbol {n}\boldsymbol{\cdot}\boldsymbol {\nabla }c=J$ with a far-field decay condition. The phoretic slip at the particle surface is defined as
\begin{equation} \boldsymbol{u}_{slip} ={-} \boldsymbol{\nabla}_s c \displaystyle \int_1^\infty (r-1)\left[\exp\left(-\phi\right)-1\right] {\rm d}r, \end{equation}
where the non-dimensional 
$\boldsymbol {u}_{slip}$ is scaled with 
$k T c_0 a / \mu$. After that, we refer to Stone & Samuel (Reference Stone and Samuel1996), also derived in (3.7), to obtain the translation velocity 
$\boldsymbol {U}_{thin}$ to be
\begin{equation} \boldsymbol{U}_{thin} ={-}\frac{1}{4{\rm \pi}} \displaystyle \int_S \boldsymbol{u}_{slip} \,{\rm d}S. \end{equation}
It is well known in the diffusiophoretic literature (Golestanian Reference Golestanian2019) that for 
$\kappa \gg 1$, 
$|\boldsymbol {U}_{thin}| \propto 1/ \kappa _{thin}^2$. Consequently, we introduce a scaled velocity expression, 
$\boldsymbol {U}_0 = \boldsymbol {U}_{thin} \kappa ^2_{thin}$ which becomes constant as 
$\kappa \to \infty$. As shown in figure 4(c), we observe that the velocity ratio 
$U\kappa ^2/U_0 = |\boldsymbol {U}|\kappa ^2/|\boldsymbol {U}_0|$ significantly decreases when the interaction limit becomes comparable to particle size, i.e. 
$\kappa \rightarrow 1$. We find that even for 
$\kappa = O(10)$, the velocity is reduced by almost a factor of two. We only reach the thin interaction limit for 
$\kappa = O(10^2)$. This observation is consistent with the passive diffusiophoretic literature where even moderately thin double layers can significantly reduce the diffusiophoretic velocity (Prieve et al. Reference Prieve, Anderson, Ebel and Lowell1984), see figure 3(b). Our analysis highlights that even for autophoretic swimmers, this effect can be observed, and using a thin interaction limit could overestimate the velocity for moderately thin interaction thickness such as 
$\kappa \rightarrow 50$. Over the past decade, there has been an increasing interest in nanoparticles or very dilute systems where relative interaction length scales are large (Leunissen et al. Reference Leunissen, Van Blaaderen, Hollingsworth, Sullivan and Chaikin2007; Shin et al. Reference Shin, Um, Sabass, Ault, Rahimi, Warren and Stone2016; Gupta et al. Reference Gupta, Shim and Stone2020b; Wu, Greydanus & Schwartz Reference Wu, Greydanus and Schwartz2021; Shi, Wu & Schwartz Reference Shi, Wu and Schwartz2023). Thus, these results could be crucial for future experimental studies. To ensure that our trends of reductions in velocity are consistent for other conditions, in figure 4(d), we observe that reduction in velocity for smaller 
$\kappa$ values is consistent for all 
$\theta _0$ values.
5. Autophoretic swimmers with external solute gradients
To further demonstrate the applicability of (2.12) in scenarios where it is difficult to use a slip velocity approach, we modify our model problem to include an external solute flux in the fluid bulk. The model geometry is preserved as shown in figure 5(a). The concentration in the far-field is 
$\boldsymbol {\nabla } c_\infty = J_{ext} \boldsymbol {e}_z$, where 
$\boldsymbol {e}_z$ is the 
$z$-coordinate basis vector in the universal Cartesian frame of reference and 
$J_{ext}$ is a free parameter used to control the direction and strength of this external field. The magnitude of the external gradient is thus 
$\left |{\boldsymbol{\nabla}}c_{\infty} \right | = J_{ext}$. We note that all additional parameters have been appropriately non-dimensionalized as per the discussion in § 4. We analyse the variations of the local surface flux (
$J$), external solute flux (
$J_{ext}$) and interaction length scale 
$\kappa$. For a given surface flux 
$J=1$, it is observed that the translation velocity is linear with the imposed solute flux, 
$\boldsymbol {\nabla } c_\infty$. The particle moves with a velocity 
$\boldsymbol {U} = U_z \boldsymbol {e}_z$. The direction of the propulsion is governed by the relative magnitudes of different contributions of the body force terms. In the absence of an external flux, the body force contributions due to the phoretic activity cause the particle to move in the positive 
$z$-direction. In our analysis, we have ignored the excess osmotic contribution (
$\boldsymbol {\nabla } c$) from the external solute gradient and only considered the phoretic contribution. This is done to ignore the motion of the particle only due to 
$\boldsymbol {\nabla } c_\infty$ when 
$\phi =0$. For 
$J_{ext}>0$, one can observe that force terms arising from the gradients in active and passive concentration fields are in the same direction, hence 
$U_z >0$ as in figure 5(b). Alternatively when 
$J_{ext}\,{<}\,0$, the phoretic force terms arising out of the passive solute concentration field compete with forcing arising due to activity. This leads to a direction reversal below a threshold concentration gradient, 
$(\boldsymbol {\nabla } c_\infty )_{crit}$, where 
$U_z < 0$. Here, 
$(\boldsymbol {\nabla } c_\infty )_{crit}$ is the external solute gradient necessary to arrest particle motion due to activity. The linearity of the results in figure 5(b) is due to the computations being performed in the weak field limit. To inspect the effects of interaction potential length, we obtained the 
$(\boldsymbol {\nabla } c_\infty )_{crit}$ values at 
$J=1$ while varying over 
$\kappa$. A qualitative agreement is observed between the results in figure 5(c) and the 
$\kappa$-dependence of velocity (
$U$). Further, we see that a stronger external flux is necessary to arrest motion as we approach the thin interaction length limit. This effect arises because, for 
$\kappa \to \infty$, the catalytically ejected solute decays more sharply, 
$\boldsymbol {\nabla } c$ increases locally near the particle surface. This necessitates the need for a larger 
$(\boldsymbol {\nabla } c_\infty )_{crit}$ to counteract activity-induced gradients to arrest propulsion.
(a) A self-phoretic particle also driven by external 
$\boldsymbol {\nabla } c_\infty$. (b) Here, 
$U_z$ versus 
$\boldsymbol {\nabla } c_\infty$ shows that both the particle speed and direction depend on the competing effects of external and self-propulsion modes. We define 
$(\boldsymbol {\nabla } c_\infty )_{crit}$ as the value when the particle motion was arrested despite the presence of concentration gradients due to surface activities. (c) Here, 
$(\boldsymbol {\nabla } c_\infty )_{crit}$ versus 
$\kappa$. We observe an increase in the magnitude of the external flux needed to arrest motion as we approach thin interaction length limits.
6. Conclusion
The main result of our work is that (2.12)–(2.13), reported in the prior literature for phoretic systems (Brady Reference Brady2021) and different physical set-ups (Brenner & Cox Reference Brenner and Cox1963; Hinch Reference Hinch1991; Leal Reference Leal2007; Datt et al. Reference Datt, Zhu, Elfring and Pak2015; Elfring & Goyal Reference Elfring and Goyal2016; Datt et al. Reference Datt, Natale, Hatzikiriakos and Elfring2017), can retrieve the mobility expressions for both electrophoretic and electrolytic diffusiophoretic motion of colloidal spherical particles at arbitrary interaction lengths at the low potential limit. Additionally, the asymptotic limit of (2.12)–(2.13) for thin interaction length recovers net translation and rotation of the particle in terms of a phoretic slip velocity as in Stone & Samuel (Reference Stone and Samuel1996). Finally, we employ these mobility relationships for self-propulsion of spherical microswimmers where we observe a peak in translational velocity as we approach a hemispherical catalytic coverage with a velocity reduction for lower or higher coverages. At moderate interaction lengths, 
$\kappa = O(10)$, a dampening in the translation velocity is observed. However, for higher 
$\kappa$ values our calculations retrieve the thin interaction length limit results arising from equivalent slip velocity calculations. Further, we add an external passive solute concentration gradient to our problem to understand the competing effects of surface-generated and externally imposed solute concentration gradients. We find configurations where the external flux arrest motion induced by surface activity 
$(\boldsymbol {\nabla } c_\infty )_{crit}$ and regions where the propulsion induced by the passive external solute gradient aids with or competes against the propulsion induced due to surface activity. A nonlinear decay of 
$(\boldsymbol {\nabla } c_\infty )_{crit}$ with 
$\kappa$ is observed from our analysis which is qualitatively similar to the 
$\kappa$-dependence of the translation velocity, 
$|\boldsymbol {U}|$.
Beyond the result described in this paper, our work can be utilized to predict 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ for an arbitrarily shaped particle. For such a calculation, one would need the appropriate expressions of the mobility (
$\boldsymbol{\mathsf{M}}$) and the disturbance tensors (
$\boldsymbol{\mathsf{D}}_T$, 
$\boldsymbol{\mathsf{D}}_R$), which might be possible to obtain analytically or numerically. While this manuscript focuses on the mobility of a single particle, a similar analysis could be extended to multiple particles and particles under confinement.
The derived expression will be particularly useful for multiphysics propulsion. For instance, it might be possible to induce propulsion of particles using multiple modes such as a combination of electric fields and other fields, such as electrodiffusiophoresis (Wang, Behdani & Silvera Batista Reference Wang, Behdani and Silvera Batista2022; Jarvey, Henrique & Gupta Reference Jarvey, Henrique and Gupta2023). Another possibility is the inclusion of chemical kinetics at the particle surface (Davis & Yariv Reference Davis and Yariv2022) which will modify the solute problem and thus consequently change 
$\boldsymbol {b}$.
Acknowledgements
The authors acknowledge H.A. Stone, A.S. Khair, F. Henrique, G. Roure and the three anonymous referees for their input in the preparation of this manuscript. A.G. acknowledges the Teets Family Endowment in Nano-Technology Graduate Fellowship and the Mukhopadhyay Graduate Research Fellowship. S.R. acknowledges the Mukhopadhyay Graduate Research Fellowship. The authors acknowledge the donors of the American Chemical Society Petroleum Research Fund.
Funding
A.G. thanks the National Science Foundation (CBET-2238412) CAREER award for financial support.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical resolution of the radial components at 
$O(\epsilon \zeta )$ and $O(\epsilon \zeta ^2)$
Order 
$O(\epsilon \zeta )$: the radial components of (3.81)–(3.85) are
$$\begin{gather} \frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2\frac{{\rm d}f_{\rho_{11}}}{{\rm d}r}\right) - \left(\frac{2}{r^2} + \kappa^2\right) f_{\rho_{11}} + \frac{2}{r^2} \frac{{\rm d}}{{\rm d}r} \left(r^2 f_{s_{10}}(r) \frac{{\rm d}\phi_{01}(r)}{{\rm d}r}\right) =0, \end{gather}$$
$$\begin{gather}\frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2 \frac{{\rm d}f_{s_{11}}}{{\rm d}r}\right)- \frac{2f_{s_{11}}}{r^2}+\frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2\rho_{01}(r)\frac{ {\rm d}f_{\phi_{10}}(r)}{{\rm d}r}\right) - \frac{2 \rho_{01}(r) f_{\phi_{10}}(r)}{r^2}=0, \end{gather}$$
$$\begin{gather}\frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2\frac{{\rm d}f_{\phi_{11}}}{{\rm d}r}\right)-\frac{2f_{\phi_{11}}}{r^2} ={-}\frac{\kappa^2}{2}f_{\rho_{11}}, \end{gather}$$
where 
$f_{\phi _{10}}(r)$ and 
$f_{s_{10}}(r)$ are the radial components of the 
$O(\epsilon )$ solutions. The appropriate boundary conditions at the particle surface are
$$\begin{gather} \frac{{\rm d}f_{\rho_{11}}}{{\rm d}r} + f_{s_{10}}(r) \frac{{\rm d}\phi_{01}(r)}{{\rm d}r}=0, \end{gather}$$
$$\begin{gather}\frac{{\rm d}f_{s_{11}}}{{\rm d}r} + \rho_{01}(r)\frac{{\rm d}f_{\phi_{10}}(r)}{{\rm d}r}=0, \end{gather}$$
$$\begin{gather}\frac{\mathrm{d}f_{\phi_{11}}}{\mathrm{d}r} = 0. \end{gather}$$
In the far-field, the radial flux of 
$f_{\phi _{11}}$, 
$f_{\rho _{11}}$ and 
$f_{s_{11}}$ all go to zero.
Order 
$O(\epsilon \zeta ^2)$: the radial dependence of 
$s_{12}(r, \theta )$ is captured by
\begin{align} &\frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2\frac{{\rm d}f_{s_{12}}}{{\rm d}r}\right)- \frac{2f_{s_{12}}}{r^2}+\frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2\left\{\rho_{01}(r) \frac{{\rm d}f_{\phi_{11}}(r)}{{\rm d}r} + f_{\rho_{11}}(r) \frac{{\rm d}f_{\phi_{01}}(r)}{{\rm d}r}\right\}\right) \nonumber\\ &\quad -\frac{2 f_{\phi_{11}}\rho_{01}}{r^2} = 0, \end{align}with the boundary condition at the particle surface being
\begin{equation} \frac{{\rm d}f_{s_{12}}}{{\rm d}r}+\rho_{01} \frac{{\rm d}f_{\phi_{11}}}{{\rm d}r} + f_{\rho_{11}} \frac{{\rm d}\phi_{01}}{{\rm d} r}=0, \end{equation}
and 
$f_{s_{12}}=0$ in the far-field. We solve (A1)–(A7) using the bvp4c() function in MATLAB for 
$\kappa \in [1, 1000]$. For each 
$\kappa$ value we solve for 
$\rho$, 
$\phi$ and 
$s$ at each order with a one-dimensional mesh, 
$r = [1, 1+100/\kappa ]$ with a thousand elements. A default relative tolerance of 
$10^{-3}$ is used. To match with (3.93) we multiply the terms proportional to 
$\zeta$ in (3.92) with a factor of 
$1/2$ and the terms proportional to 
$\zeta ^2$ with a factor of four. This is solely due to how the coefficients 
$\varTheta _1(\kappa )$ and 
$\varTheta _2(\kappa )$ are defined in (3.93).
Appendix B. Numerical solution to the autophoretic motion of microswimmers in §§ 4 and 5
Section 4: the numerical solutions to (4.1) to (4.6) were obtained from the finite element method software COMSOL. We use a non-dimensional spherical computational domain of size 
$4/3{\rm \pi} \times 20^3$ where the particle is located at the origin and possesses a radius of unity. The domain is discretized into approximately 1 162 711 elements. Around the particle, we mesh a boundary region with 12 layers and a stretching factor of 1.1, consisting of approximately 20 000 triangular elements. The simulations are performed in the reference frame of the swimmer. Upon solving for the solute concentration 
$c$, and 
$\phi$ in the domain we define 
$\boldsymbol {b} = -c\boldsymbol {\nabla }\phi -\boldsymbol {\nabla }c$ as the body force. Subsequently, (2.12) is numerically integrated over the domain to find the translation velocity 
$\boldsymbol {U}$. All results are normalized with 
$\kappa ^2$ as illustrated in figures 4 and 5.
The translation velocity at the thin double layer limit is obtained by a similar procedure by solving for 
$c$ and using (4.7) to obtain 
$\boldsymbol {u}_{slip}$. For the second simulation set-up, to calculate the velocity at the thin limit an extremely fine mesh is required to discretize the domain, with approximately 1 154 140 elements. The lumped phoretic mobility for a given surface potential 
$\phi$ is obtained by numerically integrating along the radial direction in MATLAB from 
$1$ to approximately 
$10^4$ along the radial direction. The translation velocity 
$\boldsymbol {U}_{thin}$ is obtained by solving (4.8) in COMSOL.
Section 5: to solve the results where the swimmer was subjected to an external concentration gradient, we modify the far-field boundary condition in the above numerical simulation to be 
$c_\infty = c_0 + J_{ext}z$. A Dirichlet boundary condition was used instead of a Neumann boundary condition to avoid under-specifying our model. The Dirichlet condition will satisfactorily approximate 
$\boldsymbol {\nabla }c_\infty = J_{ext}\boldsymbol {e}_z$ in the vicinity of the particle for a large enough computational domain.
Open access
