PyStokes: phoresis and Stokesian hydrodynamics in Python

We present a modular Python library for computing many-body hydrodynamic and phoretic interactions between spherical active particles in suspension, when these are given by solutions of the Stokes and Laplace equations. Underpinning the library is a grid-free methodology that combines dimensionality reduction, spectral expansion, and Ritz-Galerkin discretization, thereby reducing the computation to the solution of a linear system. The system can be solved analytically as a series expansion or numerically at a cost quadratic in the number of particles. Suspension-scale quantities like ﬂuid ﬂow, entropy production, and rheological response are obtained at a small additional cost. The library is agnostic to boundary conditions and includes, amongst others, conﬁnement by plane walls or liquid-liquid interfaces. The use of the library is demonstrated with six fully coded examples simulating active phenomena of current experimental interest.


I. INTRODUCTION
PyStokes is a Python library for studying phoretic and hydrodynamic interactions between spherical particles when these interactions can be described by the solutions of, respectively, the Laplace and Stokes equations. The library has been specifically designed for studying these interactions in suspensions of active particles, which are distinguished by their ability to produce flow, and thus motion, in the absence of external forces or torques. Such particles are endowed with a mechanism to produce hydrodynamic flow in a thin interfacial layer, which may be due to the motion of cilia, as in microorganisms [1] or osmotic flows of various kinds in response to spontaneously generated gradients of phoretic fields [2]. The latter, often called autophoresis, is a generalisation of well-known phoretic phenomena including, inter alia, electrophoresis (electric field), diffusiophoresis (chemical field) and thermophoresis (temperature field) that occur in response to externally imposed gradients of phoretic fields [3].
Hydrodynamic and phoretic interactions between "active particles" in a viscous fluid are central to the understanding of their collective dynamics [2,4]. Under experimentally relevant conditions, the motion of the fluid is governed by the Stokes equation and that of the phoretic field, if one is present, by the Laplace equation. The "activity" appears in these equations as boundary conditions on the particle surfaces that prescribe the slip velocity in the Stokes equation and flux of the phoretic field in the Laplace equation. The slip velocity and the phoretic flux are related by a linear constitutive law that can be derived from a detailed analysis of the boundary layer physics [3]. The Stokes and Laplace equations are coupled by this linear constitutive law only at the particle boundaries. The linearity of the governing equations and of the coupling boundary conditions allows for a formally exact solution of the problem of determining the force per unit area on the particle surfaces. This formally exact solution can be approximated to any desired degree of accuracy by a truncated series expansion in a complete basis of functions on the particle boundaries. This, in turn, leads to an efficient and accurate numerical method for computing hydrodynamic and phoretic interactions between active particles [5][6][7].
The principal features that set this method apart are (a) the restriction of independent fluid and phoretic degrees of freedom to the particle boundaries (b) the freedom from grids, both in the bulk of the fluid and on the particle boundaries and (c) the ability to handle, within the same numerical framework, a wide variety of geometries and boundary conditions, including unbounded volumes, volumes bounded by plane walls or interfaces, periodic volumes and, indeed, any geometry-boundary condition combination for which the Green's functions of the governing equations are simply evaluated.
The purpose of this article is to demonstrate the power of the numerical method, as implemented in Python library, through six fully coded examples that simulate experimental phenomena. Our software implementation uses a polylgot programming approach that combines the readability of Python with the speed of Cython and retains the advantages of a high-level, dynamically typed, interpreted language without sacrificing performance.
Our presentation is in the style of literate programming and draws inspiration from similar articles by Weideman and Reddy [8], Higham [9], and Trefethen [10]. The article is best read alongside installing the library and executing the example codes. The library freely is available on GitHub at https://github.com/rajeshrinet/pystokes, where detailed installation instructions can also be found. A subset of the library features is available as a Binder file which requires no installation. All software is released under the MIT license.
The section on features of the library not covered in the examples, features that can be added but have not been, and limitations of the numerical method. We end this Introduction with a brief description of each of the examples. In Section IV we compute the flows produced by the leading terms of the spectral expansion of the active slip for a single spherical particle away from boundaries. In Section V, we examine the effect of boundaries -a plane wall and a plane liquid-liquid interface -and show how the flows in the first example are altered. In Section VI we simulate the Brownian motion of a pair of hydrodynamically interacting active particles in a thermally fluctuating fluid confined by a plane wall. We identify an attractive drag force from the active flow and an unbinding transition with increasing temperature as entropic repulsion overwhelms this active hydrodynamic attraction. In Section VII we simulate the flow-induced phase separation (FIPS) of active particles at a plane wall. This provides a quantitative description of the crystallization of active particles that swim into a wall [11][12][13][14][15][16][17]. In Section IV we introduce a phoretic flux on the particle surface and compute the phoretic field that it produces, both away from and in the proximity of a no-flux wall. In Section IX, we show that a competition between the hydrodynamic and phoretic interactions of autophoretic particles can arrest the phase separation induced by flow alone [7]. In the penultimate section, we show that the cost of computation increases quadratically with the number of particles and decreases linearly with the number of computational threads.

