Instabilities driven by diffusio-phoretic flow on catalytic surfaces

We theoretically and numerically investigate the instabilities driven by diffusiophoretic flow, caused by a solutal concentration gradient along a reacting surface. The important control parameter is the Peclet number Pe, which quantifies the ratio of the solutal advection rate to the diffusion rate. First, we study the diffusiophoretic flow on a catalytic plane in two dimensions. From a linear stability analysis, we obtain that for Pe larger than 8pi, mass transport by convection overtakes that by diffusion, and a symmetry-breaking mode arises, which is consistent with numerical results. For even larger Pe, non-linear terms become important. For Pe larger than 16pi, multiple concentration plumes are emitted from the catalytic plane, which eventually merges into a single larger one. When Pe is even larger, there are continuous emissions and merging events of the concentration plumes. This newly-found flow state reflects the non-linear saturation of the system. The critical Peclet number for the transition to this state depends on Schmidt number Sc. In the second part of the paper, we conduct three-dimensional simulations for spherical catalytic particles, and beyond a critical Peclet number again find continuous plume emission and plume merging, now leading to chaotic motion of the phoretic particle. Our results thus help to understand the experimentally observed chaotic motion of catalytic particles in the high Pe regime.


Introduction
Self-propulsion at the micrometer scale frequently occurs in nature (Lauga & Thomas 2009;Lauga 2016;Bray 2000;Jeanneret et al. 2016). For example, micro-organisms selfpropel to search for nutrients, different temperatures, or sunlight. Inspired by such motile . Schematic illustration of the catalytic particles (red) with chemical reaction and diffusio-phoresis near the interface. The product originating from the catalytic reaction at the particle surface is shown in cyan. (b) shows a zoom-in of (a). If there is a concentration gradient at the interface, a slip velocity is induced (diffusio-phoresis). Beyond a critical reaction rate (expressed as a critical Péclet number), such gradient emerges through a linear instability.
biological organisms, extensive studies on artificial micro-swimmers have been done over the last one and a half decades, especially on self-propelled phoretic particles (Jiang et al. 2010;Moran & Posner 2017;Golestanian et al. 2007;Qi et al. 2020;Maass et al. 2016;Bär et al. 2020;Jin et al. 2017). Also dissolving or chemically reacting droplets can show such phenomena Li et al. 2019;Vajdi Hokmabad et al. 2019;Lohse & Zhang 2020). A typical feature of the self-propelled particles is that, instead of swimming with appendages, they can propel themselves by converting free energy from the environment into kinetic energy (Ramaswamy 2010;Ebbens & Howse 2010). The driving mechanism behind the propulsion of phoretic particles is diffusio-phoresis (Anderson 1989). Note that in some literature the terminology "diffusio-osmotic effect" is used to indicate the same mechanism. The basic feature is that whenever there exists a tangential concentration gradient on the surface of the particle, there is an induced flow within the interaction layer adjacent to the surface, as shown in figure 1. Since the layer is much thinner than the size of the object, the flow is conveniently described with a slip velocity at the surface (Golestanian et al. 2007). This effect can also be generalized to other coupled fields such as the temperature or the electric fields. The resulting flows are respectively referred to as thermo-phoretic or electro-phoretic (Piazza 2008;Squires & Bazant 2006;Long et al. 1999;Moran & Posner 2011).
The classical mathematical framework for the study of self-propelled particles has often neglected the effect of solute advection (Golestanian et al. 2007). Michelin et al. (2013), however, have revealed that the Péclet number P e is an important parameter controlling the motion of self-propelled particles. P e is the ratio of the solute advection to the diffusion rates. Through a linear stability analysis and corresponding simulations, Michelin et al. (2013) found that when P e is larger than the critical value P e cr = 4, a spherical active particle by dissolution and chemical reaction exhibits a motion in a preferred direction which breaks the rotational symmetry of the system. Later, Michelin & Lauga (2014) performed a comprehensive theoretical study on how the moving speed of the active particle depends on P e, and generalized the theory to any coverage of the reacting surface.
For large enough P e, some fascinating features can emerge. Hu et al. (2019) have numerically observed that for large enough Péclet numbers, such an active isotropic particle acquires chaotic trajectories. Analogously, in the problem of active droplets, Ruckenstein (1981) also found similar helical or chaotic motions, as caused by the interfacial Marangoni flow (Suga et al. 2018;Maass et al. 2016;Herminghaus et al. 2014). Though the phenomena of active droplets and active particles look similar, Krüger et al. (2016) explained that the helical trajectory of the active droplet is attributed to the coupling between the internal flow and the direction of the nematic field, whereas such internal flow is obviously absent in particles. In a further study, Morozov & Michelin (2019a) have considered both Marangoni and diffusio-phoretic effects into their numerical simulation, and also demonstrate that a chaotic oscillation of the droplet can occur. Very recently, Michelin et al. (2020) investigated a simplified system, namely a uniform phoretic channel. They reported spontaneous symmetry-breaking of the solute distribution which provides a route to understand the propulsion of isotropic active particle. However, it remains necessary to understand how the Schmidt number (defined as the ratio of kinematic viscosity to solute diffusivity) influences the diffusio-phoretic instability. Furthermore, the high Péclet number regime is still not fully explored, and we will see that there an interesting chaotic flow arises.
A flow very related to the diffusio-phoretic flow is Bérnard-Marangoni convection, where a spontaneous flow instability occurs too. In that system, the flow is driven by a surface tension difference caused by a variation of the temperature at the fluid surface (Pearson 1958;Davis 1987;Bergeon et al. 1998;Boeck & Vitanov 2002). These two systems share some similarities in their symmetry-breaking mechanism and with the chaotic flow motion at high enough Péclet or Marangoni number. However, the two systems are different problems with diffusio-phoretic flow being driven by the phoretic velocity at the surface and Bénard-Marangoni being driven by the difference in surface tension.
Motivated by the above mentioned recent findings, in this paper we focus on the instability due to chemical reactions and the resulting diffusio-phoretic flow near a catalytic interface, especially in the large P e regime. To start with some reduced complexity, we first consider a simplified model, namely diffusio-phoretic flow over a catalytic plane, in order to study the dynamics near the catalytic surface (see figure 2). This simplified model can reproduce the important features of the diffusio-phoretic flow, and it is also convenient to avoid the added complexities arising from the curvature of the surface. In the second part of the paper we go beyond the simplified model and numerically examine the plume emission and merging phenomena for chaotically-moving phoretic particles.
The paper is organized as follows: After a description of the problem setup and the control parameters in Section 2, the linear stability analysis for the catalytic plane system is performed in Section 3. Then the numerical method and numerical setup are provided in Section 4.1. The numerical results for the catalytic plane are presented in Section 4.2. Then we extend our research to phoretic particle in Section 5. Finally, conclusions and outlook to further work are given in Section 6. The details of the linear stability analysis is given in Appendix A.

