{'fraction_non_alphanumeric': 0.0383428208031518, 'fraction_numerical': 0.012579299219804253, 'mean_word_length': 4.695909532644227, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Weakly nonlinear dynamics of a chemically active particle near the threshold for spontaneous motion. I. Adjoint method

In this Series, we study the weakly nonlinear dynamics of chemically active particles near the threshold for spontaneous motion. In this part, we focus on steady solutions and develop an `adjoint method' for deriving the nonlinear amplitude equation governing the particle's velocity, first assuming the canonical model in the literature of an isotropic chemically active particle and then considering general perturbations about that model. As in previous works, the amplitude equation is obtained from a solvability condition on the inhomogeneous problem at second order of a particle-scale weakly nonlinear expansion, the formulation of that problem involving asymptotic matching with a leading-order solution in a remote region where advection and diffusion are balanced. We develop a generalised solvability condition based on a Fredholm Alternative argument, which entails identifying the adjoint linear operator at the threshold and calculating its kernel. This circumvents the apparent need in earlier theories to solve the second-order inhomogeneous problem, resulting in considerable simplification and adding insight by making it possible to treat a wide range of perturbation scenarios on a common basis. To illustrate our approach, we derive and solve amplitude equations for a number of perturbation scenarios (external force and torque fields, non-uniform surface properties, first-order surface kinetics and bulk absorption), demonstrating that sufficiently near the threshold weak perturbations can appreciably modify and enrich the landscape of steady solutions.


I. INTRODUCTION
Phoretic flows are surface-driven flows that result from local physico-chemical (chemical, thermal, electrical, etc.) gradients within a fluid phase adjacent to a surface [1]. Such gradients may be either externally imposed or self-generated by an 'active' surface. A prototypical example of the latter is the phoretic Janus particle, which self-propels as a consequence of its front and back sides having different properties [2][3][4][5][6][7][8]. Concomitantly, flow can modify physico-chemical gradients through advective transport. There is therefore mutual coupling between physico-chemical fields and liquid flow which under certain conditions results in unstable growth of perturbations. Together with the nonlinearity inherent to advection, such instabilities can lead to rich dynamics. A phenomenon that has received considerable attention in this context is symmetry breaking, where the physico-chemical and flow fields spontaneously form steady or unsteady structures possessing less symmetry than the governing equations. Examples include the formation of convection rolls in active channels [9][10][11], as well as the spontaneous dynamics exhibited by freely suspended and isotropic active particles [11][12][13][14][15] and drops [16][17][18][19][20][21][22][23].
A canonical model for studying the dynamics of an isotropic active particle was introduced by Michelin et al. [12]. It consists of a chemically active spherical particle that is freely suspended in an unbounded liquid solution. A single species of solute molecules is transported in the liquid bulk by diffusion and advection, approaching an equilibrium concentration at infinity. The chemical activity of the particle is represented by a prescribed flux of solute molecules at the boundary of the particle and flow is driven at that boundary by diffusioosmotic slip locally proportional to the surface gradient of the solute concentration. The prescribed solute flux and the slip coefficient are assumed uniform. This isotropic scenario allows for a steady state where both the fluid and particle are at rest and the solute distribution is spherically symmetric. It was shown by Michelin et al., however, that, depending on the signs of the solute flux and slip coefficient, the stationary-symmetric state can be linearly unstable beyond a critical 'intrinsic' Péclet number which quantifies the strength of advection relative to diffusion in the problem. At the threshold, the linear mode that becomes unstable describes steady rectilinear motion of the particle with arbitrary velocity and no rotation; the imaginary part of the growth rate vanishes there, implying a monotonic instability. Additional linear modes corresponding to higher wavenumbers become unstable at higher Péclet numbers, though these do not involve particle motion.
Michelin et al. [12] have also numerically simulated their canonical model as an initial value problem. These simulations show that, in an interval of Péclet numbers above the instability threshold, the particle approaches with time a steady state of spontaneous rectilinear motion in an arbitrary direction (determined in practice by initial conditions). As a function of the Péclet number, the speed of this spontaneous motion grows linearly away from the instability threshold, corresponding to an unconventional 'singular' pitchfork bifurcation [24]. More recent simulations [11,15] have revealed that at high Péclet numbers, at which the growth rates of higher-wavenumber linear modes dominate and the nonlinear spontaneous-motion states are expected to be unstable, the canonical particle exhibits complex and ultimately chaotic-like unsteady dynamics [11,15].
Several variations on the above canonical model of an isotropic active particle have been considered. These include variations to the chemical modeling that retain isotropy, such as accounting for first-order kinetics of the chemical reaction at the surface of the particle [13] or solute absorption in the bulk of the liquid solution [24]. Non-isotropic variations have also been considered, for example involving non-uniform surface properties [13,25], geometric confinement [26] or external force fields [27,28]. In the latter scenarios, there is a particular interest in how an imposed asymmetry influences the intrinsic spontaneous motion of the particle. The canonical model is also considered to be a 'reference model' for active drops, which despite being more complex often exhibit qualitatively (and, in certain limits, quantitatively) similar dynamics [14,23,29,30]. There have been many studies of the influence of the environment on the spontaneous dynamics of active drops, for example involving external force fields [29,31], surfactant transport [16], pair interactions [32,33], cluster dynamics [22], bi-motility associated with viscosity variations [20] and motion near boundaries [34]. Incidentally, we note that several related 'toy models' have also been proposed, involving physically inconsistent simplifications of either the canonical model of an isotropic active particle or active-drop models. These include two-and three-dimensional 'truncated' models where the concentration field (but not the flow field) is cut-off at some prescribed finite radius [21,[35][36][37], and point-particle models where advection in the vicinity of the particle is discarded [24,38] or included in an ad hoc manner [39].
Besides linear stability analysis and direct numerical simulations, several authors have applied weakly nonlinear analysis to describe various aspects of the spontaneous dynamics of chemically active particles and drops [14,16,21,24,[28][29][30][31][32]. Most of these works have focused on the limit where the bifurcation parameter, the Péclet number in the canonical model, approaches its threshold value for instability and spontaneous motion. In that limit, weakly nonlinear analysis generally leads to a nonlinear 'amplitude equation' governing the long-time dynamics, including the bifurcation of steady states, of the particle velocity vector. A useful feature of weakly nonlinear analysis is that it allows analytically treating perturbations [28,31,32], which no matter how weak can still have an appreciable effect sufficiently close to the instability threshold. Some systems, such as deformable drops, involve multiple control parameters such that two linear modes can be tuned to lose stability at nearby Péclet numbers [30,36,40]. In such scenarios, weakly nonlinear analysis allows deriving coupled amplitude equations governing the pair of modes [36,40].
For the most part, analysis near the onset of spontaneous motion of an active particle or drop follows the standard paradigm of weakly nonlinear analysis near the threshold of a monotonic instability. Specifically, the first-order terms in the weakly nonlinear expansion correspond to a linearisation about the basic stationary state at the instability threshold. At this order, one finds a homogeneous, linear and quasi-static problem which is singular, having a family of homogeneous solutions representing spontaneous rectilinear motion of the particle with arbitrary velocity. Then at a higher order one finds an inhomogeneous version of that singular problem, whose solvability yields the requisite amplitude equation. A non-standard feature is that the weakly nonlinear expansion can be spatially nonuniform, holding in the vicinity of the particle or drop but not at large distances where advection and diffusion are comparable. Following classical analyses of forced heat advection at small Péclet numbers [41], this spatial non-uniformity can be resolved using the method of matched asymptotic expansions [42], as done for steady states in [14,28,30]. This spatial non-uniformity of the asymptotics carries two important consequences. The first, explained by Farutin and Misbah [24], is the unusual singular-pitchfork bifurcation of the spontaneous-motion steady states; related to this, the amplitude equation arises from a solvability condition at second, rather than third, order of the weakly nonlinear expansion. The second consequence, on which we elaborate below, has to do with the nature of the unsteady dynamics near the instability threshold. For some models of active particles and drops the spatial non-uniformity of the weakly nonlinear expansion is absent. That is the case, for example, for the truncated models mentioned previously [21,36], models including strong absorption of solute in the bulk of the liquid solution [24] and models of 'reactive' active drops [23], where solute transport is only considered in the interior of the drop and on its interface.
While weakly nonlinear analysis is a natural avenue for studying the spontaneous dynamics of chemically active particles and drops, the existing literature suffers from two significant drawbacks which have so far limited the applicability of this approach. The first is concerned with the manner in which the solvability conditions that constitute the nonlinear amplitude equations have been derived. When faced with the question of solvability of an inhomogeneous linear problem, generally the preferred procedure is to apply the Fredholm Alternative to the forcing terms appearing in that problem [43]; that requires, however, knowledge of the relevant adjoint operator, which in the present context corresponds to the differential operator and auxiliary conditions adjoint to those at linear order of the weakly nonlinear expansion. This knowledge being absent, solvability conditions have to date been derived by solving the inhomogeneous problem in detail using separation of variables. This is a major technical complication, in particular because the linear operator (and its adjoint) is generally isotropic near the threshold for spontaneous motion -with eigensolutions that are axisymmetric about an arbitrary axis -whereas when allowing for general three-dimensional perturbations, or simply unsteadiness, the inhomogeneous problem need not be. This perhaps explains why previous weakly nonlinear analyses have assumed a priori that the motion of the particle is along a line and why genuinely three-dimensional problems, namely where such an assumption is not obvious or does not hold, have yet to be tackled. Furthermore, without employing the Fredholm Alternative the inhomogeneous problem needs to be solved for each physical scenario of interest, despite the relevant linear operator (and its adjoint) typically being common. The only related application of adjoint operators in the literature is a recent study by Farutin et al. [36], who develop a pair of coupled amplitude equations starting from a two-dimensional truncated-domain model. Farutin et al., however, first project the governing equations onto Fourier modes and then separately seek the adjoint operator that is relevant to each mode; as we shall see, that is an unnecessary complication.
The second drawback is that weakly nonlinear theories of active particles and drops near the threshold for spontaneous motion have so far been limited to steady or quasi-steady solutions, the only exceptions being models where the remote advection-diffusion region is absent [31,36]. In conventional weakly nonlinear expansions about the threshold for a monotonic instability, accounting for unsteadiness amounts to identifying the long time scale associated with the closeness to the instability threshold -often just that for which time derivatives first enter at the order of the weakly nonlinear analysis where solvability yields nonlinear amplitude equations. Unsteadiness is then manifested in the appearance of time derivatives in those equations, as is the case in [31,36]. What appears to have been overlooked in the context of active particles and drops is that, when there is a remote advection-diffusion region, the relevant long time scale is actually associated with leading-order unsteadiness in that region. As a consequence, unsteadiness is manifested in the nonlinear amplitude equations as an integral over the history of the particle's motion, physically representing interaction of the particle with its self-generated concentration wake.
In this Series, our goal is to address both of the above drawbacks. We will begin in this part by identifying the adjoint differential operator and auxiliary conditions relevant to the analysis of the canonical model of an isotropic chemically active particle near the threshold for spontaneous motion, as well as general perturbations about that model. To illustrate the application of our adjoint method, we shall, beyond recovering known results in the canonical case, consider the following four perturbation scenarios: (i) weak external force and torque fields, the force scenario being similar to [28] but without making a priori assumptions on symmetry; (ii) small arbitrary perturbations to the surface distributions of solute flux and the slip coefficient; (iii) inclusion of slow first-order kinetics at the surface of the particle, thus going beyond the model of a prescribed solute flux; and (iv) inclusion of weak absorption of solute in the liquid bulk. The first two perturbation scenarios are meant to demonstrate the efficacy of our approach in tackling genuinely three-dimensional problems. The latter two perturbation scenarios are included in order to demonstrate the less obvious applicability of our approach to perturbations involving different physico-chemical mechanisms.
The adjoint method that we develop in this part is relevant to both steady and unsteady weakly nonlinear analyses. This is because the differential operator and auxiliary conditions encountered at linear order of the weakly nonlinear expansion (in the vicinity of the particle) are identical in both cases, owing to the long time scale associated with the closeness to the instability threshold. Nonetheless, the history effects associated with unsteadiness in the remote advection-diffusion region pose a formidable challenge that is rather separate from the solvability problem. We therefore choose to restrict this part to steady solutions and address unsteady problems involving history effects in a subsequent part.
The structure of the paper is as follows. We begin in Sec. II by formulating the steady problem for the canonical model of an isotropic chemically active particle and, closely following the existing literature, constructing the weakly nonlinear expansion for this case in the limit as the Péclet number tends to its critical value. We terminate this section upon arriving at the inhomogeneous problem whose solvability yields the nonlinear amplitude equation.
In this way, we set the stage for developing the adjoint method in Sec. III, ultimately leading to an explicit solvability condition for a generalised inhomogeneous problem which includes the canonical scenario as a special case. In Sec. IV, we employ the generalised solvability condition to consider the perturbation scenarios mentioned above. We give concluding remarks in Sec. V, including a discussion of possible generalisations of our approach and a look ahead to the analysis of history-dependent unsteady dynamics in a subsequent part.

II. WEAKLY NONLINEAR EXPANSION FOR THE CANONICAL MODEL
A. Physical problem We consider a spherical particle of radius a * that is freely suspended in an unbounded fluid of viscosity η * . (An asterisk subscript indicates a dimensional quantity.) The fluid is a liquid solution in which a solute undergoes diffusion and advection, approaching an equilibrium concentrationc * at infinity. The solute also undergoes chemical reactions at the surface of the particle, represented by a prescribed flux j * (positive into the liquid) that is assumed to be uniform. The latter flux is associated with the formation of a solute cloud in the vicinity of the particle, where the concentration perturbation fromc * is of order c * = a * |j * |/D * , D * being the solute diffusivity.
Variations of solute concentration over the surface of the particle drive fluid flow through an effective diffusio-osmotic slip mechanism. In the standard continuum model, the local slip velocity is given by the local product of the concentration surface gradient and a scalar slip coefficient b * , here assumed uniform. This slip mechanism, in conjunction with an order c * concentration perturbation varying on the particle scale, implies the intrinsic velocity scale u * = |j * b * |/D * . In turn, fluid flow can affect the distribution of solute through advection in the liquid bulk. The importance of solute advection relative to diffusion is quantified by the intrinsic Péclet number Pe = a * u * /D * .
Since the particle is spherical and freely suspended, and the prescribed flux and slip coefficient uniform, there exists for all Pe a steady state where the particle and fluid are stationary and the solute cloud surrounding the particle is spherically symmetric. In the case j * b * > 0, this stationary-symmetric state is known to be unstable for Pe > 4 [12].
Numerical simulations of initial-value problems [11,12,15] suggest a supercritical pitchfork bifurcation at the threshold into steady states of spontaneous rectilinear motion of the particle in an arbitrary direction, without rotation. Those numerical simulations also give the speed of that spontaneous motion as a function of Pe, though only in a finite interval above the threshold where the dynamics indeed approach such spontaneous motion. More recently, several authors [14,28] have employed weakly nonlinear analysis to derive the local asymptotic behavior of the spontaneous speed near the bifurcation: assuming a priori that the motion is along a line and that the concentration and flow fields are axially symmetric.
Our goal in this section is to review the derivation of the local bifurcation relation (2.1) by a steady-state weakly nonlinear analysis in the limit Pe → 4. Specifically, we will carry through the analysis up to the final key step, in which (2.1) arises from solvability of an inhomogeneous problem at second order of the weakly nonlinear expansion. This will set the stage for the adjoint method developed in the subsequent section, which furnishes a solvability condition for a generalised inhomogeneous problem relevant to a wide range of physical scenarios, including but not limited to the perturbation scenarios considered in Sec. IV and the extension in the next part to unsteady dynamics.
We stress that the analysis in the present section closely follows [14,28] except that the solvability condition is quoted based on the adjoint method to be developed instead of a detailed solution of the inhomogeneous problem. We also avoid a priori assumptions on symmetry, a generality which will be crucial in Sec. IV when considering non-isotropic perturbations and in the next part when considering unsteady dynamics.
As mentioned above, spontaneous motion occurs only if j * b * > 0. We shall assume, without loss of generality, that j * and b * are both positive.

B. Dimensionless formulation
We adopt a dimensionless formulation where lengths are normalised by a * , concentration by c * , velocities by u * , stresses by η * u * /a * , forces by η * u * a * and torques by η * u * a 2 * . With these conventions, the position vector measured from the center of the particle is denoted by r = rê r , with r = |r| andê r a radial unit vector such that the surface of the sphere is at r = 1 and the fluid domain is r > 1; the solute concentration is denoted by c; the fluid velocity field, in a 'co-moving' frame of reference that moves with the particle's centroid without rotating, is denoted by u; the pressure field associated with that flow is denoted by p; the velocity of the particle centroid is denoted by U; and the particle angular velocity is denoted by Ω. We look for steady-state solutions such that the particle velocities U and Ω are constant vectors and the fields c, u and p are constant in time in the co-moving frame.
The complete formulation is composed of two problems that are mutually coupled. The first can be thought of as governing the concentration field c. It consists of the steady advection-diffusion equation The second problem can be thought of as governing the flow field u and pressure field p, which without loss of generality is taken to decay at infinity. The flow problem consists of the Stokes equations is the stress tensor, I being the identity tensor and † denoting the tensor transpose; the boundary condition wherein ∇ s denotes the surface-gradient operator (see Appendix A); the far-field condition u → −U as r → ∞; (2.8) and the integral conditions where the force F and torque T are defined as wherein dA is an infinitesimal area element.

C. Asymptotic expansions
As previously mentioned, there exists for arbitrary Pe a solution where the particle is stationary, the fluid is at rest, and the concentration field is spherically symmetric. Indeed, the flow problem is then trivially satisfied while the advection-diffusion equation ( In what follows we look for additional steady solutions as Pe → 4, the limiting value being known to be the threshold above which the stationary-symmetric state is linearly unstable. To be consistent with the perturbation scenarios considered later in the paper, it is convenient to write Pe = 4 + ǫχ (2.12) and then consider the limit 0 < ǫ ≪ 1, with χ independent of ǫ. In this section, only the product ǫχ has meaning so χ can, in principle, be chosen as plus or minus unity without loss of generality.
Quantities associated with the flow field are expanded similarly to (2.14). These include the stress tensor σ, force F and torque T, as well as the particle velocities U and Ω.
Expansion (2.13) for the concentration, with r fixed, represents a 'particle region' corresponding to order unity distances from the particle. Allowing r to vary, this expansion fails at large, order 1/ǫ, distances from the particle, corresponding to a 'remote region' where advection and diffusion are comparable. Indeed, if the concentration field decays like 1/r and the velocity field scales like ǫ and approaches a uniform stream, the advection and diffusion terms in (2.2) are of order ǫ/r 2 and 1/r 3 , respectively, for large r. Accordingly, it is necessary to supplement the particle-region expansion (2.13) with a remote-region expansion in which ǫr, rather than r, is held fixed. Asymptotic matching [42] of the particle and remote regions will be employed in order to derive far-field conditions as r → ∞ at subsequent orders of the particle-scale concentration expansion, effectively replacing the decay condition (2.4).
We will find that the assumed decay of the leading term c 0 in the particle-scale expansion is consistent with asymptotic matching, despite that term not holding as a leading-order approximation in the remote region.

D. Remote region
To analyze the remote region, we define the strained position vectorr = ǫr, withr = |r| = ǫr, make the change of variablesc(r; ǫ) = c(r; ǫ) and consider the limit ǫ ց 0 with r > 0 fixed. The 1/r decay of the particle-region concentration, along with the scaling of the remote region, suggests thatc is of order ǫ in that region. We therefore pose the expansioñ c(r; ǫ) ∼ ǫc 1 (r) as ǫ ց 0. (2.15) Unlike the concentration field, the expansion for the flow field is clearly uniformly valid for r > 1. Hence the far-field condition (2.8) and the expansion for U (cf. (2.14)) together imply that u ∼ −ǫU 1 holds in the remote region. It then follows from (2.2) thatc 1 satisfies an advection-diffusion equation with a uniform flow, in which∇ is the gradient operator with respect tor. Furthermore, the decay conditioñ follows from (2.4), while the singular boundary conditioñ follows from asymptotic matching between the one-term remote-region expansion and the particle-region expansion taken to leading order. The solution to (2.16)-(2.18) is [41] Higher-order matching between the one-term remote-region expansion and the particleregion expansion to orders ǫ and ǫ 2 will be seen to provide sufficient information on the far-field behavior of the particle-region fields c 1 and c 2 , respectively.

E. Homogeneous problem
Returning to the particle region, we consider the homogeneous problem governing the order-ǫ concentration and flow field. The concentration problem (2.2)-(2.4) gives the coupled advection-diffusion equation and the boundary condition ∂c 1 ∂r = 0 at r = 1, (2.21) while asymptotic matching with the remote region implies the far-field condition The flow problem (2.5)-(2.9) gives the Stokes equations the boundary condition the far-field condition u 1 → −U 1 as r → ∞; (2.25) and the force and torque balances

F. Spontaneous motion
The order-ǫ problem is homogeneous; it is also linear, except for the nonlinear constant term in the matching condition (2.22), which merely determines a uniform reference value for c 1 . Despite being homogeneous, the order-ǫ problem possesses a family of non-trivial solutions corresponding to steady rectilinear motion of the particle with an arbitrary velocity vector U 1 and without rotation, i.e., Ω 1 = 0. The existence of these solutions is no coincidence but an expected consequence of perturbing about the threshold value of the Péclet number. These homogeneous solutions have been calculated in [14] using separation of variables in spherical coordinates. We here rewrite these solutions as are linear in U 1 , with g(r) being the radial function

G. Inhomogeneous problem
We next consider the inhomogeneous problem governing the order-ǫ 2 concentration and flow fields, whose solvability will be seen to restrict the particle velocity U 1 .

The concentration problem (2.2)-(2.4) gives the inhomogeneous coupled advectiondiffusion equation
and the boundary condition ∂c 2 ∂r = 0 at r = 1. (2.31) The inhomogeneous far-field condition is derived as follows. First, matching the particle-scale expansion to order ǫ 2 with the remote-scale expansion to order ǫ determines the order-r terms. Second, a local analysis of (2.30) as r → ∞, using (2.25), (2.27) and (2.35), is consistent with those order-r terms and further implies the form of the order-unity terms. In particular, the constant term in (2.32), which sets the uniform reference value of c 2 , could be derived by higher-order matching involving a first correction in the remote region. As we shall see, however, only the order-r forcing terms are relevant to the matter of solvability.
The flow problem (2.5)-(2.9) gives the Stokes equations the boundary condition the far-field condition and the integral constraints

H. Solvability condition yields nonlinear amplitude equation
The order-ǫ 2 problem is an inhomogeneous version of the order-ǫ homogeneous problem.
Since the latter problem is singular, we anticipate that the inhomogeneous order-ǫ 2 problem is solvable only under a certain condition on the forcing terms, suggesting a relation between the bifurcation parameter χ and the particle velocity U 1 . In order to derive this solvability condition, it is not necessary to construct detailed solutions to the inhomogeneous problem.
Rather, in the next section we shall derive a solvability condition for a generalised inhomogeneous problem with the help of an adjoint linear operator. The specific solvability condition that is relevant here, which is obtained in Sec. III E, reads The solvability condition (2.37) constitutes the requisite nonlinear amplitude equation, which in the present steady formulation serves as a local bifurcation relation. For arbitrary χ there is always the trivial solution, U 1 = 0, which is consistent with the stationarysymmetric state. For χ > 0, there are also non-trivial solutions corresponding to steady rectilinear motion without rotation, in an arbitrary direction and with speed Rewriting (2.38) in dimensional notation implies the local bifurcation relation (2.1) stated at the beginning of this section. The singular-pitchfork bifurcation implied by (3.25) is schematically depicted in Fig. 1a, recalling that in the present scenario it suffices to consider χ = ±1. As indicated in the figure, the axisymmetric linear-stability analysis in [12] and three-dimensional numerical simulations in [11,12,15] suggest that, for χ > 0, the trivial solution loses stability in favour of the non-trivial spontaneous-motion states.

A. Linear operator and the direct problem
To develop the adjoint method we start by formally defining the 'matrix' differential which acts on concentration-flow 'pairs', namely 'column vectors' ψ ∼ = (c u) T in which c is a concentration field and u an incompressible flow field that has associated with it a pressure field p. In (3.1), S represents the stress-tensor operator and the components of L act on the components of ψ ∼ analogously to standard matrix multiplication. We also define a set of auxiliary conditions, involving particle velocities U and Ω, and we shall say that a pair satisfying these conditions is included in the 'natural domain' of L. These auxiliary conditions consist of the boundary conditions the far-field conditions where E = o(1) and ∇E = o(1) as r → ∞; and the integral constraints where the torque T is defined as in (2.10b) while, henceforth, we adopt the following definition for the force F: . The modified force definition (3.6) agrees with the conventional definition (2.10a) in the case where the stress is divergence-free but is otherwise more specific.
With the above definitions, we define the 'direct' problem where 0 ∼ is the zero pair and ψ

B. Adjoint operator
We now introduce another matrix differential operator, which is essentially the transpose of L, except for the omission of the 'dot' operator from the off-diagonal term (implied by that term now operating on a scalar to give a vector). Similar to L, the operator L * acts on concentration-flow pairs where the flow is incompressible and has associated with it a pressure field. We also define a set of adjoint auxiliary conditions, and we shall say that pairs satisfying these conditions are included in the natural domain of L * . We specify these adjoint auxiliary conditions considering a primed pair (c ′ u ′ ) T , with associated pressure p ′ , stress tensor σ ′ = S[u ′ ] and particle velocities U ′ and Ω ′ . The adjoint auxiliary conditions consist of the boundary conditions where ∇ s · is the surface-divergence operator (see Appendix A); the far-field conditions as r → ∞; and the integral constraints We shall now verify that the operator L * with its natural domain is formally the adjoint of the operator L with its natural domain, with respect to the inner product It is this requirement that has guided us in defining the differential operator L * and its natural domain as we have.
We can write (3.13) explicitly as where we note that the terms involving the off-diagonal coupling terms in (3.1) and (3.8) have cancelled out. The volume integrals in (3.14) can be transformed into surface integrals.
For the first integral, we use the divergence theorem applied to c ′ ∇c − c∇c ′ (i.e., we use Green's second identity). Analogously, for the second integral we use the divergence theorem applied to σ · u ′ − σ ′ · u, the definitions of the stress tensors and the fact that the flows are incompressible (as in the proof of the Lorenz Reciprocal Theorem for Stokes flow [44]).
Overall, we find where ∂D R is the boundary of D R andn is the normal unit vector in the direction outward from D R , with ∂/∂n =n · ∇.
The boundary ∂D R is composed of the surface of the unit sphere and the surface of the sphere r = R. Separating these contributions and applying the natural boundary conditions at r = 1 (cf. (3.3) and (3.9)), we find The first line of (3.16) vanishes by 'integration by parts' over the surface of the unit sphere.
Indeed, we have that 17) and the integral over a closed surface of the surface divergence of a tangential vector field trivially vanishes (see Appendix A). The second line can be seen to vanish using the tripleproduct rule and the fact that the torques T and T ′ vanish for flow fields in the natural domains (cf. (3.5b) and (3.11b)). It remains to show that the last line of (3.16) vanishes.
The limit of the concentration integrals is readily seen to vanish on account of the natural far-field conditions (3.4a) and (3.10a). The stress integrals also vanish, since where we have used the natural far-field conditions and force constrains. We note that this last step rationalises our force definition (3.6).

C. Adjoint problem and adjoint spontaneous motion
Consider now the 'adjoint problem' for a concentration-flow pair ψ ′ ∼ = (c ′ u ′ ) T , say, which is defined by with ψ ′ ∼ restricted to the natural domain of L * .
In this adjoint problem, the coupling between concentration and flow is physically transposed. Thus, while in the direct problem stress is divergence-free and the divergence of concentration flux is proportional to the flow field, in the adjoint problem it is the concentration flux that is divergence-free whereas the divergence of stress is proportional to the concentration. Explicitly, (3.19) gives Laplace's equation and, together with the condition that u ′ is incompressible, the concentration-coupled Stokes The adjoint homogeneous problem, like the homogeneous problem, possesses non-trivial solutions that represent spontaneous rectilinear motion with arbitrary particle velocity U ′ , without rotation, i.e., Ω ′ = 0. In Appendix B, we find these solutions as A property of these solutions that will be useful later is the projection of the traction tangent to the sphere: It will be useful to derive this condition for a class of inhomogeneous problems, formulated below, which is more general than the inhomogeneous problem in Sec. II G. The resulting solvability condition we shall find could then be applied to a wide range of scenarios, including but not limited to the perturbation scenarios we shall consider in Sec. IV.
To be consistent with the notation in Sec. II G, we formulate the generalised inhomogeneous problem for the pair ψ ∼ = (c 2 u 2 ) T , with associated pressure p 2 , stress tensor σ 2 = S[u 2 ] and particle velocities U 2 and Ω 2 . We consider the forcing column vector We consider the boundary conditions where A and B are scalar and vector functions of angular position, respectively, with the latter satisfying the impermeability constraint B ·ê r = 0. The far-field condition on the concentration is prescribed by writing where R is a scalar function of angular position and the remainder E satisfies E = o(r) and ∇E = o(1) as r → ∞. The far-field condition on the flow field is prescribed as Lastly, we pose the integral constraints in which F and T are arbitrary constant vectors. We note that the far-field condition (3.27) is not specific enough to close the generalised inhomogeneous problem. As we shall see, it nevertheless suffices for the purpose of deriving a solvability condition.
It remains to substitute the general form of the adjoint homogeneous solutions (3.22) in order to make the above condition explicit. Note that the term involving the inhomogeneous torque T vanishes trivially, since Ω ′ = 0, and that the term involving B can be simplified using the result (3.23) upon recalling that B has no radial component. Using the fact that U ′ is arbitrary, we arrive at the explicit condition Note that the behavior of C as r → ∞ must be such that the limit on the right-hand side of (3.31) exists.
Condition (3.31) constitutes the main result of our adjoint method, which we shall apply in the next subsection to the canonical model and in the next section to a number of perturbation scenarios. We have here only shown that (3.31) is a necessary condition for existence of a solution to the generalised inhomogeneous problem, by direct calculation. This represents one direction of the Fredholm Alternative theorem for differential operators [43], which if valid in the present setup would imply that (3.25), and so (3.31), also constitute a sufficient condition for solvability of the generalised inhomogeneous problem. Rigorously proving the Fredholm Alternative in our setting, which involves an unbounded domain, as well as integral and incompressibility constraints, is outside the scope of our formal analysis.
From a pragmatic point of view, a necessary condition suffices for the purpose of formally deriving nonlinear amplitude equations.

E. Solvability condition for the canonical model
We now apply the general solvability condition (3.31) to the canonical model of an isotropic chemically active particle considered in the preceding section. Referring to the inhomogeneous problem of Sec. II G, we have in that case that The amplitude equation (2.37) quoted in the preceding section readily follows from (3.31).

IV. PERTURBATION SCENARIOS
To further illustrate the adjoint method developed in the previous section, we go beyond the canonical model of an isotropic chemically active particle to consider the effects of perturbations about that model involving either the particle or its environment. The magnitude of each perturbations will be represented by some small positive parameter, which for the sake of discussion we here denote by δ. For δ arbitrarily small, some perturbations have a leading-order effect sufficiently close to the instability threshold. In accordance with the notation used in Sec. II, we shall still use the small parameter ǫ defined by (2.12) to quantify the closeness of the Péclet number to the unperturbed critical value Pe = 4. Whereas in Sec. II we have carried out a local analysis in the limit Pe → 4, i.e., ǫ ց 0, in the present section we shall analyze distinguished limits where both ǫ and δ are small, the smallness of ǫ relative to δ being such that the perturbation is just strong enough to influence the leading-order amplitude equations. In studying these distinguished limits, we shall slave δ to ǫ such that the parameter χ appearing in definition (2.12) will represent a rescaled and shifted bifurcation parameter that takes on arbitrary real values. This should be contrasted with the canonical isotropic scenario of Sec. II, where χ merely indicated the sign of Pe − 4.
A. Uniform force and torque fields

Uniform force field
The first perturbation we consider is that of a uniform force field [27]. While this scenario has already been analysed in [28], our derivation here -based on the adjoint methodhas two advantages: (i) we do not need to solve the order-ǫ 2 inhomogeneous problem, and (ii) we do not make the a priori assumption that the particle motion is co-linear with the force field. The latter assumption, which is argued in [28] on the basis of the bifurcation curves needing to match with intervals of Pe away from its critical value, follows here from a general three-dimensional analysis. As discussed in [28], this result corresponds to a steadystate alignment of the nominally isotropic spontaneous motion with the direction parallel or anti-parallel to the force field.
We shall represent the external force field by a dimensionless force 6πV acting on the particle, where V is the dimensionless velocity that would be attained by the particle if it was chemically inert. With this convention, the force constraint (2.9a) becomes Let V = Vî, where V > 0 andî is a unit vector, and consider the weak-force limit V ≪ 1 (in the present scenario, V represents the small parameter δ discussed at the beginning of Sec. IV). Inspecting the weakly nonlinear analysis in Sec. II, we see that the force first has a leading-order effect for ǫ = O( √ V ). We thus set V = ǫ 2 , without loss of generality.
The only change to the weakly nonlinear analysis in Sec. II is in the order-ǫ 2 inhomogeneous problem. The homogeneous force constraint (2.36a) becomes the inhomogeneous force constraint anti-parallel to the torque field; note that the spontaneous speed is identical to that in the isotropic canonical case.
The general solvability condition (3.31) applies as in the canonical scenario of Sec. III E, now with F = −6πî. We find the amplitude equation It is readily seen that the vector equation (4.3) reduces to the one-dimensional bifurcation relation obtained in [28]. Indeed, it follows from (4.3) that U 1 = U î, with U satisfying 16U |U | − χU = 2, (4.4) in agreement with [28]. The directionally restricted imperfect-pitchfork bifurcation relation implied by (4.4) is depicted in Fig. 1b. As discussed in [28], the 'parallel' branch (U > 0), given by exists for all χ. As χ → −∞, this branch matches with solutions representing the linear response of the particle to the weak force field at Péclet numbers away from the threshold.
As χ → ∞, it matches with the spontaneous-motion solutions in the canonical scenario when restricted to the direction parallel to the force field. For χ > 8 √ 2, there are additionally two 'anti-parallel' branches (U < 0), given by which are degenerate at χ = 8 √ 2. As χ → ∞, one anti-parallel branch matches with solutions representing the linear response of the particle to the weak force field at Péclet numbers away from the threshold, while the other matches with the spontaneous-motion solutions in the canonical scenario when restricted to the direction anti-parallel to the force field. Stability of the above solution branches will be addressed in the subsequent part of this series.

Uniform torque field
Consider now the effect of an external uniform torque field, represented by a dimensionless torque 8πW, in which W is the dimensionless angular velocity that would be attained by the particle if it was chemically inert. It is clear from the general solvability condition (3.31) that an order-ǫ 2 torque cannot influence the amplitude equation, in contrast to a force field at that order. An order-ǫ torque, however, would generate an order-ǫ rotational flow which could modify the order-ǫ 2 inhomogeneous problem and thence the amplitude equation. To investigate this possibility, we set W = ǫî, withî a unit vector. Accordingly, the torque constraint (2.9b) becomes The weakly nonlinear analysis of Sec. II is first modified at order ǫ of the particle-scale expansion. In light of (4.7), the homogeneous order-ǫ condition (2.26b) is replaced by the inhomogeneous condition The general solution (2.27) to the order-ǫ particle-scale problem therefore needs to be supplemented by a particular solution that accounts for the right-hand side of (4.8). A suitable particular solution is simply provided by the classical Stokes flow corresponding to the fluid velocity and pressure fields, respectively, that would be induced by the external torque if the sphere was chemically inert; this solution is associated with linear and angular particle velocities with c L , u L and p L provided by (2.28), U 1 arbitrary and The order-ǫ 2 particle-scale problem is the same as for the canonical isotropic model, only that u 1 is now provided by (4.11b). In applying the general solvability condition (3.31), the modification to u 1 enters solely through the quantity C, which is still defined by (3.32a). We instead of (3.33a). Otherwise, the calculation is the same as in Sec. III E, leading to the amplitude equation (4.14) It follows from (4.14) that implying the directionally restricted singular-pitchfork bifurcation depicted in Fig. 1c: we have the stationary-symmetric state, U = 0, for all χ, as well as parallel and anti-parallel states, U = ±16/χ, for χ > 0. Thus, the external torque results in alignment, as in the force scenario, but without introducing imperfection, in contrast to the force scenario.

Parallel and perpendicular force and torque fields
Consider now the scenario where the particle is subjected to both an external force field 6πV, of order ǫ 2 , and external torque field 8πW, of order ǫ. Without loss of generality, we set V = ǫ 2î , whereî is a prescribed unit vector, and W = ǫw, where w is a prescribed vector. Combining the analyses of the force and torque scenarios, we find In what follows we assume that neither the force nor torque vanish. In the case where the torque is parallel or anti-parallel to the force, it is readily seen that the torque has no effect, i.e., the steady states are the same as in the force scenario. We henceforth focus on the case where the force and torque fields are perpendicular. Let where w > 0 and we introduce a right-handed orthogonal basis of unit vectors (î,,k). To simplify expressions, we define W, with U = |U|.
We look for solutions U of (4.19) as a function of the real parameter X and the positive parameter W. There are two families of solutions to consider.
1. U = X . In that case, (4.19) degenerates tô 20) to be solved together with U = X > 0. It is readily seen from (4.20) that U ⊥î, namely that the velocity is perpendicular to the force, and thatk · U = −1/W. As depicted schematically in Fig. 2a, the constraint U = X then implies that there are zero solutions for X W < 1, one for X W = 1 and two for X W > 1. Explicitly, those solutions are given by Note that U ∼ ±X as W ր ∞, i.e., these solutions limit to the non-trivial states in the torque scenario as the force magnitude vanishes relative to the torque magnitude.
2. U = X . In that case, the component of (4.19) shows that U ⊥, namely that the velocity is perpendicular to the torque. Inspecting the remaining components of (4.19), we find where U > 0 satisfies the quartic equation which exists for arbitrary X , as well as the pair of solutions 25) which, in the limit, exist for χ > 2 and are degenerate for χ = 2. In the opposite limit, W ր ∞, we find one positive root which vanishes in the limit like U ∼ 1/W.
In Fig. 2b.i, we plot the the positive roots of (4.23) as a function of W, for several values of W, alongside the small-W approximations (4.24) and (4.25).
Combining the small-and large-W limits of the magnitude U and angle φ, we note the following regarding the torque-perpendicular solutions. As W ց 0, there are between one and three solutions that approach the states in the force scenario. As W ր ∞, the only solution approaches the trivial state.
To summarise, for perpendicular force and torque fields, there are between one and five solutions for the particle velocity vector. There are up to two that are perpendicular to the force, as depicted in Fig. 2a. These limit to the solutions in the torque scenario as the magnitude of the force is made negligible in comparison to that of the torque. Furthermore, there are between one and three that are perpendicular to the torque, as depicted in Fig. 2b.
These limit to the solutions in the force scenario as the magnitude of the torque is made negligible in comparison to that of the force; in the opposite limit, there is one such solution that approaches the trivial state.
B. Non-uniform surface properties

Steady formulation and amplitude equation
We next consider weak, generally non-uniform, perturbations to the surface properties of the particle. Specifically, the dimensional prescribed flux and slip coefficient are modified as where α and β are functions of position on the particle boundary, fixed in surface coordinates attached to the particle; the flux and slip-coefficient perturbations have been scaled such that both have a leading-order effect as ǫ ց 0. In this part, we limit ourselves to solutions which appear steady in the co-moving frame introduced in Sec. II B, which does not rotate with the particle. In the present scenario, this implies that the particle does not rotate, or that its rotation is such that it leaves the surface properties fixed in the co-moving frame.
There is also the possibility that the particle rotates sufficiently slowly such that the steady weakly nonlinear theory developed in this part holds in a quasi-static sense; this, however, requires an extremely small rotation rate, Ω = o(ǫ 2 ), to ensure that the associated time variation of U 1 is sufficiently slow such that the steady remote-region analysis in Sec. II D remains valid in a quasi-static sense. In what follows, we shall apply our weakly nonlinear framework assuming a steady state in the co-moving frame, and then derive conditions, associated with the rotation of the particle, for a given candidate solution to be consistent as a steady, or quasi-steady, solution.
Given and R and C as in Sec. III E, we find the amplitude equation where we have used (4.29b), with (2.27a), to show that We see from (4.30) that the effect of the flux perturbation is similar to that of an external force field, the equivalent force field being proportional to a dipole moment of the flux perturbation. The slip-coefficient perturbation modifies the homogeneous linear term, effectively adding a general second-order tensor to the bifurcation parameter χ.
To check whether a solution U 1 to the amplitude equation (4.30) is consistent, we must consider the particle's angular rotation.
Since Ω 1 vanishes trivially in the present scenario, and since Ω = o(ǫ 2 ) is permissible, we need only consider Ω 2 . Fortunately, it is possible to calculate Ω 2 without solving the generalised order-ǫ 2 inhomogeneous problem in detail.
Following [45], we apply the Lorenz Reciprocal Theorem to the Stokes problem included in the order-ǫ 2 inhomogeneous problem, obtaining Ω 2 as a functional of the relative fluid velocity at the surface (cf. (4.28b)): While we do not know c 2 , the associated contribution can be shown to vanish trivially using a Stokes-type integral theorem (see Appendix A; the physical fact that slip uniformly proportional to the surface gradient of a scalar field cannot drive particle rotation was pointed out to me by Ehud Yariv). Accordingly, substitution of (4.29b), with (2.27a), yields A solution U 1 is consistent only if Ω 2 vanishes or is such that it leaves the surface distributions α and β fixed in the co-moving frame. Otherwise, the first term in the amplitude equation (4.30) is wrong, as it originates from matching with an inconsistent steady solution in the remote region.

Axisymmetric perturbation
For the sake of illustration, we henceforth focus on surface perturbations that are symmetric about an axis that points in the direction of the unit vectorp and passes through the particle's centroid, such that α and β are functions of cos θ =p ·ê r , with 0 ≤ θ ≤ π.
We begin by considering the angular rotation in this case. As a consequence of the axial symmetry, we can write 9 16π r=1 dAê r β = β rp , where we define β r = 9 8 π 0 dθ β(θ) sin θ cos θ. We then find from (4.33) that The fact that Ω 2 is perpendicular top rules out the possibility of rotation that leaves the axisymmetric surface distributions fixed in the co-moving frame. A consistent steady (or quasi-steady) state therefore requires Ω 2 = 0, implying that non-longitudinal solutions (U 1 ×p = 0) only represent consistent steady states if β r = 0. We note that β r vanishes trivially for fore-aft symmetric slip-coefficient perturbations, i.e., when β(θ) = β(π − θ).
In Fig. 3, we depict the above solution branches obtained forβ = 0, for bothβ = ±1 and for several values ofᾱ. In the subsequent part, we will investigate the stability of the above solution branches, and explore the dynamical ramifications of having non-longitudinal solutions that cease to be consistent steady states for β r = 0.

C. First-order surface kinetics
In the following two subsections we consider perturbations to the chemical model of the particle and the liquid solution, respectively. Unlike the first two perturbations considered in this section, these perturbations retain the isotropy of the canonical model of Sec. II.
In this subsection, we consider a generalised model for the chemical activity at the surface of the particle. We assume that, in addition to the constant and uniform supply of solute molecules at the particle boundary, represented by the positive flux j * , solute molecules are also absorbed at the surface according to a first-order chemical reaction [13]. The total surface flux can be written as where k * is a rate constant. Withc * the concentration at infinity as in Sec. II B, we assume that j ′ * = j * − k * c * is positive and normalise the concentration deviation fromc * by c ′ * = a * j ′ * /D * , which modifies the characteristic concentration c * defined in Sec. II B. The dimensionless problem is then identical to that formulated in Sec. II B, except that the In the above regime, the weakly nonlinear expansion is modified already at order ǫ.
Recalling that c 0 = 1/r, it follows from (4.47) that the homogeneous boundary condition instead of (2.27).
We are now ready to consider the inhomogeneous problem at order ǫ 2 of the weakly nonlinear expansion. As a final example we consider the effect of solute absorption in the liquid bulk [46].
We assume that the solute is absorbed in proportion to the deviation of the concentration from the equilibrium valuec * . We denote the absorption rate by κ * and adopt the same dimensionless notation as in Sec. II. The dimensionless formulation of the problem is then the same as in Sec. II B except that the advection-diffusion equation (2.2) becomes [46] ∇ 2 c − Peu · ∇c = Da b c, (4.55) where we define a 'bulk' Damkohler number Da b = κ * a 2 * /D * . Naively, the leading-order absorption ∼ −Da b c 0 in the particle region suggests that bulk reactions enter the inhomoge- to be solved together with the decay condition (2.17) and the matching condition (2.18), as in the canonical isotropic scenario. The latter condition is not modified since it is relies on matching with the leading-order particle-scale concentration c 0 , which is not affected by the weak bulk absorption. The solution to the modified remote problem is readily found as which replaces (2.19).
The modified leading-order solution in the remote region implies, through asymptotic matching, modified far-field conditions at orders ǫ and ǫ 2 of the particle-region expansion.
At order ǫ we find that (2.22) is generalised as while at order ǫ 2 we find that (2.32) is generalised as As a consequence of (4.58), the general solution to the order-ǫ particle-region problem is modified from (2.27) to The only difference is in the reference value of the concentration c 1 ; we will see that this has no effect.
Turning to the inhomogeneous problem at order ǫ 2 , the inhomogeneous coupled advectiondiffusion equation (2.30) becomes The quantity C appearing in the solvability condition (3.31) is defined as the right-hand side of (4.61) -this is just (3.32a) plus the isotropic contribution 1/r, which is readily seen to have no effect on the solvability condition. In contrast, the far-field condition (4.59) modifies R from (3.32b) to the expression in the curly brackets of (4.59); we accordingly find that (3.33b) generalises as Substituting the above results into the solvability condition (3.31), we find the amplitude equation and arbitrary direction. In unscaled notation, (4.63) implies that for Pe > 4 + 8 The bifurcation relation (4.63) is depicted in Fig. 4. As already noted in [24], based on a point-particle model where advection is discarded in the vicinity of the particle, we find that bulk absorption 'regularises' the pitchfork bifurcation, viz., the speed initially grows away from the bifurcation like a square-root function rather than linearly as in the canonical isotropic scenario. Since our analysis is limited to weak bulk absorption, the regularisation we find does not fundamentally modify the structure of the weakly nonlinear analysis. In particular, the remote region remains essential and as a consequence the bifurcation still appears 'singular' for √ Da b ≪ Pe − 4 ≪ 1. Stronger bulk absorption would eliminate the remote region.

V. CONCLUDING REMARKS
Our main contribution in this part has been to identify an adjoint differential operator and auxiliary conditions that facilitate the derivation of nonlinear amplitude equations governing the steady velocity of a chemically active particle near the threshold for spontaneous motion. Our adjoint method circumvents the need to directly solve the inhomogeneous problem at quadratic order of a weakly nonlinear expansion, valid as the threshold is approached, making it relatively straightforward to analyze a wide range of perturbation scenarios; the technical simplification stems from the fact that the only problems that need to be explicitly solved, namely the linearised homogeneous problem at the threshold and its adjoint, are axisymmetric (about an arbitrary direction) and common to all scenarios, whereas the inhomogeneous problem at quadratic order is generally not axisymmetric and scenario-specific.
To illustrate our approach, we have derived and then solved steady amplitude equations for a number of perturbation scenarios, demonstrating that sufficiently near the threshold weak perturbations can appreciably modify and enrich the landscape of steady solutions.
While we have introduced the adjoint homogeneous problem of Sec. III C as an auxiliary mathematical tool, we find it intriguing from a physical perspective that the transposed flow-solute coupling associated with that problem (cf. (3.8)-(3.11)) can support spontaneous rectilinear motion of a particle, as represented by the adjoint homogeneous solutions (3.22).
This may suggest an alternative physical mechanism for spontaneous self-propulsion, where a scalar field is associated with a body force on the fluid (hinting to buoyancy or electrostatics) and stress variations give rise to a surface flux of that scalar. As far as we know, such a mechanism for spontaneous motion has not yet been encountered.
We have only considered steady-state solutions in this part. As we shall see in the subsequent part, our adjoint method is equally useful, without modification, when studying the unsteady weakly nonlinear dynamics of a chemically active particle. This is because the weakly nonlinear dynamics evolve on a sufficiently long time scale such that the linear operator at first order of the weakly nonlinear expansion remains quasi-static (in fact, identical to that herein). Inspecting the steady problem formulation and form of the weakly nonlinear expansions in Sec. II, the appropriate long time scale can be deduced as t * = ǫ −2 a 2 * /D * . With t denoting time normalised by t * , generalising the problem formulation of Sec. II B to allow for unsteadiness amounts to the addition of the term ǫ 2 ∂c/∂t to the advection-diffusion equation (2.2); the particle-scale expansion remains quasi-steady at all relevant orders, while the leading-order remote-region equation (2.16) becomes unsteady. As a consequence, the dynamics involve a history effect associated with the particle interacting with its own concentration wake. In the subsequent part, we shall develop unsteady nonlinear amplitude equations which include that history effect. This will allow us to study the stability of the steady solutions found herein, explore transient dynamics, e.g., the alignment of the particle velocity vector with the direction of an external force field, as well as study inherently unsteady scenarios, such as: perturbation scenarios where there are no stable steady states (as we shall see, this can occur in the case of non-uniform slip-coefficient perturbations, leading to stable circular motion); inter-particle and particle-wall interactions; and particles subjected to shear flow or an unsteady force field.
By suitably adapting the adjoint operators found here, it may be possible to develop weakly nonlinear theories for closely related scenarios where physico-chemical activity gives rise to spontaneous dynamics. In particular, it would be of interest to apply a similar approach to experimentally relevant models of so-called 'solubilising' drops, whose activity can be modeled similarly to the canonical isotropic model, but with the flow typically being driven by a Marangoni, rather than diffusio-osmotic, effect [23]. Despite the different physics, we expect the form of the amplitude equations to be similar -both the steady ones developed herein and the unsteady, history dependent, ones to be derived in the subsequent part -such that solutions found for active particles could be adapted to such active drops.
the coordinate-invariant relation ∇ s f = (I −nn) · ∇f between the surface-gradient and gradient operators. In Sec. IV B, we utilise a Stokes-type integral theorem saying that, for a closed surface, we have [47] dAn × ∇ s f = 0. (A2) The surface divergence of a vector field, say P, is defined as with P evaluated on the surface and considered as a function of ν 1 and ν 2 ; comparing with the expression for the divergence operator in the extended curvilinear coordinates, we see that ∇ s · P = [(I −nn) · ∇] · P. Let P = P + P ⊥n , with P = (I −nn) · P and P ⊥ = P ·n.
For a tangential vector field P = P 1û1 + P 2û2 , we find from (A3), using orthogonality and the definitions of the scale factors, the form In Sec. III B, we employ two identities involving the surface-divergence operator. The first, readily follows from (A3). In particular, in (3.17) we have used this identity in the form ∇ s · (f P ) = ∇ s f · P + f ∇ s · P . The second states that, for a closed surface, we have [47] dA ∇ s · P = dA (∇ s ·n)P ⊥ .
In particular, tangential vector fields satisfy the Divergence-like law dA ∇ s · P = 0, which we have used in Sec. III B to carry out integration by parts over the unit sphere.

Appendix B: Adjoint spontaneous motion
Consider the adjoint problem formulated in Sec. III C for the concentration field c ′ and flow field u ′ , with associated pressure p ′ , stress tensor σ ′ and particle velocities U ′ and Ω ′ .
The problem consists of Laplace's equation (3.20), the concentration-coupled Stokes equations (3.21) and the natural auxiliary conditions associated with L * , namely the boundary conditions (3.9), the far-field conditions (3.10) and the integral constraints (3.11). Recall that the force is defined as in (3.6); given that stress is not divergence-free in the adjoint problem, that definition differs from the conventional one (2.10a). We shall constructively show that, like the direct problem defined in Sec. III A (and the order-ǫ homogeneous problem of Sec. II E), this adjoint problem possesses a family of non-trivial solutions describing steady rectilinear motion of the particle with an arbitrary particle velocity, without rotation.
Let θ be the angle between r and U ′ . We introduce spherical coordinates (r, θ, φ), with r defined as in the main text and φ an azimuthal angle, and associated unit vectors (ê r ,ê θ ,ê φ ).