II. MATHEMATICAL UNDERPINNINGS
Our method relies on the reduction of linear elliptic partial differential equations (PDE) to systems of linear algebraic equations. The four key mathematical steps underpinning it are illustrated in this diagram: The first step is the representation of the solution of an elliptic PDE in a three-dimensional volume V as an integral over the boundary of the volume S [5,[18][19][20][21][22][23][24]. For the Laplace equation, this is the classical theorem of Green [25]; for the Stokes equation, it is the generalization obtained by Lorentz [18,19,26]. The integral representation leads to a linear integral equation that provides a functional relation between the field and its flux on S. Thus, if the surface flux in the Laplace equation is specified, the surface concentration is determined by the solution of the Laplace boundary integral equation. Similarly, if the surface velocity in the Stokes equation is specified, the surface traction is determined by the solution of the Stokes boundary integral equation. This transformation of the governing PDEs is the most direct way of relating boundary conditions (surface flux, slip velocities) to boundary values (surface concentration, surface traction). It reduces the dimensionality of the problem from a three-dimensional one in V to a two-dimensional one on S. The second step is the spectral expansion of the field and its flux in terms of global basis functions on S. We use the geometry-adapted tensorial spherical harmonics, which provide an unified way of expanding both scalar and vector quantities on the surface of a sphere. These functions are both complete and orthogonal and provide representations of the three-dimensional rotation group [27]. Thus, symmetries of the active boundary conditions can be represented in a straightforward and transparent manner. The third step is the discretization of the integral equation using the procedure of Ritz and Galerkin [28,29], which reduces it to an infinite-dimensional self-adjoint linear system in the expansion coefficients. This exploits the orthogonality of the basis functions on the sphere. The matrix elements of the linear system can be evaluated analytically in terms of the Green's functions of the respective elliptic equations. The fourth step is the truncation of the infinite-dimensional linear system to a finite-dimensional one that can be solved by standard methods of linear algebra adapted for self-adjoint systems [30]. Analytical solution can be obtained by Jacobi iteration, which is equivalent to Smoluchowski's method of reflections. Numerical solutions can be obtained by the conjugate gradient method, at a cost quadratic in the number of unknowns. From this solution, we can reconstruct the field and the flux on the boundary, use these to determine the fields in the bulk, and from there, compute derived quantities. These steps have been elaborated in several papers [5][6][7]31, 32] and we do not repeat them in detail here. Instead, we show below how the method is applied to problems of experimental interest.