Problem setup and control parameters
We start with the two-dimensional system sketched in figure 2. The domain has periodic boundary conditions on both sides and a catalytic plane at the bottom. The width and height of the domain are denoted by L and H. The physical variables to describe the system are the concentration of the productĉ(x, y, t) and the velocity of the fluidû(x, y, t). Note that all dimensional physical variables are marked with hat (e.g.ĉ, u), while the dimensionless ones without (e.g. c, u). At the catalytic surface, chemical reactions take place which convert the reactant into the product. By assuming a constant reaction rate, the concentration boundary condition of the product at the bottom plane is given by where D is the diffusivity of the product in the fluid and α measures the strength of the reaction activity at the catalytic surface, i.e. the generation of solute by the reaction. The tangential concentration gradient induces a slip velocity at the surface of the plane. This is the so-called diffusio-phoretic flow, which is parallel to the surface and its magnitude is proportional to the tangential concentration gradient. The relationship between the induced slip velocity and the tangential concentration gradient is given bŷ where M is the phoretic mobility. The sign of M can either be positive or negative, depending on the type of the solute-surface interaction (Anderson 1989). Michelin et al. (2013) prove that the diffusio-phoretic system is unstable only if M α is positive. In this work, we study the case M > 0 and α > 0. The time evolution of the concentration field c(x, y, t) and the velocity field u(x, y, t) = (u(x, y, t), v(x, y, t)) are governed by the Navier-Stokes equations and the convectiondiffusion equation. The characteristic scales for non-dimensionlization are M α/D for velocities, L for lengths and αL/D for concentrations. The dimensionless form of the governing equations can then be written as: where P e is the Péclet number, characterizing the ratio of the solutal advection rate to the diffusion rate and Sc the Schmidt number, characterizing the ratio between the momentum and mass diffusivities: Note that for the case of fixed P e and Sc → ∞, the pressure term in equation (2.4a) should be rescaled as p = p Sc/P e , since the pressure gradient always exists, even in Stokes flow, where it then balances the viscous forces.
In dimensionless form, the concentration boundary condition at the catalytic plane becomes ∂c ∂y y=0 = −1. (2.6) The tangential velocity is proportional to the tangential concentration gradient, and its dimensionless form is while the normal component of the velocity vanishes at the plane surface v| y=0 = 0. (2.8) Both the velocity and the concentration boundary conditions at the top wall are zero, In Section 3, we have conducted the linear stability analysis with a semi-infinite domain. We note that for aspect ratios H/L > 0.8 (discussed in Section 4), the growth rate of the instability becomes insensitive to the aspect ratio. Therefore, we can compare the results on a linear stability analysis for a semi-infinite domain to the numerical results for a finite domain with H/L = 1.

