Numerical simulations of self-diffusiophoretic colloids at fluid interfaces

The dynamics of active colloids is very sensitive to the presence of boundaries and interfaces which therefore can be used to control their motion. Here we analyze the dynamics of active colloids adsorbed at a fluid-fluid interface. By using a mesoscopic numerical approach which relies on an approximated numerical solution of the Navier-Stokes equation, we show that when adsorbed at a fluid interface, an active colloid experiences a net torque even in the absence of a viscosity contrast between the two adjacent fluids. In particular, we study the dependence of this torque on the contact angle of the colloid with the fluid-fluid interface and on its surface properties. We rationalize our results via an approximate approach which accounts for the appearance of a local friction coefficient. By providing insight into the dynamics of active colloids adsorbed at fluid interfaces, our results are relevant for two-dimensional self assembly and emulsion stabilization by means of active colloids.

The dynamics of active colloids is very sensitive to the presence of boundaries and interfaces which therefore can be used to control their motion.Here we analyze the dynamics of active colloids adsorbed at a fluid-fluid interface.By using a mesoscopic numerical approach which relies on an approximated numerical solution of the Navier-Stokes equation, we show that when adsorbed at a fluid interface, an active colloid experiences a net torque even in the absence of a viscosity contrast between the two adjacent fluids.In particular, we study the dependence of this torque on the contact angle of the colloid with the fluid-fluid interface and on its surface properties.We rationalize our results via an approximate approach which accounts for the appearance of a local friction coefficient.By providing insight into the dynamics of active colloids adsorbed at fluid interfaces, our results are relevant for two-dimensional self assembly and emulsion stabilization by means of active colloids.
Figure 1.Illustration of the system of a self-diffusiophoretic colloid (sphere) at a fluid-fluid interface (red, solid line) between two fluids with dynamic viscosities η1 and η2 and reactant number densities ρ1 and ρ2.The blue line represents the inert region of the colloid surface, while the orange, solid line represents the catalytic region with opening angle θo.The colloid is partially wetted by both fluids and its position with respect to the interface is captured by the contact angle θc.The orientation of the colloid with respect to the interface is characterized by the angle θ between the axis o of the colloid and a line parallel to the interface (black, solid line).Here the orientation of the particle corresponds to θo > 0.
In this contribution we analyze the dynamics of self-diffusiophoretic colloids (i.e., particles inducing local gradients in the chemical potentials of certain suspended species) which are adsorbed at a fluid-fluid interface.The fluid flow induced by the local stresses is caused by the chemical reaction on the surface of the particle and it will be affected by the presence of two, phase separated fluid phases.The description of such a system via the standard coarse-grained approach, within which the relative velocity between the particle and the fluid is accounted for by the so-called phoretic slip velocity on the surface of the particle [9,[25][26][27][28], might be insufficient.Indeed, the phoretic slip velocity has been invoked in those cases in which the imbalance in the local chemical potential is confined to a thin shell around the particle where the Stokes equation is solved analytically [25].This approach becomes more complicated if the active colloid is adsorbed at a fluid interface in the presence of a three-phase contact line.In this context, we present a novel approach, based on numerical simulations, in which the motion of the self-diffusiophoretic colloid is obtained by using the lattice Boltzmann method in order to construct approximate solutions of the Navier-Stokes equation directly.In our scheme this hydrodynamics solver is combined with an advection and diffusion equation for the reactants.Such an approach allows us to discuss the reliability of the slip-velocity approach by comparison with previous approximate analytical results [19].In particular, in the present study we focus on the case in which the two fluids have the same viscosity, for which the approximate analytical model [19] predicts the absence of any torque on the particle.Interestingly, our results show that, even in this case, self-phoretic colloids trapped at fluid interfaces reorient their symmetry axis.This reorientation occurs whenever the axis of symmetry of the particle is not perpendicular or parallel to the interface.Indeed, for these cases the presence of the interface affects the velocity profile and leads to net torques on the particle.
The presentation of our study is organized as follows.In Sec.II we describe our numerical method based on lattice Boltzmann simulations.In Sec.III we report our results for the dynamics of self-diffusiophoretic colloids trapped at fluid interfaces, and in Sec.IV we summarize our main findings.

II. NUMERICAL METHODS
Our system is composed of two phase separated fluids (e.g., oil and water) acting as solvents, the reactants, and the products of the chemical reaction (such as the decomposition of hydrogen peroxide into water and oxygen), and the colloid (see Fig. 1).In order to determine the dynamics of the system, we put forward diverse numerical approaches for describing each of these components.

A. The lattice Boltzmann method
In order to solve the dynamics of the fluids we use the lattice Boltzmann method (LBM) as it is implemented in the LB3D package [29,30].Within the LBM the fluid phases are described by their discretized single particle distribution functions f σ i (r, t), which give the probability of finding a fluid particle of component σ at position r with velocity c i .Our system consists of two species, such as σ 1 = oil and σ 2 = water, forming two fully segregated phases.Here, we use a so-called D3Q19 lattice with 19 discrete velocities c i , i = 1, ..., 19, in three dimensions.We measure times in units of the integration time step ∆t and lengths in units of the lattice spacing ∆x.These microscopic auxiliary quantities have no physical meaning and their values are chosen to be smaller than any other physically relevant length or time scale.Eventually, it turns out to be convenient to fix the magnitudes of ∆t and ∆x to unity and to measure times and lengths in units of ∆t and ∆x, respectively.Accordingly, the actual dimensional values of ∆t and ∆x in actual units follow from the smallest length and time scale with physical meaning1 .The particle distribution functions f σ i (r, t) evolve in time due to advection from the neighboring lattice sites and due to collisions among particles at the same lattice site.In the following we use the so-called Bhatnagar-Gross-Krook (BGK) collision operator [31].After some algebra involving the discretization of space, time, and velocities [32] the time evolution of the distribution functions f σ i (r, t) follows as where the rhs of Eq. ( 2) is the BGK collision operator, τ σ is the relaxation time of the fluid component σ, and f σ,eq i (r, t) is the local equilibrium distribution function which, in the small Mach number limit, is given by [33] where ζ i are the lattice weights [29,30] and m σ is the mass of species σ.The relaxation time is related to the kinematic viscosity of the fluid as ν σ = (c σ s ) 2 τ σ − ∆t 2 , where c σ s = kB T m σ is the speed of sound in the phase dominated by species σ.Our numerical approach, in the present form, requires all species to have the same mass m σ = m, and hence the same speed of sound c σ s = c s .In the following we choose to fix the lattice time step ∆t and the lattice spacing ∆x to unity.In these units, it is common to choose c s = 1 √ 3 ∆x ∆t as the lattice speed of sound in terms of the time step ∆t and the lattice spacing ∆x [32].Once the distribution functions f σ,eq i (r, t) are known, it is possible to compute the local mass density of the fluid: where m σ is the mass of a single particle of species σ.The barycentric velocity of the fluid mixture is Finally, the velocities of the individual fluid components are given by In order to account for multiple solvent phases we follow the method introduced by Shan and Chen [33].Within this method the interaction force density acting among distinct species has the form where g σσ ′ denotes the interaction strength between the components σ and σ ′ ; ψ σ is a dimensionless pseudo-potential, which is a functional of the mass density.Here, the functional form of ψ σ is chosen as2 where ρ σ 0 is a reference density which is related to the bulk properties of the phase dominated by species σ.The force density in Eq. ( 6) is applied to the fluid by adding a shift to ū: Accordingly, in the expression for the equilibrium distribution functions f σ,eq i (r, t) , ū is replaced by ū′ (see Eq. ( 8)).Values of the interaction strengths g σσ ′ between distinct species, i.e., σ = σ ′ , which exceed a threshold, eventually lead to their separation, whereas values of g σσ ′ for σ = σ ′ exceeding the threshold give rise to the separation of the liquid and the vapor phases of a certain species [33].In the following we take g σσ ′ = 0.1 and g σσ = 0 which leads to an interface with a thickness of ca. 5 lattice units and to a surface tension of the order of 0.1 in lattice units.
In order to study the behavior of a colloid suspended at a fluid-fluid interface, a scheme is needed for treating objects, which are large compared with the particles forming the fluids.The separation of length scales between the mesoscopic colloidal size and the molecular size of the fluid particles allows one to keep the coarse-grained description for the fluid (via LBM) while simultaneously treating the colloid as a spherical object characterized completely by its size, position, orientation, and its linear and angular velocities.The interaction between the colloid and the fluid gives rise to forces and torques acting on both the colloid and the fluid.The technical details of the implementation of the coupling between the colloid and the fluid within LBM are discussed in Refs.[34][35][36][37][38].In the following we shall outline only the basic features of this method.
The colloid occupies those lattice cells which are inside the spherical colloid (see Fig. 2).As the colloid moves, the configuration of lattice cells occupied by the colloid is updated.Solid impenetrability is accounted for by bouncing back those contributions to the fluid flow which attempt to invade the solid boundaries [37,39].This is implemented by updating the distribution function after the streaming step according to for all i, if the lattice site r − c i ∆t is occupied by fluid particles, whereas if r − c i ∆t is occupied by the colloid, the fluid particles are bounced back, i.e., their velocity is flipped: where i ′ is defined as the index corresponding to c i ′ = −c i .This procedure leads to a no-slip boundary condition at the surface of the colloid and to a momentum transfer between the fluid species σ and the colloid, which induces a local force density and a torque density where r(t) is the vector pointing from the center of the colloid to the site where bounce-back occurs.In order to be consistent the above mentioned bounce-back rule (see Eq. ( 10)) has to be modified by accounting for the motion of the colloid [37,39].In order to do so, a correction is added to the fluid distribution functions [37,39]: where u surf is the local velocity at the surface of the colloid, r is the position of the fluid lattice site.(Note that Eq. ( 13) holds explicitly in the limit ∆t → 0 due to Eq. ( 10).)Consequently, the force density acting on the colloid is modified, too: Finally, when the colloid moves it occupies pristine cells at its front and it releases cells at its back.In the case of newly occupied cells, the fluid located therein is deleted and its momentum is transferred to the colloid by adding to the rhs of Eq. ( 14).In the case of cells being released by the colloid, fresh fluid is created at the corresponding sites with velocity u surf (r, t) and density ρσ , where ρσ is obtained by averaging the fluid composition in the direct neighborhood [37]: where the sum ′ i is restricted to those values of i for which r + c i ∆t is a fluid site; N is the number of these sites.In addition, in order to account for the unknown exact density profile in the close vicinity of a colloid, we apply a density correction algorithm as explained in Refs.[37,38].In order to conserve momentum, a contribution is added to the force density (Eq.( 14)): In order to avoid artifacts during the computation of the Shan-Chen forces acting on the colloid and to tune the wetting properties of the colloid, the outermost layer of lattice sites occupied by the colloid is filled with a virtual fluid which is only used during the computation of the Shan-Chen forces acting in the direct vicinity of the colloid.Its density is obtained similarly to Eq. ( 16), but can be tuned by adding an offset ∆ρ to Eq. ( 16), which controls the wettability of the colloid.Finally, as usual the Shan-Chen forces are used to compute the force acting on the colloid: Accordingly, by tuning the magnitude of ∆ρ, it is possible to control the contact angle of the colloid [37,38].

B. Dynamics of the reactant and the reaction product
We assume that the two fluids forming the adjacent phases act as reservoirs of the reactant and of the reaction product such that the mass density ρ r of the reactant is regarded as to be homogeneous in the bulk of both fluid phases.Close to the interface the mass density of the two fluid phases varies.In the following we assume that the ratios of the number densities of the reactant, f r , and those of the fluid phases, are kept constant throughout both fluid phases.Accordingly, the mass density of reactants at a given position is defined as: where the sum runs over all fluid species σ and m r is the mass of a reactant molecule.Equation ( 21) implies that ρ r (r) = κm r f r (r) where κ is the number of fluid species.Concerning the interaction between the reactant molecules and the colloid, we assume that the surface of the colloid is covered by a catalyst with axial symmetry (see Fig. 2(a)).
The strength of the chemical reaction is controlled by the surface activity ξ, which determines the rate at which reactant molecules are converted into solute molecules.A simple choice of the surface activity ξ is given by The prefactor ξ 0 is the base activity, θ ′ is the azimuthal angle (see Fig. 2), θ o is the opening angle which defines the size of the catalytic cap (see Fig. 2(a)), and ∆ is an interpolation length.This corresponds to full activity for angles well within θ o , zero activity outside, and a linear interpolation in between.Without this linear interpolation (i.e., for ∆ = 0), small rotations of the colloid may not change the activity, due to the roughness of the colloid surface 3 .Optimally, the value of ∆ should correspond to ca. one lattice cell, which is accomplished if ∆ ≈ ∆x/R, where R is the radius of the colloid.This is the quantity which we have used in our simulations.
The local flow J of the reaction product, produced by the chemical reaction at the surface of the colloid, into a fluid cell at r per time is given by where α is the reaction rate, n is the local normal at the colloid surface pointing towards the fluid phase, ρ rp is the reaction product mass density, and Ξ(r) denotes the activity of the colloid cells neighboring the fluid cell at r: where the sum is performed only over neighboring cells r p inside the colloid, and Θ is the Heaviside function.This recipe prevents fluid cells in contact with more than one catalyst cell to be exposed to a larger flux, induced by the discretization of the shape of the colloid.Equations ( 23) and ( 24) enforce a constant flux per area even in the discrete representation of the colloid.
In the following we specialize on the case of a single reaction product (our approach can be easily extended to an arbitrary number of them).Moreover we focus on the small Péclet limit, in which the advection contribution to the time evolution of the density of the reaction product is negligibly small as compared with the diffusion contribution.Accordingly, the dynamics of the mass density of the reaction product is decoupled from that of the two fluid phases.
The former appears as an additional scalar field defined on the LBM lattice, such that a cell at lattice position r contains a mass density ρ rp (r): where, as for the reactant (see Eq.( 21)), we have introduced the number density ratios C rp σ (r) and we have accounted for their spatial variation.In order to ensure that ρ rp (r) reaches a steady state within the simulation box with periodic boundary conditions along the three spatial directions, we introduce a homogeneous sink term with constant decay rate4 χ.Hence the mass density of the reaction product is governed by the partial differential equation with the boundary condition Here |r − r c | is the distance of a fluid cell from the center of mass of the colloid, located at r c (see Fig. 2).In Eq. ( 26) we have introduced the interaction potential between the colloid and the reaction product where F 0 and l are the strength and the range of the potential, respectively.In the following we assume l = 4 lattice units and βF 0 l = 5 × 10 −5 (∆x) 3 where β is the inverse thermal energy.This potential has no angular dependence, i.e., it is the same for the catalytic and the inert side of the colloid surface.Equation ( 26) is solved via a finite-difference scheme on the same grid as the one used by the LBM.Finally, the potential U , together with the non-equilibrium mass density profile ρ rp of the reaction product, induces a pressure gradient on the fluid phases as [25] ∇p

III. RESULTS
After validating our numerical scheme against analytical predictions for the velocity of a self-diffusiophoretic colloid in a homogeneous fluid [40] (see Appendix A), we study the dynamics of a self-diffusiophoretic colloid adsorbed at a fluid interface.All simulations are initialized by equilibrating the interface after the colloid has adsorbed; then the chemical reaction at the surface of the colloid is turned on.In the following, we analyze the dynamics of the colloid in the case that the viscosities of the two fluids are equal, η 1 = η 2 = η, and for various values of the opening angle θ o of the catalytic cap and of the contact angle θ c (see Fig. 1), which characterizes the adsorption of the colloid at the interface.Concerning the mass density of the reactant we consider two situations.In the first one, both fluids have the same mass fraction of reactant molecules: i.e., their mass fractions are assumed to be equal In the second case one has so that one fluid does not contain any reactant molecules.Finally, we study two different scenarios: one, in which the colloid can move freely along the interface, and another one, in which the lateral colloid position is fixed by an external force.

A. Equal reactant mass fractions
First, we consider the case in which the reactants are suspended with equal mass fraction in both fluid phases (Eq.( 30)).In this case certain symmetries can be identified, depending on the relative orientation θ of the axis of the colloid with respect to the plane of the interface: i) The fore-aft symmetry for θ = ±π/2, i.e., if the axis of the colloid is parallel to the normal of the interface.(ii) For c 2 = c 1 = 0 an additional mirror symmetry appears about θ = 0, i.e., if the axis of the colloid lies within the plane of the interface.
As a first case, we study the dynamics of a colloid which is partially covered by catalyst, i.e., here θ o = π/4, and has a contact angle θ c = π/2.Figure 3 shows the orientation θ(t) of the colloid as a function of time for a free colloid (Fig. 3(a)) as well as for a colloid the center of mass of which is kept at a fixed position by an external force (Fig. 3(b)).Interestingly, the dynamics in the two setups are quite similar, in that the catalytic reaction induces a net torque on the colloid which leads θ to approach a steady state with θ(t = ∞) ≡ θ ∞ = 0. We note that, while for the moving colloid θ ∞ = 0 is the only steady state, for the fixed colloid θ ∞ = π/2 is also a steady state.This difference is due to the fact that if the colloid is not fixed, its center of mass will move also along the direction normal to the interface.Hence, the distinction between the two cases emphasizes that θ ∞ = π/2 is not stable with respect to fluctuations of the position of the center of mass of the colloid.
Next, we study the dependence of the dynamics of the colloid on the areal size of the catalytic cap, characterized by θ o , and on the contact angle θ c . Figure 4 shows the stable steady state orientations θ ∞ as a function of the opening angle θ o for various contact angles of a free colloid (Fig. 4(a)) and of a colloid the center of mass of which is kept fixed (Fig. 4(b)).Interestingly, both panels of Fig. 4 show that for opening angles θ o = π/2 the steady state orientation of the colloid is θ ∞ = π/2, i.e., the colloid attains a steady translation along the interface, because the driving force points into the axial direction of the colloid and thus provides a lateral component.In particular, θ ∞ grows upon increasing θ o for both θ c = 0.55π (red upward triangles) and θ c = 0.6π (green downward triangles).Finally, for θ o = π/2, the stable steady state is θ ∞ = ±π/2 for all contact angles θ c which leads to a vanishing velocity along the interface.We note that the case of contact angle θ c = π/2 (blue circles) is peculiar because in this case the steady state orientation is θ ∞ = 0 which implies maximum velocity along the interface.This result holds for all opening angles θ o < π/2.For θ o = π/2 with θ c = π/2 we face numerical difficulties due to the symmetry of the problem.In fact, for θ c = π/2 the center of mass of the colloid lies at the interface and if the colloid is half-covered (θ o = π/2) the only non-motile state requires θ = 0.However, when θ → 0, the typical size of the active site being wetted by one of the two fluid phases becomes comparable to the lattice constant and, for θ c = π/2 and θ o = π/2, our discrete numerical approach would require larger ratios of the colloid size and the lattice constant which, however, we could not explore.
In the case in which the colloid can move freely, at steady state its lateral velocity along the interface depends on both the contact angle and the opening angle.Figure 5 shows that for θ c = π/2 the maximum speed is obtained for θ o = π/2 and it equals the one (v 0 ) obtained in a homogeneous fluid.This tells that, under the condition of equal viscosity among the two fluid phases, due to the symmetry of the problem, the interface does not affect the fluid flow and hence the motion of the colloid.In contrast, Fig. 5 shows that, for θ c = π/2, v has a non-monotonous dependence on the opening angle.Indeed, v vanishes for θ o = 0 (i.e., for a passive colloid) and for θ o = π/2, and it attains a maximum for θ o ≃ π/4.This non-monotonous dependence of v on θ o can be understood by relating the data in Fig. 5 to those shown in Fig. 4(a).Indeed, for θ o → 0 at steady state the axis of the colloid lies in the plane of the interface (i.e., θ(t = ∞) → 0), whereas upon increasing θ o the steady state orientation, θ ∞ , also grows and so does the steady state velocity.For even larger values of θ o the value of θ ∞ increases until, for The rotation of the symmetry axis of the colloid we have just described is rather counter-intuitive, in particular for θ c = π/2 and η 1 = η 2 .Indeed, the equality of the viscosity of the two fluids and the lack of accumulation of the reaction product at the interface enforce the symmetry of the mass density profile of the reaction product about the symmetry axis of the colloid and therefore of the local pressure gradients acting on the fluid.Hence, at first glance, in this scenario one might expect no net torque.Actually, for this case a simplified analytical model [19], which disregards the fluid flow in the vicinity of the three-phase contact line, predicts that there is no reorientation of the colloid at all.In the following, we show that the observed torque is due to the fact that the mobility is not homogeneous along the surface of the colloid.Indeed, the boundary conditions imposed on the fluid velocity by the presence of the interface lead to an effective local mobility the spatial variation of which generates the torque.In order to highlight this relationship, we determine the net velocity, which a passive colloid attains if it is pushed by an  33)) as a function of the position in the x-y plane through the colloid center (black dot).The interface lies in the x-z plane located at y = 0 (black horizontal line).The blue color tells that a force at this position leads to a negative angular velocity ωz < 0, whereas a force in cells with red color leads to a positive angular velocity ωz > 0. The inverse mobility coefficient is only shown for positions which are at most l = 4 lattice units away from the colloid (R < r < R + l), where l is the range of the interaction potential (see Eq. ( 28)).The blue circle represents the colloid.
external local force density f r0 (r), localized at r 0 near the colloid surface and acting towards the center of the colloid: where δ (r − r 0 ) is the Dirac delta function, F is a force, and r c is the position of the center of mass of the colloid.In these corresponding simulations, the colloid attains a steady state in which it rotates with constant angular velocity ω.We are interested only in its component ω z (t = ∞) = e z • ω(t = ∞) along the z-direction5 .We introduce the dimensionless inverse mobility coefficient as the ratio of the angular velocity and the magnitude of the applied force 6 , where R is the radius of the colloid, γ 0 = 6πηR is the friction coefficient of the colloid; γ(r 0 ) depends on r 0 via ω(t = ∞).
In Fig. 6, the dimensionless inverse mobility coefficient is shown as a function of the x and y positions in a plane perpendicular to the interface (black horizontal line in Fig. 6) and passing through the colloid center of mass.This highlights the point symmetry of the dimensionless inverse mobility and that forces closer to the interface give rise to larger torques and hence higher values of the steady angular velocity ω z .In the small Reynolds number regime, which is valid for the motion of the present diffusiophoretic colloid, the Navier-Stokes equation reduces to the (linear) Stokes equation.Hence, the local mobility allows us to determine the torque on the colloid in z-direction for any given force density f (r), as long as the force acts on the colloid along the radial direction.For the case of a self-diffusiophoretic colloid, the force density distribution, due to phoresis, is given by f ph (r) = ∂U(r) ∂r ρ rp (r) m rp (see Eq. ( 29)) where here r indicates the distance from the center of mass of the colloid.Accordingly, the local mobility allows us to determine the torque τ z = τ • e z exerted on the colloid by the density profile of the reaction product: where B is a dimensional fitting parameter, and L is the size of the simulation box.Actually, Eq. ( 34) approximates the total torque in that we are considering solely the contributions to the torque stemming from the plane perpendicular to the interface and passing through the center of mass of the colloid (see Fig. 6).Contributions stemming from the rest of the surface of the colloid are accounted for by the fitting parameter B. We emphasize that the same value of the parameter B, i.e., B = 5.2 ∆x, fits well the data for colloids with various opening angles θ o as shown in Figs. 7, 10, 11, and 12.The validity of this approach will be discussed a posteriori, i.e., by comparing it with the values extracted from the corresponding lattice Boltzmann simulations.
Figure 7 provides a comparison of the torque calculated via Eq.( 34) with the torque obtained directly from the lattice Boltzmann simulations for a colloid with opening angle θ o = π/4 and for three contact angles.(For a discussion of these kind of results for various opening angles see Appendix B.) Interestingly, Fig. 7(a) shows that for a contact angle θ c = π/2 the prediction of Eq. ( 34) agrees very well with the torque obtained from the simulations for both a free and a fixed colloid.As expected, Fig. 7(a) shows that the torque on the colloid is zero for the steady states characterized by θ = 0, ± π 2 .For contact angles θ c = π/2 (see Figs. 7(b) and 7(c)), the torques on the free and on the fixed colloid are no longer equal.This is expected, because for these contact angles the lateral motion along the interface leads to an additional contribution to the torque [41].This explains the mismatch between the prediction of Eq. ( 34) and the torque of the free colloid obtained from the simulations shown in Figs.7(b) and 7(c).However, this additional torque is absent for the fixed colloid.This explains the good agreement between the prediction of Eq. ( 34) and the torque on the fixed colloid calculated from the simulations shown in Figs.7(b) and 7(c).