III. LIBRARY STRUCTURE
The overall organization of the library is show in Table I and will be referred to throughout the remainder of the paper. The PyStokes library solves, respectively, the Stokes and Laplace equations, using the reduction method explained in the previous section. The library takes as input a set of expansion coefficients J (l) i for the prescribed active flux j A on the surface of the i-th particle and computes the expansion coefficients of the surface concentration. Table I. This schematic shows the governing equations that determine the hydrodynamic and phoretic interactions between active particles in a three-dimensional domain V . The equations are coupled by the active boundary conditions on the surface Si of the i-th particle. The library takes as input surface fluxes j A , specified in terms of coefficients J The active slip velocity v A is obtained from this using the linear coupling relation v A = µ c ∇ s c [3]. The library outputs the expansion coefficients V (lσ) i of the active slip. The PyStokes library takes as input these expansion coefficients, which may also be specified independently, and any body forces F B i and body torques T B i acting on the particles and returns their rigid body motion in terms of the velocities V i and angular velocities Ω i . In addition to the joint computation of phoretic and hydrodynamic interactions, the PyStokes library can be used to compute the hydrodynamically interacting motion of squirming particles where the slip is specified independently of a phoretic field, or the dynamics of passive sus-pensions where the slip vanishes and forces and torques are prescribed. The PyStokes library can also compute hydrodynamically correlated Brownian motion, and thus, allows the study of the interplay between passive, active, and Brownian contributions to motion. The library optionally computes the corresponding field in V , necessary for insight and visualization. Additionally, PyStokes computes the dissipation of mechanical energy and the rheological response of the suspension.

IV. EXAMPLE 1 -IRREDUCIBLE ACTIVE FLOWS
Our first example shows how to plot the irreducible parts of an active flow around a spherical particle, which we take to be far removed from boundaries. For this example, we will take the reader step by step from the governing PDE to the expressions for the irreducible flows that the PyStokes library evaluates and plots.
1. Elliptic PDE . The flow field v(r) satisfies the Stokes equation in the region V exterior to the sphere (radius b, centered at R, oriented along unit vector p). On the sphere surface S, it satisfies the slip boundary condition where V and Ω are the rotational and translational velocities of the sphere, v A (ρ) is the slip velocity and and ρ is the radius vector from the center to S. The fundamental solution of the Stokes equation is given by the system where G αβ is the Green's function, P α is the pressure vector, K αβγ is the fundamental solution for the stress tensor, and η is the fluid viscosity. The Green's function may need to satisfy additional boundary conditions which we keep unspecified for now.
2. Boundary integral. The fundamental solution, together with the Lorentz reciprocal relation gives the boundary integral representation which expresses the flow field in V in terms of a "single-layer" integral involving the traction and a "double-layer" integral involving the boundary velocity. The latter is specified by the boundary condition; the former must be determined in terms of it.
3. Spectral expansion. The analytical evaluation of the two integrals is possible if the slip and the traction are expanded spectrally in terms of tensorial spherical harmonics, where The orthogonality of the tensorial harmonics implies that The integral of the traction, F (1) , is the net hydrodynamic force and the integral of the cross-product of the traction with the radius vector, the antisymmetric part of F (2) , is the net hydrodynamic torque [6].
4. Ritz-Galerkin discretization. Letting the point r approach S from V and matching the flow in the integral representation with the prescribed boundary condition leads to an integral equation for the traction. Multiplying both sides of the integral equation by the l-th harmonic and integrating yields an infinite-dimensional system of linear equations for the coefficients of the traction, 1 2Ṽ where repeated indices are summed over, for l > 2, and G (l, l ) and K (l, l ) are the matrix elements of the linear system given in terms of integrals of the Green's function and the fundamental stress solution. These can be evaluated analytically [6].