Linear stability analysis for catalytic plane
In this section, the linear stability analysis is performed to investigate the stability of the system. In the linear stability analysis we add small amplitude perturbations to the basic state: u =ū(y, t) + ũ(x, y, t), p =p(y, t) + p(x, y, t), c =c(y, t) + c(x, y, t), (3.1) whereū,p andc are the basic state of the velocity, pressure and concentration fields, and ũ, p and c are small perturbations with the coefficient 1. A trivial solution to the basic configuration is a static state with zero velocity and pressure. Substitute zero velocity into (2.3), the concentration field is (see Appendix A, see also Wu et al. (2006), p. 144) For t → ∞ and any finite y, we obtain the following concentration gradient   Figure 3. (a) Stability diagram for the catalytic plane in the P e vs Sc parameter space for wavenumber n = 1. An eigenvalue s > 0 indicates instability. The color represents the actual value of s, i.e., the strength of the exponential growth. When P e > 8π, s is positive and the system is unstable, independently of Sc. (b) s as a function of P e at various wavenumber for Sc = 1 by linear stability analysis. The wavenumber of the curve increases from left to right. For wavenumber n, when P e < 8πn, s is negative and the system is stable. When P e > 8πn, s becomes positive and the system becomes unstable towards this mode n. The function of the maximum growth rate curve (dashed line in (b)) is equation (3.10).
The boundary conditions become: We now assume as ansatz a separation of variables and periodic behavior in the lateral direction, such that the perturbation can be written as: ( u(x, y, t), v(x, y, t), p(x, y, t), c(x, y, t)) = (ǔ(y),v(y),p(y),č(y))e ikx+st , (3.7) where k = 2πn and n ∈ N is the wavenumber. This is the standard normal mode analysis (see, for example, Drazin & Reid (2004)). Note that the sign of s determines the flow stability of the system: s > 0 means exponential growth, or instability (the larger, the more unstable), whereas s < 0 indicates stability. Combining the above equations (3.4)-(3.7) and the boundary conditions, we obtain the relation which allows us to calculate how the stability depends on P e and Sc (detailed derivations are in Appendix A): (3.8) Assuming s = 0 in equation (3.8), we get the critical P e for transition from stability to instability for different wavenumber n, P e cr = 4k = 8πn, (3.9) Note that P e cr is independent of Sc.
If we combine equation (3.8) and its derivative with respect to n, we obtain the function of the maximum growth rate at different wavenumber for Sc = 1 (for a detailed derivation Simulation with non-linear terms Simulation with linear terms only Figure 4. Time evolution of the kinetic energy E k for the case P e = 125 and Sc = 1 with random perturbation from simulations with only linear terms (blue solid curve) and with both linear and non-linear terms (red solid curve). The kinetic energy E k is in log scale. For the case with both linear and non-linear terms, the growth of E k levels off near time instant b, compared to that with only linear terms, because of non-linear saturation. The process is divided into two subprocesses: plume generation and plume growth and merging. During the first subprocess, E k grows exponentially E k ∼ e 2st . The points at the curve represents four states in the process, of which the concentration fields are shown in panels a to d, respectively. (a) Plume generating: Triggered by a perturbation, the kinetic energy increases exponentially. (b-d) Plume growing and merging. As E k reaches around 0.02, the kinetic energy reaches a plateau; at the same time the plumes emerge in the concentration field. The plumes grow and merge with each other. In the end, only one major plume remains in the field (d).
see the Appendix A): (3.10) and the corresponding wavenumber is (3.11) where the symbol represents the calculation of the nearest integer. The exponential growth rate s as obtained from equation (3.8) as function of P e and Sc for the case of n = 1 is shown in figure 3(a). Moreover, s as a function of P e for different wavenumbers and Sc = 1 is plotted by the solid curves in figure 3(b). The dashed line shows the maximum growth rate curve, which is equation (3.10). The way to calculate figure 3 from equation (3.8) is explained in Appendix A.
The linearized diffusion-convection equation (3.4) helps us to understand the physical mechanism of the diffusio-phoretic instability. If there is local concentration variation at the surface, the diffusion term 1 P e ∇ 2 c smoothes out the local concentration difference, which makes the system stable. In contrast, dominance of the advection term −u·∇c will increase the concentration difference, such that the system becomes unstable. Thus it can be seen that the competing mass transport by diffusion and advection determines the instability, which is quantified by P e. If P e is above a critical value, the advection term results in positive feedback, which amplifies the disturbance and leads to the instability.