Inhomogeneous reactant densities
Up to now, we have studied the case in which both fluid phases contain reactant molecules.Here, we consider the system in which one of the two fluids does not contain the reactant (Eq.( 31)), yet the product of the chemical reaction can diffuse, with equal diffusivity, in both fluid phases.Clearly, in this case the point symmetry of τ z (θ) at θ = 0 is broken even for the contact angle θ c = π/2, whereas its mirror-symmetry at θ c = ±π/2 remains.Since we have not observed major discrepancies between the case of a moving colloid and that of a colloid fixed in space, in the following we focus on the case of a fixed colloid, because it is computationally easier.If the catalytic cap is fully immersed in that fluid phase which contains the reactant, the dynamics of the colloid is identical to that observed if the reactant is dissolved in both fluid phases.From geometric considerations (see Fig. 1), it follows that the catalytic cap is exposed to the fluid without reactant if θ ≥ θ 1 with and it is fully submerged in the fluid without reactants for θ ≥ θ 2 with Accordingly, as long as θ < θ 1 the catalytic cap does not get into contact with the fluid without reactant and hence the colloid will behave exactly as in the case in which both fluids contain the reactant.If θ > θ 1 , a large gradient of the solute concentration occurs close to the interface due to the lack of reaction in the other fluid phase.This drives the rotation of the colloid further towards the fluid without reactant, until the catalytic cap becomes fully submerged in this fluid, i.e., θ ≥ θ 2 .At this point the system becomes "passive" because there is no longer production of solute and thus the colloid does not move anymore.
Figure 8 shows the steady states θ ∞ of the colloid as a function of the opening angle.We remark that Fig. 8 does not show all stable steady states with θ > θ 2 , i.e., when the system becomes passive due to the lack of chemical reactions.Interestingly, the steady states shown in Figure 8 are in good agreement with our simple geometrical estimate (Eq.( 36)) except for a constant offset.We speculate that this offset is a numerical artifact because a small concentration of reactant molecules can occur in the interface region.(We recall that in the lattice Boltzmann scheme the fluid-fluid interface has a finite thickness, which in the present simulations is ca. 5 lattice constants.)In particular, within the Shan-Chen model the density of the fluid, which contains the reactant, is not zero everywhere because the two fluids do not completely demix.Therefore, when the colloid protrudes far into the fluid without reactant, there is a torque rotating it back towards the fluid with reactant, in line with the results from Section III A. Finally, we note that for a Janus colloid, i.e., θ o = π/2, and with θ c > π/2 the steady state θ ∞ = π/2 is independent of the contact angle, because the Janus colloid is always in contact with the fluid with reactant.