Truncation:
The infinite-dimensional linear system has to be truncated to a finite-dimensional one for tractability. We truncate the system at l = 2 and decompose the coefficients into their irreducible symmetric (s), antisymmetric (a) and trace (t) parts, so that slip and the traction are Here and δ are the Levi-Civita and Kronecker tensors and ∆ αβµν = 1 2 (δ αν δ βµ +δ αµ δ βν − 2 3 δ αβ δ µν ) symmetrises and detraces second-rank tensors. This truncation includes all long-ranged contributions to the active flow and is sufficient to parametrize experimentally measured active flows around microorganisms, active drops, and autophoretic colloids [16,31,32]. The solution of the finite-dimensional linear system yields a linear relation between the irreducible force and velocity coefficients, the "generalized Stokes laws", which are F (lσ) = −γ (lσ,l σ ) ·Ṽ (l σ ) where σ = s, a, t and repeated indices are summed over. In an unbounded domain the friction tensors γ (lσ,l σ ) take on a particularly simple form: they are diagonal in both the l and σ indices, γ (lσ,l σ ) ≡ δ ll δ σσ γ lσ ∆ (l) so that a single scalar γ lσ determines them. For lσ = 1s and lσ = 2a these are the familiar coefficients 6πηb and 8πηb 3 that appear in Stokes laws for the force and torque. Inserting the truncated spectral expansions for the slip and traction in the boundary integral, eliminating the unknown traction coefficients in favour of the known slip coefficients, expanding the Green's function about the center of the sphere and finally using the orthogonality of the tensorial harmonics, we obtain the flow due to each irreducible slip mode as We emphasise that these expressions are valid for any Green's function of the Stokes equation, provided they satisfy the additional boundary conditions that may be imposed. We also note that no discretization of space, either in V or on S is involved, making the result "grid-free". The flow class in PyStokes computes these expressions for supplied values of the force, the torque and slip coefficients. For this example, we choose an unbounded domain with the flow vanishing at infinity, for which the Green's function is the Oseen tensor, The code listed in Fig.(1) computes the irreducible flows 1s, 2s, and 3t for radius b = 1, viscosity η = 1/6, location R = (0, 0, 0) and orientation p = (0, −1, 0). The coefficients are parametrised as F This information is supplied to the flow class which is instantiated for an unbounded fluid. The lσ irreducible component of the flow is computed by the calling the function flow.flowfieldlσ. This is passed to a generic plotting function to compute the streamlines in a plane of symmetry. These are shown in the code output where the polar symmetries of the 1s and 3t modes and the nematic symmetry of the 2s mode can be seen clearly. For vanishing radius, these are the Stokeslet, potential dipole and stresslet singularities [33].

V. EXAMPLE 2 -EFFECT OF PLANE BOUNDARIES
Our second example illustrates how irreducible flows are modified by the proximity to plane boundaries. This is of relevance to experiments, where confinement by boundaries is commonplace [16,34]. This also illustrates the flexibility of our method, as the only quantity that needs to be changed is the Green's function. The Green's function for a no-slip wall is the Lorentz-Blake tensor Here is the image of the j-th colloid at a distance h from plane boundary and M = I − 2ẑẑ is the reflection operator. The Green's function for a no-shear plane air-water interface is The plane boundary is placed at z = 0 and the flows are plotted in the half-space z > 0. The irreducible flows for each boundary condition are obtained by evaluated Eq. (8) with the corresponding Green's function. The irreducible flows for the modes in Example 1 are shown in Fig.(2), with no-slip wall in the top panels and no-shear interface in the bottom panels. We use the same initialization as in Example 1, but instantiate wall-bounded and interfaciallybounded classes wflow and iflow respectively. Though arbitrary viscosity ratios are allowed, PyStokes assumes an air-water interface as default, setting the viscosity ratio between the two fluids to zero, as in this example. The lσ irreducible component of the flow is computed by the calling the function wflow.flowfieldlσ and passed to a generic plotting function for streamline computation. Note that the streamlines near an interface do not close, unlike those near a plane wall. The perpendicular component of the flow near a wall is about an order of magnitude larger than that near an interface. These features have been recently used to understand the control by boundaries of flow-induced phase separation of active particles [16].