Simulation of catalytic plane
We now numerically study the diffusio-omostic instability. The objective of the numerical simulaiton is to understand the effect of the non-linear terms and random initial perturbation which are ignored in the linear stability analysis.

Numerical setups
The fluid motion and concentration field are solved using direct numerical simulation (DNS) of the Navier-Stokes equations and diffusion-convection equation in Cartesian coordinates. Equations (2.3)-(2.4) are spatially discretized using the central secondorder finite difference scheme. Along both horizontal and vertical directions, homogenous staggered grids are used. The equations are integrated by a fractional-step method with the nonlinear terms computed explicitly by a low-storage third-order Runge-Kutta scheme and the viscous terms computed implicitly by a Crank-Nicolson scheme (Verzicco & Orlandi 1996;van der Poel et al. 2015). The simulations are then conducted with the concentration and the velocity boundary conditions written in equations (2.6)-(2.9). We first examine how the growth rate responds to the domain size. Figure 5(a) shows that when the aspect ratio H/L 0.8, the growth rate of the instability becomes insensitive to the aspect ratio. Besides, the mesh refinement test is given in figure 5(b), from which we see the convergence of the growth rate when the number of grid points in one direction has reached roughly 300. Therefore, we chose H/L = 1 and the mesh 401 × 401 for all our phoretic channel simulations.
The initial condition is the fluid at rest and a constant concentration gradient along the y-direction (see equation (3.3)). Then a small sinusoidal perturbation is added to the concentration field to trigger the instability: where n is the wavenumber of the perturbation.

Nonlinear saturation
To quantify the long term growth of the instability, we examine how the kinetic An example time series of E k is shown in figure 4, which corresponds to the case of Sc = 1 and P e = 125. The result suggests that after the initial perturbation, there is a transient stage during which the kinetic energy grows exponentially, i.e. E k ∼ e 2st . As a consistency test, our simulation confirms that the involvement of non-linear terms in the simulation does not change the initial growth rate s. However, later the growth of E k begins to level off after some time because of non-linear saturation. Such non-linear saturation is common in most linearly unstable non-linear systems, such as Rayleigh-Bénard convection (Greenside & Coughran Jr 1984), Taylor-Couette flow for inner cylinder rotation (Grossmann et al. 2016), or Rayleigh-Taylor instability (Haan 1989). The concentration fields at different time show that during the saturated stage, the emitted plumes merge into a larger one, and eventually the flow structure develops into the state with a single large plume. Next, we compare the exponential growth rate s of the instability for various wavenumber cases (n = 1, 2, 3) during the initial stage with exponential growth as shown in figure  6. For the benchmark cases with only the linear terms, the data points (circles) agree excellently with the linear stability analysis (solid curves). This result can be regarded as further validation for our numerical code.
We further examine the situation with random initial perturbation solved with the full equations, including nonlinear terms. The above theoretical analysis has shown that for higher P e, the larger wavenumber mode can be triggered. The concentration fields in figure 7(a) provide more insight into the triggering of higher-order modes for larger P e. Different time instants of the concentration fields for different P e are shown in the figure.
For P e = 50, there is a single concentration plume generating initially. However, for larger P e = 125, multiple plumes are initially emitted. They undergo a merging process to form a single large plume. After formation of the single large plume, E k reaches the asymptotic value shown in the time series in figure 7(b). Interestingly, for even larger P e = 628, E k has spiky signals within a statistically steady state as shown in figure  7(c). The corresponding concentration fields in figure 7(a) reveals that small plumes are continuously generated from the reacting wall, and the merging of the plumes occurs simultaneously. Such continuous plume emission and merging can also clearly be seen from the Supplementary Movie.
Finally, we classify the four regimes based on the three criteria: • Growth rate of the instability.
• Number of plumes generated initially.
• Fluctuation of the kinetic energy (E k ) after reaching the statistically steady state. To quantify the number of generated plumes in the initial stage, we perform a Fourier transformation of the concentration field along the reacting wall (y = 0) at the instant when the plumes emerge (instant b in figure 4). The wavenumber, i.e., the initial number of plumes (circles), is compared with the dominant wavenumber as obtained from linear stability analysis (red dashed line) in figure 8(b). Both are in good agreement. The dominant wavenumber is that of the maximum growth rate s at a certain P e. Regarding  I  II  III IV Regime: Figure 8. The simulation result for the catalytic plane with non-linear terms and random initial perturbation for Sc = 1 and different P e: (a) Theoretical (dashed curve, which is equation (3.10)) and numerical result (circle) of s as a function of P e, which indicates that the system becomes unstable when P e > 8π. (b) Theoretical (dashed line, equation (3.11) without rounding operation) dominant wavenumber and numerical wavenumber n calculated by Fourier transform (circle) as a function of P e. The result indicates that when P e > 16π, multiple plumes are generated. (c) Standard deviation σ of the kinetic energy for different P e, which indicates that when P e 603, the kinetic energy eventually fluctuates because small plumes are continuously generated. Thus four regimes are classified, marked with different colors: stable (I: blue), a single wave (II: green), multiple waves which merge with each other (III: orange), and multiple waves with small plumes continuously being regenerated (IV: pink).
the fluctuation of E k , we evaluate the standard deviation of E k after reaching the statistical steady state in figure 8(c).
The four regimes are as follows: • Regime I (P e 8π): the system is stable.
• Regime II (8π < P e 16π): the system becomes unstable. Single plumes generate as can be seen in figure 8(b), and thus the dominant wavenumber is 1.
• Regime III (16π < P e 603): the initial wavenumber n becomes larger than one, and it increases with P e. The trigger of higher-order mode can be explained by the linear stability curve in figure 6. As P e > 16π, the perturbation of wavenumber n > 1 becomes unstable. For high enough P e, a higher wavenumber mode can grow even faster than the single wavenumber mode. After a while, the individual plumes merge into a single large one, and the system reaches an asymptotic state with constant E k (σ = 0 shown in figure 8(c)); • Regime IV (P e 603): the plume emission and merging happen continuously even after reaching statistically steady state, and therefore E k fluctuates with time (σ > 0). Figures 8(a) and (b) indicate that the exponential growth rate and the number of plumes generated initially can be approximately predicted by linear stability analysis. However, at high P e, there is a small deviation between the theory and our simulation. An explanation is that at high P e, various wavenumbers are excited simutaneously, such that the average growth rate becomes lower than the maximum growth rate predicted by linear stability analysis (equation (3.10)).