IV. CONCLUSIONS
We have presented a novel numerical approach which is capable of capturing the dynamics of self-phoretic colloids even in the presence of a fluid interface which affects the boundary condition of the fluid velocity close to the surface of the colloid.In particular, we have characterized the dynamics of a self-diffusiophoretic colloid adsorbed at a fluid-fluid interface for the case in which the colloid is let free to move as well as when the colloid center of mass is kept at a constant position by an external force.We have found that a rotation of the axis of symmetry of the colloid arises even for fluids with equal viscosity.In particular, we have found that the steady state orientation depends on both the opening angle θ o of the catalytic cap and the contact angle θ c .
In order to understand the origin of such a reorientation we have calculated the local mobility by applying locally external, constant forces on a passive colloid.We have extracted the local value of the mobility by taking the ratio between the applied force and the steady-state angular velocity.This local mobility has been used to determine the effective torque acting on the self-diffusiophoretic colloid.Interestingly, we have found that, for θ c = π/2, the mobility matrix predicts quite well the effective torque on a free colloid as well as on a colloid the center of mass of which is kept spatially fixed.If θ c = π/2 we have found good agreement for the case of the fixed colloid, whereas an additional torque arises in the case of a mobile colloid.Finally, we have discussed the case in which only one fluid contains the reactant.In this case we have found that the stable steady-state orientation is the one with the cap immersed in the fluid without reactant so that the colloid becomes inactive.