VI. EXAMPLE 3 -ACTIVE BROWNIAN HYDRODYNAMICS
Our third example shows how to simulate the dynamics of active particles including active, passive and thermal forces. This example assumes that orientational degrees of freedom are not dynamical, as is often the case with strongly bottom-heavy active particles. The PyStokes library computes the rigid body motion of the active particles consistent with the overdamped Langevin equation where F H i = F (1s) i , F B i andF i are the hydrodynamic, body and thermal forces respectively and i is the particle index. We now consider a minimal model of the slip, retaining only the lowest modes of vectorial symmetry: We parametrize the coefficients of the slip uniaxially, in terms of the orientation p i of the i-th colloid and its selfpropulsion speed v s , as The hydrodynamic forces are then given by the generalized Stokes laws, where sub-dominant rotational contributions have been neglected. The explicit forms of the tensors γ T T ij and γ (T,lσ) ij follows from the solution of the many-particle version of the linear system given in Eq. (6). The off-diagonal terms, with i = j, represent hydrodynamic interactions and are obtained as an infinite series in the Green's function and its derivatives [6]. With this, force balance becomes This shows that in the absence of slip modes with l > 1, external forces, thermal fluctuations and hydrodynamic interactions, the translational velocity of the i-th particle is given by −V which is the basis for our Brownian dynamics algorithm. Here ξ T j is a vector of zero-mean unit-variance independent Gaussian random variables, the mobility matrix µ T T ij is the inverse of the friction tensor γ T T ij and the propulsion tensor π . We now look at the dynamics of a pair of active particles, i = 1, 2, near a no-slip wall. We assume a truncated Lennard-Jones two-body interaction between the particles and a different truncated Lennard-Jones one-body interaction with the wall, designed to prevent particle-particle and particle-wall overlaps. The orientation is taken to point into the wall and a sufficiently large torque is added to prevent re-orientation. Vertical motion ceases when repulsion from the wall balances active propulsion into it. Then, the net external force on the particle points normally away from the wall. This implies that the leading contribution to active flow is due to the F (1s) mode, as shown in Fig.(2). This flow drags neighboring particles into each other and leads to the formation of a bound state [32]. This is the basis for numerous aggregation phenomena of active particles near walls and interfaces [11][12][13][14][15][16][17]. In terms of equations of motion, the z−component of force balance now reads where the first term is the active propulsive force in the z-direction, the second term is the z−component of the net force, and the third term is the noise. The solution of this equation implicitly gives the mean heighth above the wall at which the particles come to rest. The x−component of force balance gives where the friction tensors are now evaluated at the mean height and the instantaneous separation between the particles. This effectively decouples the vertical and horizontal components of motion. The first two terms are the self-and mutual-Stokes drags, the third term is the drag from the active hydrodynamic flow, the fourth term is the body force and the last term is the noise. The above overdamped Langevin equation can be written in standard Ito form as where the matrix of variances containing σ ij is the Cholesky factor of the matrix of mobilities µ xx ij . These are the equations we simulate in Example 3. The truncated Lennard-Jones interactions are U = ( rmin r ) 12 − 2 ( rmin r ) 6 + , for separation r < r min , and zero otherwise [36]. In the absence of thermal fluctuations, the pair form a bound state due to the attractive active drag but at sufficiently large temperatures, they "unbind" due to entropic forces. Remarkably, the distribution of their in-plane separation can be expressed in Gibbsian form, with a non-equilibrium potential, as shown in Fig. (3). We now explain why this is so.
Denoting the x-component of the active hydrodynamic drag as F A x and using the leading form for γ xz 21 we have where γ and γ ⊥ are the self-friction coefficients in the directions parallel ( ) and perpendicular (⊥) to the wall [37]. This shows that the active drag force can be written as the gradient of a potential whose strength depends on the propulsion speed. Potentials of identical functional form, but with different prefactors, have been obtained before for electrophoresis [38] and thermophoresis [39] without associating them to an active drag force, as we have done here. In spite of the non-equilibrium origin of the potential, it leads to a Gibbs distribution for the particle positions, P ∼ exp[−(Φ + U )/k B T ], as the Ito equation satisfies potential conditions whenever γ xx 12 , and by Onsager symmetry γ xx 21 , is a gradient. This understanding of the active Brownian hydrodynamics of a pair of bottom-heavy active particles near a plane wall rationalizes the ubiquitously observed crystallization of active particles near plane boundaries [11][12][13][14][15][16][17], which is our next example.