Dependence on Schmidt number
Based on the same classification criteria, we work out the full phase diagram in the (P e, Sc) parameter space, for 0.1 Sc 10. Figure 9 shows the four different regimes, namely the stable regime (I), the single plume regime (II), the multiple plume regime with a steady final state (III), and the regime with an unstable final state (IV). The transition points between the stable and the unstable regime (P e = 8π), and between the single plume and the multiple plume regime (P e = 16π) are insensitive to Sc. This can be understood from the linear stability analysis where the onset s for the nth wavenumber is P e = 8πn, independent of Sc. However, the onset of regime IV occurs at smaller P e, provided Sc < 1. When Sc 1, the onset P e of regime IV becomes independent of Sc.
To further understand why the onset of regime IV behaves differently for Sc < 1 and Sc 1, we have a close inspection on the event of the plume emission and merging for Sc = 0.1 and Sc = 1 shown in figure 10. First, for both cases when Pe is large enough, chaotic plume emissions are observed near the catalytic surface. However, the dynamics of the concentration plume are different for large and small Sc: For Sc = 1 as shown in figure 10 (a), the emitted small plumes gradually merge into the domain-sized plume, and this large plume is relatively stable. Thus, the velocity and concentration fluctuations are limited to near the vicinity of the catalytic surface without penetrating into the bulk region. In contrast, for Sc = 0.1 as shown in figure 10 (b), separate plumes merge and eventually be energetic enough to penetrate into the bulk, causing strong fluctuations in the bulk region.
To quantify this effect, we compute the fluctuation strength, once the system has reached the statistically steady state. It is characterized by the standard deviation of the horizontally averaged horizontal velocity u std (y): where u(x, y, t) is the instantaneous horizontal velocity and represents the average over time or x-direction, which is denoted by the subscript. The result is plotted in figure 11. From the figure, we find that the fluctuation is maximum at the bottom wall y = 0 since the diffusiophoretic flow at the wall drives the fluid flow. Moreover, for Sc < 1, the strong velocity fluctuations are not limited to the near-wall region, but also penetrate into the bulk. In contrast, for Sc 1, there are only large fluctuations in the near-wall region and u std (y) is monotonically decaying with wall distance y.
We now understand that the chaotic fluctuations observed in regime IV originate from different physical mechanisms for small and large Sc. For small Sc, as the fluctuations are mainly contributed from the bulk, one expects that the bulk viscous dissipation plays a role, and thus lower onset P e should be obtained for smaller Sc. However, it does not hold for the situation of large Sc since the fluctuations are mainly contributed by the chaotic plume emission close to the catalytic surface. To work out the details of the chaotic plume emission, non-linear stability is worthy to be conducted in the future. 16π Figure 9. The phase diagram for the case of the catalytic plane with different Sc and P e: For P e < 8π, the system is stable; for 8π < P e < 16π, the system becomes unstable and a single plume is generated; finally, for P e > 16π, multiple plumes are generated. For the last regime, there are two sub-regimes: for low P e (orange triangle), multiple plumes eventually merge to a single one and for higher P e (red circle), there is a newly found regime where the smaller plumes are continuously regenerated. The underlay colors are to guide the eyes.
As already mentioned in the introduction, our results share some similarities with those of the Bénard-Marangoni instability. For both cases, if the Péclet or Marangoni number is above a critical value, the system becomes unstable. Bergeon et al. (1998) comprehensively studied the Marangoni convection and found that as the Marangoni number increases, the plume will develop into single-roll or multiple-roll structures, which is similar to the single wavenumber or higher-order wavenumber modes observed in regime II and III, respectively.
As a final remark, Michelin et al. (2020) have shown the diffusio-phoretic instability in a confined phoretic channel, from which they also observe the generation of the plumes. Note that Michelin et al. (2020) have also considered the non-linear terms in the advection-diffusion equation, however, for the momentum equation, they consider the case of Stokes flow, such that the non-linear terms and effects of Schmidt number have not been considered. The chaotic plume emission observed in regime IV is the unique feature for high Péclet numbers, which, however, has not been focused on in most of the previous studies. Moreover, with an analytical calculation, we obtain the dominant wavenumber and its growth rate for different Sc and P e, which agrees with our simulation.

