Phoretic self-propulsion of a slightly inhomogeneous disc

Abstract In two dimensions, the problem governing a homogeneous phoretic swimmer of circular cross-section is ill-posed because of the logarithmic divergence associated with a purely diffusive solute transport. We address here the well-posed problem that is devised by introducing a slight inhomogeneity in the interfacial chemical activity. With the radial symmetry being perturbed, phoretic motion is animated by diffusio-osmosis. Solute advection, associated with that motion, becomes comparable to diffusion at large distances. The singular problem associated with that scale disparity is analysed using matched asymptotic expansions for arbitrary values of the Damköhler number $\textit {Da}$ and the intrinsic Péclet number $\textit {Pe}$. Asymptotic matching provides an implicit equation for the particle velocity in terms of these two parameters. The velocity exhibits a non-trivial dependence upon the sign $M$ of the slip coefficient. For $M=-1$, we observe the appearance of several solutions beyond a $\textit {Da}$-dependent critical value of $\textit {Pe}$. We also address the respective limits of small and large $\textit {Da}$ for fixed $\textit {Pe}$ and arbitrary inhomogeneity, and illuminate their linkage to the limit of weak inhomogeneity.


Introduction
Motivated by experimental realisation of catalytic swimmers, there is a surging interest in phoretic self-propulsion of colloidal particles (Aubret, Ramananarivo & Palacci 2017).The mechanism underlying phoretic self-propulsion includes two ingredients: production and consumption of solute molecules at the particle boundary, and short-range interactions between these molecules and the rigid boundary.A spatial asymmetry in the interfacial activity may lead to particle motion in a preferred direction.
The preceding discussion may give the impression that a non-zero Péclet number guarantees the presence of solute advection, which in two dimensions regularises the ill-posed problem.To see that this is not necessarily the case, it is expedient to consider a highly symmetric configuration, involving a spherical particle with a homogeneous interfacial activity.The governing equations in that configuration are clearly solved by a quiescent state, where the velocity field vanishes trivially and the particle is stationary.For any value of Pe, then, advection is absent in this quiescent solution.This observation can be traced back to the meaning of the Péclet number Pe.Since the velocity scale U associated with diffusio-osmosis is not externally imposed (as in forced convection), the Péclet number Pe is a native property of the system (see (2.4)).As such, it does not necessarily represent the ratio between advection and diffusion.
If the sphere is replaced by a cylindrical particle, then the associated two-dimensional problem is ill-posed -for any value of Pe.A possible regularisation is provided by a deviation from isotropy, say via a slightly inhomogeneous interfacial activity.Because of the deviation from radial symmetry, a flow field is triggered by diffusio-osmosis, resulting in turn in particle motion.Solute advection then renders the problem well-posed.If the degree of inhomogeneity is quantified by the small parameter , then it is evident that the problem is singular in the limit → 0.
Our goal is to address this singular limit, where, following the spirit of Yariv (2017), we allow for an arbitrary distribution of the deviation from homogeneity.The analysis of this problem involves several attractive features.The first is its distinctive generality: as will become evident, the problem can be solved in closed form for arbitrary values of Pe and Da -a desirable attribute in this class of self-propulsion problems (Boniface et al. 2019).As will also become evident, the ability to cover the entire range of these parameters results in the prediction of anomalous phenomena, such as velocity reversal beyond a critical value of Pe.
The second attractive feature is a mathematical one.Unlike the previous two-dimensional investigations (Yariv 2017;Yariv & Crowdy 2020), where the interfacial inhomogeneity was appreciable, the present problem requires going beyond leading-order analysis in the asymptotic paradigm.It therefore provides a vivid illustration of the way by which generalised Poincaré expansions allow for a systematic use of Van Dyke matching in problems involving an inherent logarithmic singularity (Van Dyke 1964;Hinch 1991).
The third attractive feature has to do with spontaneous motion of isotropically active particles, a phenomenon made possible by the nonlinear mechanism of solute advection.The symmetry breaking leading to spontaneous motion has been discussed extensively in the context of the above-mentioned three-dimensional highly symmetric problem (Michelin, Lauga & Bartolo 2013;Saha, Yariv & Schnitzer 2021), and the possibility of spontaneous motion in two dimensions is of interest (Hu et al. 2019).The linkage of the present problem to spontaneous motion in two dimensions is discussed in § 11.