VII. EXAMPLE 4 -FLOW-INDUCED PHASE SEPARATION AT A WALL
Our fourth example demonstrates the flow-induced phase separation (FIPS) of active particles at a plane no-slip wall. We use the same model of active particles as described in the Example 3 but consider a large number of them.  . Flow-induced phase separation of active particles at a plane no-slip wall. Starting from a non-crystalline distribution, active particles crystallize into a single cluster due to long-ranged active hydrodynamic interactions between them [32].
As before, the slip is truncated to contain only 1s and 3t modes of vectorial symmetry, which are parametrized in term of the orientation p i of the colloids, see Eqs. (13) and (14). All colloids are oriented along the normal to the wall, such that p i = −ẑ. In experiment, the active force ∼ 6πηbv s , is of the order 10 −13 N in [12,40], while 10 −11 N in [15]. Thus, it is orders of magnitude larger than the Brownian force k B T /b ∼ 10 −15 N, which we ignore for this example.
The code and snapshots from the simulations are shown in Fig.(4). The initial random distribution of positions is provided by pystokes.utils.initialConditionRandom(Np). The rigid body motion class is instantiated as pystokes.wallBounded.Rbm, while the Lennard-Jones inter-particle and particle-wall forces are obtained from classes pystokes.forceFields.Forces. A short function computes the velocity that is passed to a standard Python integrator which advances the system in time and saves positional data. This is used to plot the snapshots at specified time instants.

VIII. EXAMPLE 5 -IRREDUCIBLE AUTOPHORETIC FIELDS
Our fifth example introduces phoretic fields and, as in Example 1, shows how to plot the irreducible parts of an active phoretic field around a spherical particle. We take the reader step by step from the governing PDE to the expression for the irreducible phoretic fields that the library evaluates and plots. The notation for the particle size, # ex5 . py : chemical field of an autphoretic colloid -unbounded domain ( first row ) and near a wall ( second ) import pystokes , numpy as np , matplotlib . pyplot as plt where j A (ρ) is the active flux and D is the diffusivity. The fundamental solution of the Laplace equation is by where H is the Green's function which may need to satisfy additional boundary conditions. The normal derivative of the Green's function is L = Dρ α ∇ α H.
2. Boundary integral. The fundamental solution, together with Green's identities, gives the boundary integral representation which expresses the phoretic field in V in terms of a "single-layer" integral involving active flux and a "doublelayer" integral involving the boundary phoretic field. The former is specified by the boundary condition; the latter must be determined in terms of it.
3. Spectral expansion. The analytical evaluation of the two integrals is possible if the phoretic field and the flux are expanded spectrally in terms of tensorial spherical harmonics, where C (m) and J (m) are l-th rank tensorial coefficients, symmetric and irreducible in all their indices. The orthogonality of the tensorial harmonics implies that 4. Ritz-Galerkin discretization. Letting the point r approach S from V and requiring the phoretic field to attain its value on the boundary leads to an integral equation for it. Multiplying both sides of the integral equation by the l-th harmonic and integrating yields an infinite-dimensional system of linear equations for the coefficients of the phoretic field, where the matrix elements H (l, l ) and L (l, l ) are give in terms of the Green's function and its normal derivative. These can be evaluated analytically [7].