Simulation of the phoretic particle
Given the analysis of the catalytic plane, we now conduct three-dimensional simulations of a spherical phoretic particle to study the effect of plume emission and merging on the particle motion.

Numerical setup
The set-up is as follows: a phoretic particle is positioned at the center of the domain, and then due to diffusio-phoresis, the particle will self-propel. The governing equations consist of two parts. The first is the same as that in section 4.1, which is to solve the three-dimensional version of equations (2.3) and (2.4), except for the characteristic length which now is the radius of the particle. The second part involves the governing equation for the dynamics of the phoretic particle. However, one faces the challenge of dealing with a moving immersed boundary condition. To deal with it, we make use of moving least squares (MLS) based immersed boundary (IB) method, where the particle interface is represented by a triangulated Lagrangian mesh. For details of our MLS-based IB method, we refer to Spandan et al. (2017). The concentration boundary condition is that the wall normal concentration gradient is a constant ∂c ∂n = −1, which can be achieved by forcing the concentration at the particle surface based on the concentration interpolated at the probe located at a short distance (1 grid size) from the surface of the particle. The velocity boundary condition is where u s is the surface gradient (∇ s ) of the concentration. The domain size is L x × L y × L z = 20R × 20R × 40R, in terms of the particle radius R. We use uniform grids N x × N y × N z = 201 × 201 × 401. Mesh refinement tests are done at P e = 15, 16 and 20 with doubled grid numbers in each dimension. Figure 12(a) indicates that the result for the grid 401 × 401 × 801 is nearly indistinguishable from that for 201 × 201 × 401.
For the spherical particle, the radius R is used as the length scale in Péclet number: We will present the result of phoretic particles for different P e from 3 to 20 with Sc = 1.