Figure 2 .
Figure 2. The solid sphere (a) represents the idealized colloid, while the cubic cells represent the lattice Boltzmann cells (b).Fluid cells are white, colloid cells are blue, and catalytically active cells are orange.The center of the cells is indicated by a black dot.

Figure 3 .
Figure 3. Angle θ of the orientation of the colloid with respect to the interface (see Fig.1) as a function of time, normalized by the time t0 it takes for a half-covered colloid to move a distance which equals its own radius, t0 = R/v0, for the case of a free colloid (panel (a)) and for the case of a fixed colloid (panel (b)).The velocity v0 is the velocity the very same colloid attains in a homogeneous fluid characterized by the same viscosity (see Appendix A).The lines correspond to various initial orientations θ(t = 0) of the colloid [θ(t = 0) = 0, π/8, π/4, 3π/8, π/2] for θo = π/4.In the case of the free colloid, the trajectories have been smoothed by integrating over a shifting time-window 10 4 ∆t wide.The lines are symmetric with respect to θ = 0.

Figure 6 .
Figure 6.Inverse mobility coefficient (red/blue color) γ (Eq.(33)) as a function of the position in the x-y plane through the colloid center (black dot).The interface lies in the x-z plane located at y = 0 (black horizontal line).The blue color tells that a force at this position leads to a negative angular velocity ωz < 0, whereas a force in cells with red color leads to a positive angular velocity ωz > 0. The inverse mobility coefficient is only shown for positions which are at most l = 4 lattice units away from the colloid (R < r < R + l), where l is the range of the interaction potential (see Eq. (28)).The blue circle represents the colloid.

Figure 7 .
Figure 7. Component τz in the z-direction of the torque on the colloid as a function of the rotation θ.Blue (green) solid lines correspond to the LBM measured torque for a fixed (free) colloid; the error bars indicate the standard error.(Note that for certain values of θ the blue, green, and red lines totally overlap and thus are not visible.)Red solid lines correspond to the torque calculated via Eq.(34) for θo = π/4 with (a) θc = 0.5π, (b) θc ≈ 0.55π, and (c) θc ≈ 0.6π.The torque components are normalized by τ0 ≡ τ (θ = −π/4, θc = π/2) as obtained from Eq. (34) with θc = π/2 (τ0 = 0.04 in LBM units).The error bars are due to the discretization of the colloid surface (for both fixed and mobile colloids) and due to the motion normal to the interface (mobile colloid, green lines).

Figure 9 .
Figure 9.The velocity v of a self-diffusiophoretic colloid as a function of the opening angle θo, normalized by v0 = v(θo = π/2) (see Fig.5).The symbols are simulation data whereas the solid line provides the analytical prediction[40].