Truncation:
The infinite-dimensional linear system has to be truncated to a finite-dimensional linear system for tractability. We truncate the system at l = 2, so that the phoretic field and its flux are The solution of the finite-dimensional linear system yields a linear relation between the coefficients of the field and its flux, the "elastance relations", C (m) = −ε (m,m ) · J (m ) . In an unbounded domain, the elastance tensors ε (m,m ) take on a particularly simple form: they are diagonal in the m indices, ε (m,m ) ≡ δ mm ∆ (m) /4πDw m so that the single scalar 4πDw m determines them. For m = 0 this is the familiar coefficient 1/4πb for the inverse capacitance of a spherical conductor.
Inserting the truncated spectral expansions for the phoretic field and the active flux, eliminating the unknown phoretic coefficients for the known flux coefficients, expanding the Green's function about the center of the sphere and finally using orthogonality of the tensorial harmonics, we obtain the phoretic field due to each irreducible flux mode as We emphasise that these expressions are valid for any Green's function of the Laplace equation, provided they satisfy the additional boundary conditions that may be imposed. For this example we consider the Green's function in an unbounded domain, where the flux vanishes at infinity, and in a domain confined by an infinite planar wall at which the flux vanishes, Here, as earlier, r * = M · r is the mirror image of the point, where M = I − 2ẑẑ is the reflection operator.
The code listed in Fig.(5) computes the irreducible concentration profile for m = 0,1,2 modes of the active flux for radius b = 1, diffusion constant D = 1, location R = (0, 0, 5) and orientation p = (0, 0, 1). The coefficients are parametrised as J (0) = 1, J (1) α = p α and J (2) αβ = p α p α − δ αβ /3. These are supplied to the ufield class which is instantiated for an unbounded domain. The m-th irreducible component of the field is computed by the calling the function ufield.phoreticFieldm. This is passed to a generic plotting function to compute the pseudocolor plot in a plane of symmetry. These are shown in the code output where the spherical , polar, and nematic symmetries of the m = 0, 1 and 2 modes can be seen. The second row shows the same fields near a plane wall, obtained by changing the flow instantiation to wfield class. The fields have reduced symmetry in z−direction due to the introduction of the plane wall.

IX. EXAMPLE 6 -AUTOPHORETIC ARREST OF FLOW-INDUCED PHASE SEPARATION
Our sixth example describes the modification of active slip by phoretic interactions between active particles and the resulting arrest of flow-induced phase separation seen in Example 4. The slip, from being an apriori specified quantity for each particle, must now be computed from the interacting phoretic fields of all particles. In the previous example, we obtained the phoretic field for a specified surface flux. Here we show how to coupled it to the Stokes equation and obtain both the hydrodynamic and phoretic interactions of active colloids. The coupling of the Stokes and Laplace equation and Laplace equations is through the following expression for the active slip [3], In the previous section, we showed that C . This solution can be used to obtain the coefficients of the slip as V (l) , where χ (l,m) is a coupling tensor of rank (l + m) that depends on the phoretic mobility µ c [7]. Thus the problem is fully specified once the coefficients of the flux and phoretic mobility on the surface of all the particles is given.
For this example, we consider an active surface flux j A with spherical and polar modes and a phoretic mobility that is constant, Phoretic interactions are computed with he no-flux condition at the plane wall using the Green's function of Eq. (32). Hydrodynamic interactions are computed with the no-slip on the plane wall using the Lorentz-Blake tensor of Eq.(10).

X. CODE PERFORMANCE
Our final example is on performance and benchmarks. We show that our codes scales linearly with the number of CPU cores and quadratically with N , the number of particles simulated. In the current implementation, the velocities of about 10 5 particles can be computed in a few seconds for a mode of the active slip, as shown in Fig.(7).
With current many-core architectures a dynamic simulation of about N ∼ 10 5 is within reach. For larger number of particles, accelerated summation methods are desirable, which can reduce the cost to O(N log N ) [41,42] or even O(N ) [43][44][45]. These methods can be implemented in the present numerical architecture as an improvement over the direct kernel sum, while maintaining the overall library structure.

XI. WHAT ELSE ?
In the six examples above we demonstrated the use of PyStokes library for simulating hydrodynamic and phoretic phenomena. These do not exhaust the capabilities of the library and much else can be done with them. We conclude by listing implemented, implementable and unimplementable features. Implemented but not shown: the library supports periodic geometries [46] and parallel plane walls [16,47]. Polymers [48], membranes and other hierarchical assemblies of active particles can be simulated. Can be implemented: other boundary conditions, for instance flows interior to a spherical domain such as a liquid drop, near-field lubrication interactions, and numerical solutions of the linear system can be implemented in the current design, with treecode or fast-multipole accelerations. Cannot be implemented: irregular geometries for which analytical forms of the Green's functions of the Laplace and Stokes equations cannot be evaluated analytically and/or irregularly shaped particles on whose boundaries globally defined spectral basis functions are not available. A shorter version of this paper has been published in JOSS [49].