Result of the phoretic particle
Similar to the case of the catalytic plane, the diffusio-phoretic instability breaks the rotational-symmetry of the phoretic particle. It has been shown by Michelin et al. (2013)  Sc=0.4 Figure 11. The standard deviation u std (y) as function of the wall distance y for different Schmidt numbers with Péclet number P e = 754. Averaging was done of time and over the x-direction. All cases belong to regime IV.
that the phoretic particle breaks the symmetry when P e is larger than 4. Therefore, as a validation, we first simulate cases with small P e and compare to the results obtained from Michelin et al. (2013). In figure 12(a), we plot the numerical terminal velocity U ∞ . For P e > 4, indeed symmetry breaking occurs (e.g. figure 13 (a) for P e = 10) and the particle moves along a straight trajectory. The terminal velocities agree with those obtained from Michelin et al. (2013). Furthermore, in figure 12(b) we check whether the terminal velocity is sensitive to the Schmidt number Sc. Interestingly, the figure suggests that when Sc 1 for P e = 8 and 12, the terminal velocity converges to a constant. However, to understand why terminal velocity levels off, further study is needed in the future. Regarding the grid resolution requirement, the necessary grid resolution will increase dramatically for very large Sc. In order to run as many cases as possible to fully explore the phase diagram, we stick to Sc = 1 for the rest of our simulations. After this validation we now extend the calculations to higher P e. Multiple plumes emission and merging occur at the surface of the phoretic particle (e.g. figure 13 (b) for P e = 60), which is similar to that observed for the catalytic plane. The continuously emitted plumes change the direction of the phoretic particle and lead to chaotic motion.
To characterize the motion of the particle, we calculate the mean temporal autocorrelation of the particle direction: where e u (t) = u(t)/|u(t)| is the unit direction vector of the particle velocity at t. The integral upper limit T is chosen large enough to achieve statistical stationary. The autocorrelation for different P e is shown in figure 12 (c). When P e = 20, the auto-correlation becomes considerably less than 1 with increasing ∆t, which means that the particle starts to meander in different directions. The chaotic behavior of the particle has also been shown by the fluctuation of velocity in figure 12(a), which is represented by the solid bars. Interestingly, for P e 15, the average velocity (the red circle) still lies near the result by Michelin et al. (2013), but for larger P e, the velocity shows larger fluctuation, i.e. more chaotic behavior.
Thus we can classify the motion for a phoretic particle into three regimes: The motion of the phoretic particle is divided into three different regimes: stable, symmetry breaking, and chaotic motion due to plume generation. (b) The nomalized terminal velocity of different Sc for the case P e = 8 and 12. The velocity is normalized by the terminal velocity at Sc = 40. The result shows that when Sc > 1, the terminal velocity converges to a constant. (c) The temporal auto-correlation function of the unit direction vector for three different P e for Sc = 1. The temporal auto-correlation indicates whether the particle performs chaotic motion or not.
• 4 < P e 15, symmetry breaking occurs and the particle moves straight.
• P e 15, the particle moves chaotically. Figure 13 (b) shows that the plumes are continuously generated and merge, which alters the concentration distribution and steer the moving direction of the phoretic particle. This shows that our newly found regime IV in figure 9 leads to the chaotic motion for the case of phoretic particles.

Comparison between phoretic particle and catalytic plane
We now compare the various regimes for the catalytic plane (Section 4) with those for the phoretic particle (Section 5). The similarity between the two setups is that both the instability and chaotic flow can be observed for both setups. A major difference between them is that the regimes for the catalytic plane are classified by the distinct plume dynamics whereas the regimes for the phoretic particle are classified by the distinct particle motions. However, both regimes II and III for the catalytic plane lead to a steady final state with a single plume, which for the case of the phoretic particle implies t=160 t=200 Figure 13. The concentration cross-section from three-dimensional simulations of an isotropic catalytic particle for P e = 10 and 60. Again, Sc = 1. The simulation is at a domain Lx ×Ly ×Lz = 20R×20R×40R, in terms of the particle radius R. The grids are 201×201×401.
To better demonstrate the chaotic trajectory, the motion of the particle in (b) is projected to x-z plane. (a) P e = 10, the particle moves straightly. (b) P e = 60, plumes are generated at the surface of the particle, which starts to move irregularly.
straight motion. For classification of the particle motion, the same fate of particle motion leads one to define only one regime despite the distinct plume dynamics during the initial transient stage.
Another difference between the catalytic plane and phoretic particle is that the onset Péclet numbers are different. However, this difference simply reflects that the characteristic length scales in the definition of P e are different for both systems.
We note that also Hu et al. (2019) numerically observed the chaotic motion of phoretic particles. However, the plume generation and merging, which could provide a route to understand the chaotic motion of the particle at high enough P e, were not studied in that paper. Besides, in the experiments of active droplets, it is also observed that the droplet can move in a helical or even chaotic trajectory at high P e (Suga et al. 2018;Maass et al. 2016). Morozov & Michelin (2019b) also observed the helical and chaotic motion of the catalytic particle.Recently, the stochastic dynamics of active particles was analysed (Gaspard & Kapral 2018;Chamolly & Lauga 2019). Based on the stochastic approach using Langevin equations, the active particle motion is split into a diffusive part and a ballistic part. In our work here, it is the deterministic plume emission that is the source of the diffusive motion, and an one-to-one comparison is difficult due to the quite different natures of the approaches. Very recently, Vajdi Hokmabad et al. (2021) observed plume generation and merging at the surface of a meandering chemically active droplet. This recent finding reflects the importance of the plume dynamics in determining the droplet motion. Here, for the diffusio-phoretic particle, we have also revealed such plume generation and merging phenomena and have related it to the instability of the flow near the surface.