Problem formulation
Consider a two-dimensional configuration where a circular disc of radius a is suspended freely in an unbounded liquid solution (solute diffusivity D).The reference solute concentration at large distances from the disc is denoted by c ∞ .The chemical activity at the particle boundary is modelled as a first-order chemical reaction (Ebbens et al. 2012 where the (positive) rate constant k generally varies along the circumference.We assume that the distribution of k is symmetric about a diameter of the circle.The circumferential average of k is denoted by k.
We employ a macroscale description, where the short-range interaction between the solute molecules and the particle is expressed by an effective diffusio-osmotic slip (Anderson 1989).Thus in a particle-fixed system, the standard no-slip condition is replaced by slip velocity = b × surface gradient of solute concentration. (2.2) Following common practice (Michelin & Lauga 2014;Sondak et al. 2016), we assume that the slip coefficient b is uniform.Note that b is a signed quantity, positive for repulsive interactions and negative for attractive ones.The velocity scale associated with (2.2) is (2.4) It follows from the problem symmetry that the force-and torque-free cylinder reacts to the interfacial slip by translating along its symmetry diameter at a constant velocity, without rotating.The translation velocity, denoted by s, is reckoned positive when the particle propagates in the direction of its more active hemicircle.(This definition will be made precise later on; see (3.12).)Our goal is the determination of this velocity.

Dimensionless formulation
In what follows, we normalise all length variables by a.The analysis is carried out in a particle-fixed reference system with origin at the disc centre.In that system, we use (x, y) Cartesian coordinates defined such that the x-axis is aligned along the symmetry diameter of the particle, pointing in the direction of the more active hemicircle.Additionally, we utilise (r, θ) polar coordinates, with θ measured in the counterclockwise direction from the x-axis.Accordingly, the distribution of the rate constant k is denoted by k(θ ).The dimensionless geometry is portrayed in figure 1.
We employ a dimensionless notation where the solute concentration is normalised by c ∞ and the velocity field is normalised by U. Our interest is in the velocity U(= s/U ) of the particle relative to the otherwise quiescent liquid; in the particle-fixed reference frame, this velocity is manifested as the uniform streaming −Uî at infinity, î being a unit vector in the x-direction.
To obtain this velocity, we consider the coupled problem governing the solute-concentration distribution c and velocity field u.The solute-transport problem is governed by: (i) the advection-diffusion equation, (3.1) (ii) the approach to the reference concentration at large distances, and (iii) the kinetic condition at the particle boundary, where is the dimensionless distribution of rate constant, and is the Damköhler number, representing the ratio of diffusive (a 2 /D) to reactive (a/ k) time scales.
The velocity field u = êr u + êθ v is governed by: (i) the continuity and Stokes equations (the former tacitly employed in (3.1)); (ii) the impermeability condition where M = b/|b|; (iv) the far-field approach to a uniform stream, lim r→∞ u = −Uî; (3.8) and (v) the requirement that the particle is force-free.The latter, in conjunction with (3.7) and (3.8), provides the particle velocity as a quadrature (Squires & Bazant 2006), The inhomogeneous circumferential reaction is described in condition (3.3) by the distribution F(θ ), which is required to be strictly positive: (3.10) By definition of k it follows that F satisfies the integral constraint Also, recalling the sense of direction of the x-axis, we impose the constraint (3.12) Now, the symmetry of the k distribution about the x-axis implies that F is even, i.e.
so that (3.11) and (3.12) may be written as With activity higher on the right hemicircle, one expects the solute concentration about that hemicircle to be lower than the corresponding concentration about the left hemicircle.The associated solute gradients are then expected, in some average sense, to point from right to left.For M = 1, this would apply slip at the same direction, whence particle motion to the right, whereby U > 0 (see indeed (3.9)).For M = −1, one would then expect U < 0. The preceding crude arguments may seem to suggest that MU > 0. (3.15) A stronger version of the odd dependence (3.15) holds when U is proportional to M: This version of (3.15) holds in the linear problem that applies (in three dimensions) in the absence of solute advection.It was also observed in a small-Pe analysis in two dimensions (Yariv & Crowdy 2020).Since we did not prove property (3.15), it remains questionable for the time being.

Simplifications
The nonlinear problem formulated in § 3 clearly admits the symmetries It is therefore sufficient to analyse the problem in the domain 0 < θ < π.

Phoretic self-propulsion of a slightly inhomogeneous disc
A further simplification follows by making use of a streamfunction ψ, defined by The continuity equation is then satisfied trivially, while the Stokes equation governing u is transformed to the biharmonic equation governing ψ: The impermeability condition (3.6) implies that ψ is a constant at r = 1.Since ψ is defined to within an additive constant, we simply write In terms of ψ, the slip and far-field conditions (3.7) and (3.8) read, respectively while the symmetries (4.1a-c) imply that Finally, we note that in term of ψ, the advection-diffusion equation (3.1) becomes (4.8)

Weak inhomogeneity
We write where > 0. (5.2) From (3.13), (3.14a,b) and (5.2), we find that f satisfies the respective properties In addition, we impose the constraint (without any loss of generality) max (5.7) In terms of and f , the absorption condition (3.3) becomes With the norm of f set to ord(1) by (5.4), the parameter may be interpreted as an effective measure of interfacial inhomogeneity.To capture the case of weak inhomogeneity, we consider the asymptotic limit 0. We expand all pertinent field variables in the generic form (5.9) This results in the respective expansion of U: (5.10)

Breakdown of the weak-advection approximation
At ord(1), condition (5.8) becomes radially symmetric: In the absence of flow at ord(1), the right-hand side of (4.8) vanishes at this order.It then follows that c 0 is governed by Laplace's equation We anticipate that the solution of (6.1) and (6.2) is radially symmetric for all r.Conversely, with a radially symmetric c 0 , no flow is animated by the slip condition (3.7).We therefore proceed assuming that c 0 is a function of r alone, (6.3) and that u 0 (and hence ψ 0 ) vanish trivially.Consequently, U 0 = 0. (6.4) Substituting (6.3) into (6.1) and ( 6.2), we obtain c 0 = a 0 (Da ln r + 1).(6.5) The integration constant a 0 cannot be determined from the ord(1) balance of the far-field condition (3.2), lim which is incompatible with (6.5).
The above incompatibility is associated with the breakdown of the dominance of diffusion, associated via (6.2) with the small-limit.Indeed, the solution (6.5) implies that the presumably ord( ) advective term in (3.1) becomes comparable to the diffusive term at distances r = ord( −1 ).As a consequence of the above non-uniformity, we need Phoretic self-propulsion of a slightly inhomogeneous disc to employ two separate asymptotic expansions of c.The first, holding at ord(1) distances from the disc (the 'nearby' region), is provided by the generalisation of (5.9) where we allow for a weak (i.e.logarithmic) dependence upon at each order.The second expansion holds at the 'remote' region, where r = ord( −1 ).Towards that end, we employ the stretched radial coordinate r = r, (6.8) which is ord(1) in the remote region.Employing the change of variables c(r, θ; ) = c (r , θ; ), (6.9) the remote region expansion is where again we allow for a logarithmic dependence upon at each asymptotic order.The use of generalised Poincaré expansions, where the terms are allowed to depend upon ln , enables a straightforward use of the Van Dyke matching rule, where logarithmically large terms are considered on a par with ord(1) terms (Hinch 1991).Two points are worth noting: First, the far-field decay (3.2) now applies on c : lim No far-field condition applies on the nearby field c, which describes the concentration only for r = ord(1).Thus the distribution (6.5), obtained prior to the identification of two length scales, remains valid (with a 0 being allowed to depend weakly upon ).Since (3.9) clearly applies in the nearby region, (6.4) is also retained.Instead of the far-field condition (3.2), the nearby concentration c must satisfy appropriate large-r conditions that follow from the requirement of asymptotic matching with the remote concentration.These conditions determine uniquely the integration constants appearing in the various {c i } (e.g. a 0 in (6.5)).
Second, the velocity field is affected by the solute concentration only through the slip condition (3.7), which applies at r = 1.Accordingly, there is no need for a separate expansion of the velocity field, akin to (6.10), at the remote region.The only modification required in the expansion of the velocity field is the transition from (5.9) to (6.7).
Defining the remote-region gradient operator ∇ = −1 ∇, solute transport in that region is governed by (cf.(3.1)) wherein u is to be evaluated at r = ord( −1 ).

Leading-order remote solution
Given condition (3.8) and observation (6.4), the velocity in the remote region is simply given at leading order by constant vector − î U 1 .Consequently, the leading-order balance 940 A24-9 Writing we find that G, which decays at large r , also satisfies (7.1).Substituting we find that H satisfies the modified Helmholtz equation The solution of that equation that decays at large r and is least singular at r = 0 is a radially symmetric screened source, K 0 (Pe |U 1 | r /2), in which K 0 is the modified Bessel function of the second kind.While all spatial derivatives of that source also satisfy (7.5), they are too singular at the origin and are excluded by the requirement of asymptotic matching with (6.5).We therefore obtain where the source magnitude Q remains to be determined.Using the small-argument behaviour of K 0 (Abramowitz & Stegun 1965), we find wherein γ is the Euler-Mascheroni constant.The asymptotic error in (7.7) is 'algebraically small' in r , i.e. smaller than some positive power of r .We now perform Van Dyke ord(1) : ord(1) matching, using (6.5).Comparison of the ln r terms gives Comparison of the r-independent terms then yields We have therefore evaluated c 0 in terms of the (as yet unknown) leading-order particle velocity.To obtain that velocity, we need to consider deviations from radial symmetry in the nearby region.Towards formulating the far-field condition governing c 1 , we expand c 0 at small r up to ord(r ), thus obtaining the following refinement of (7.7): where we have made use of (7.8).It is seen readily that the asymptotic error in (8.3) is 'algebraically smaller' than r , meaning that its ratio to r is algebraically small in r .Making use of Van Dyke ord(1) : ord( ) matching, we obtain the following condition on c 1 : It is evident from both the far-field behaviour (8.4) and the forcing terms in (8.1) and (8.2) that c 1 lacks radial symmetry.The integral expression (3.9) then implies a non-zero U 1 : consistently with the analysis of solute transport in the remote region.Our first objective is to determine the general form of ψ 1 .This field is governed by: (i) the biharmonic equation (see (4.Given the general structure of single-valued solutions of the biharmonic equation in polar coordinates (Leal 2007), we find that the force-free solution that satisfies (8.6)-(8.8) is In principle, the velocity U 1 and the set {α (n) } are determined from the slip condition (8.9).We therefore head to the calculation of c 1 .Substituting (6.5) and (8.10) into (8.1),we obtain the Poisson equation Also, substitution of (6.5) into (8.2) yields the condition (8.12)In addition to (8.11) and (8.12), c 1 must also satisfy the matching condition (8.4).We note that condition (5.3a) implies a cosine series representation for f , say (8.13) wherein are the associated Fourier coefficients.We write c 1 as a sum of corresponding Fourier modes: Substituting into (8.5) and making use of the orthogonality of the set {cos nθ} ∞ n=0 , we obtain Given the linear problem governing c 1 , we can exploit that orthogonality further and seek only c (1) (r).This function satisfies the inhomogeneous equation (cf.(8.11)) together with the matching condition (cf.(8.4)) (8.18) and the boundary condition (cf.(8.12)) dc (1)  dr − Da c (1) = Da a 0 f (1)

Phoretic self-propulsion of a slightly inhomogeneous disc
The most general solution of (8.17) and (8.18) is (8.20) The coefficient N 1 is obtained readily from (8.19).It then follows that Substitution into (8.16)eventually gives, upon making use of (7.9) and (8.14), The implicit equation ( 8.22) provides U 1 as a function of the governing four parameters.
The appearance of ln results in a non-analytic dependence upon at = 0, as was to be expected from the singular nature of the limit → 0.

The limits of slow and fast kinetics
By replacing U 1 with U/ in (8.22) and making use of (5.1), we obtain an implicit equation for U itself: with an asymptotic error that is algebraically small in .Note that the dependence upon has disappeared; this was to be expected, as this parameter is not required for the problem specification.
When Da becomes small, we find from (9.1) the approximation which is no longer singular at U = 0.In the other limit, as Da → ∞, we obtain which essentially implies an ord(Da −1 ) scaling of U. To understand the above limits, it is expedient to consider the limits Da 1 and Da 1 from the outset, allowing for to attain ord(1) values (subject to both (5.2) and (5.5)).It therefore proves useful to employ the original problem formulation, in terms of F(θ ).
We start with the small-Da limit.With interfacial reaction being ord(Da) weak, it results in an ord(Da) perturbation to the unity equilibrium concentration specified by (3.2).The slip condition (3.7) then implies that the velocity field is ord(Da) small too, as is U itself.Writing c = 1 + Da c, we find that at leading order c satisfies Laplace's equation in the vicinity of the disc, together with the Neumann condition (cf.(3.3)) The constant ã cannot be determined from the analysis in the particle vicinity.Rather, it may be obtained by asymptotic matching with the leading-order approximation for c in a remote region, which is now at distances ord(Da −1 ).(This matching procedure excludes terms diverging like r and faster in (9.5).)However, this constant is irrelevant to the calculation of the particle velocity.Thus, writing U = Da Ũ, we obtain readily from (3.9) that in agreement with (9.2).With the remote region being irrelevant, it is hardly surprising that approximation (9.2) is independent of Pe.While the original problem is ill-posed for Da = 0, it appears as though the singularity in the limit Da → 0 is a 'removable' one.With c 0 uniform at r = 1, there is no slip to drive ord(1) fluid velocity, so we anticipate that u = ord(Da −1 ).The leading-order concentration then satisfies Laplace's equation.
Seeking the solution that diverges least rapidly at large r, we obtain c 0 = a 0 ln r. (9.9) The constant a 0 is determined from matching with a remote-region solution, now at distances Da −1 .Writing U = Da −1 U 1 + • • • and repeating an analysis similar to that in § 7 we find that (cf.(7.9)) At ord(1), we now find from (3.3) that Substituting into (3.9)then gives Introducing f via (5.1) and forming the small-limit yields Upon using (9.10) and the Fourier expansion (8.13)-(8.14),we finally obtain Reverting to U and F, we then retrieve (9.3).
The preceding analyses clarify why U becomes small for both small and large Da.For small Da, the interfacial reaction is so slow that c deviates only by ord(Da) values from the unity equilibrium values.These deviations result in an ord(Da) particle speed.For large Da, the reaction is so fast that c nearly vanishes at the interface, reducing there to ord(1/Da) values (up to logarithmic terms).Since it is the interfacial distribution that determines the particle speed (recall (3.9)), we again find that U is small, now essentially ord(1/Da).For fixed and Pe, both limits are associated with weak advection.

Dependence upon M
We now address the validity of (3.15) and (3.16).In the limits Da 1 and Da 1, we observe from (9.2) and (9.3) that U satisfies (3.16).This proportionality does not apply, however, in the general case: indeed, it is verified readily from (8.22) that U 1 is not an odd function of M. This is hardly surprising: a non-trivial dependence upon M was already observed by Michelin & Lauga (2014), who analysed self-propulsion in three dimensions at finite Pe.The dependence upon M in the present problem is illustrated in figure 2, where the product MU 1 is portrayed, for both M = 1 and M = −1, as a function of Da for Pe = 1 and = 0.1, using the distribution (5.7).
Having found that U does not satisfy (3.16), we now address the weaker version (3.15).With the logarithmic expression ln(4/ Pe |U 1 |) being asymptotically large and positive for 1, it is evident that U 1 > 0 for M = 1, in agreement with (3.15).For M = −1, the situation is more delicate.Assuming again that the logarithmic expression is large and positive, we find that U 1 < 0 for Pe that is smaller than 2 + 2Da.For Pe = 2 + 2Da, For Pe larger than 2 + 2Da, however, the term multiplying the logarithmic expression becomes negative, suggesting in turn a positive U 1 .
To understand the possibility for a velocity reversal, we assess the critical value of Pe at which the expression within the brackets appearing in the left-hand side of (8.22) would vanish approximately.Towards this end, we temporarily abandon our generalised asymptotic approach, where logarithmically large terms are considered on a par with ord(1) terms, and consider an expansion of that expression in inverse powers of ln(1/ ).In performing that expansion, we take the difference between the critical value Pe c and 2 + 2Da to be of order 1/ ln(1/ ).Thus writing with Δ being ord(1), we find that at ord(1) the expression within the brackets becomes The above analysis may appear to suggest a transition from negative to positive U 1 at a critical value at which the expression within the brackets appearing in the left-hand side of (8.22) vanishes, giving rise to a diverging U 1 .This is, however, only a crude estimate: as the expression within the brackets approaches zero, U 1 (which is negative) becomes large, and the above balance becomes invalidated due to the ln |U 1 | term in the brackets.The numerical solution of (8.22) for M = −1 actually reveals that a bifurcation occurs at a critical value of Pe, where two positive solutions appear.This is illustrated in figure 3, where we portray the velocity U 1 as a function of Pe for = 0.1 and Da = 1, using the distribution (5.7).At these values, the bifurcation appears at Pe ≈ 17.3.We conclude that (3.15) is invalid.

Concluding remarks
We have analysed the phoretic self-propulsion of a circular disc whose interfacial activity is slightly inhomogeneous.The analysis was based upon a linearisation scheme that exploited the slight inhomogeneity.This scheme implies that advection is relatively weak even for ord(1) values of the intrinsic Péclet number Pe.The dominance of diffusion over advection is not uniformly valid.The analysis of this singular limit requires the use of matched asymptotic expansions.As is customary in singular problems in two dimensions (Frankel & Acrivos 1968), we have employed generalised Poincaré expansions, allowing for a logarithmic dependence of the various asymptotic orders upon the expansion parameter .This has allowed for a straightforward use of the Van Dyke matching rule.
Our key result is (9.1), which provides, for a given distribution of interfacial reaction rate, the dimensionless velocity of the particle as an implicit function of the slip sign M, the intrinsic Péclet number Pe, and the Damköhler number Da.The functional dependence upon M is non-trivial, with a possibility of velocity reversal existing for M = −1.We have also obtained simpler approximations in the limits of small and large Da, which are not restricted by the assumption 1.The velocity reversal appearing for M = −1 beyond a critical value of Pe brings to mind the possibility of spontaneous motion at finite Pe.This possibility was investigated in detail in three dimensions, using the highly symmetric configuration discussed in § 1.Thus Michelin et al. (2013) observed that in addition to the 'trivial' quiescent solution, with a stationary particle, the governing equations also admit a non-trivial solution beyond a critical value of Pe; this solution corresponds to self-sustained motion.For an absorbing boundary (in tune with (3.3)), such spontaneous motion can take place only for M = −1.Recently, Saha et al. (2021) studied the linearised response of the same system to a weak external force and found a bifurcation at the same critical value.In fact, it was shown by Saha et al. (2021) that the study of the forced problem provides a convenient alternative to the stability analysis of the unforced problem.
Can spontaneous motion also occur in two dimensions?This is a subtle question: with the reference state being ill-posed (Hu et al. 2019), there is no two-dimensional counterpart to the stability analysis of Michelin et al. (2013).One can imagine investigating the weakly-forced problem, following Saha et al. (2021).In that problem, the advection triggered by the motion of the disc renders the solute-transport problem well-posed.Unfortunately, a forced problem cannot be studied under the framework of Stokes equations in two dimensions, due to the Stokes paradox.We propose that the present analysis, where the highly symmetric configuration is perturbed by a weak inhomogeneity, may provide a first step in the right direction.Indeed, the bifurcation of U 1 at a critical value of Pe, which approaches 2 + 2Da as → 0, may be a signature of spontaneous motion.
Finally, it is illuminating to provide here the dimensional counterpart of (9.1), namely associated with the particle speed.While our analysis is valid provided that the latter number is small, we anticipate that (11.1) may provide a rough approximation even when this restriction is violated.

Figure 1 .
Figure 1.Schematic of the dimensionless geometry, indicating the dimensionless differential equations and boundary conditions.
.3) Defining the intrinsic Péclet number Pe as a U/D thus gives Pe = |b| c ∞ D .

Figure 2 .
Figure 2. Variation of MU 1 with Da for Pe = 1 and = 0.1, using the distribution (5.7).The solid curves are obtained from the implicit equation (9.1).The dashed curves show the small-Da approximation (9.2) and the large-Da approximation (9.3).

940Figure 3 .
Figure3.Variation of U 1 with Pe for M = −1, Da = 1 and = 0.1, using the distribution (5.7).We show both the negative branch and the multi-valued positive branch, which exists for Pe > Pe c .