Concluding remarks
In summary, we have studied the instability driven by diffusio-phoretic effects at the interfaces for two different systems: a catalytic plane and a spherical phoretic particle. The Péclet number (P e) and Schmidt number (Sc) are the parameters that determine the states of the system.
For a catalytic plane, via linear stability analysis, we quantitatively studied the growth of various wavenumber perturbation. With assistance of the simulation, we have classified four regimes for different P e and Sc based on the exponential growth rate of the instability, number of plumes generated initially and fluctuation of the kinetic energy after reaching the statistically steady state (E k ). For P e 8π, the system is stable. For 8π < P e 16π, the system becomes unstable, a single plume is generated and the system reaches a steady state eventually. For P e > 16π, multiple plumes are generated initially, which merge into a single one to attain a stable state eventually due to nonlinear saturation. However, for even higher P e (P e 603 for Sc = 1), small plumes are continuously generated and merge with each other, the system remains unstable and therefore E k fluctuates in time.
Based on the linear stability analysis, we understand that the onset P e between regime I and II, and regime II and III are independent of Sc. However, there is noticeable effect of Sc on the transition to regime IV, which is associated with different flow structures for different Sc. For small Sc, the strong fluctuations of concentration and kinetic energy also occur in the bulk region, whereas for large Sc, the fluctuations are only contributed by the chaotic plume emission close to the catalytic plane. As the viscous dissipation in the bulk plays a role for small Sc cases, lower onset Péclet numbers of regime IV are obtained for lower Sc. However, it does not hold for large Sc.
Then we extend our research to three-dimensional simulations of the spherical phoretic particle. Despite the geometric difference, an analogous phenomenon happens at the surface of the particle which triggers different particle motions. Similar to the case of the catalytic plane, for the case at Sc = 1, when P e > 4, the particle starts to break the symmetry. For higher P e 15, also similar to the observation of the catalytic plane, the small plumes start to be generated continuously at the surface of the particle, which will steer the particle and lead to meandering motion. The analogous phenomenon indicates that the chaotic motion of the phoretic particle results from the instability at the interface driven by diffusio-phoretic effects.
The present work makes a contribution to the understanding of the diffusio-phoretic instability. First, the study reveals the existence of a highly unstable regime at high P e. We not only study the onset P e of the unstable mode, but also analytically work out the dominant wavenumber as a function of P e. For high enough P e (regime III), multiple plumes are emitted into the surrounding fluid. For even higher P e (regime IV), smaller wavelength perturbations are dominant, which leads to continuous plume generation and subsequently to chaotic flow. Second, our results show that the diffusio-phoretic instability at the catalytic surface can eventually lead to chaotic motion of the phoretic particle. Through simulations of the phoretic particle at high P e, we not only see its chaotic motion, but also observe the plume emission and merging events near the surface of the particle, which is similar to the situation of regime IV for the case of the catalytic plane. The study of the phoretic plane thus provides a framework to understand the motion of the phoretic particle.
Many questions remain open. For example, how does the particle motion and flow field change for phoretic particles in a complicated environment, such as phoretic particle near a wall? How does the plume generation and merging change the collective behavior of phoretic particles? How about the effect of plume generation on rod particles rather than spheres? Building on the here obtained insight into the mechanism behind the chaotic motion of phoretic particles, it is worthwhile to further explore the effects of the plume generation on the motion of particles in the more complicated setups as mentioned above, in particular, on collective effects. The solution of (A 6) that vanishes at infinity is Recall that (Gradshteyn & Ryzhik (2007), p. 1110) and use the convolution theorem. Then the concentration solution in physical space can be obtained as (Wu et al. (2006), p. 144), The corresponding derivative is ∂c ∂y = − ∞ ξ0 2 √ π e −ξ 2 dξ, with ξ 2 = P ey 2 4(t − τ ) , ξ 2 0 = P ey 2 4t .

(A 21)
Thus, for given k and Sc, we obtain the curve s vs P e in Figure 3(b) by varying δ. Similarly, for given k, we obtain the contour s in the (P e, Sc) plane in Figure 3(a) by varying δ and Sc.

A.3. Determination of the dominant wavenumber
The maximum growth rate, as well as the dominant wavenumber, can also be determined from (3.