Active Brownian motion of emulsion droplets: Coarsening dynamics at the interface and rotational diffusion

A micron-sized droplet of bromine water immersed in a surfactant-laden oil phase can swim (S. Thutupalli, R. Seemann, S. Herminghaus, New J. Phys. 13 073021 (2011)). The bromine reacts with the surfactant at the droplet interface and generates a surfactant mixture. It can spontaneously phase-separate due to solutocapillary Marangoni flow, which propels the droplet. We model the system by a diffusion-advection-reaction equation for the mixture order parameter at the interface including thermal noise and couple it to fluid flow. Going beyond previous work, we illustrate the coarsening dynamics of the surfactant mixture towards phase separation in the axisymmetric swimming state. Coarsening proceeds in two steps: an initially slow growth of domain size followed by a nearly ballistic regime. On larger time scales thermal fluctuations in the local surfactant composition initiates random changes in the swimming direction and the droplet performs a persistent random walk, as observed in experiments. Numerical solutions show that the rotational correlation time scales with the square of the inverse noise strength. We confirm this scaling by a perturbation theory for the fluctuations in the mixture order parameter and thereby identify the active emulsion droplet as an active Brownian particle.

There are various methods to construct a microswimmer. One idea is to generate a slip velocity field close to the swimmer's surface using a phoretic mechanism. A typical example of such an artificial swimmer is a micron-sized spherical Janus colloid, which has an inherent polar symmetry. Its two faces are made of different materials and thus differ in their physical or chemical properties [32]. For example, a Janus particle with faces of different thermal conductivity moves if exposed to heat. The conversion of thermal energy to mechanical work in a self-generated temperature gradient is called self-thermophoresis [33]. Janus colloids also employ other phoretic mechanisms to become active [34,35,36,37].

Send offprint requests to:
A different realization of a self-propelled particle is an active emulsion droplet. The striking difference to an active Janus particle is the missing inherent polar symmetry. Instead, the symmetry between front and back breaks spontaneously, for example, in a subcritical bifurcation [38]. The self-sustained motion of active droplets is due to a gradient in surface tension, which is usually caused by an inhomogeneous density of surfactants. The resulting stresses set up a solutocapillary Marangoni flow directed along the surface tension gradient that drags the droplet through the fluid. An active droplet generates a flow field in the surrounding fluid typical for the "squirmer" [39,40,41,42,43]. Originally, the squirmer was introduced to model the locomotion of microorganisms that propel themselves by a carpet of short active filaments called cilia beating in synchrony on their surfaces. The squirmer flow field at the interface is then a coarse-grained model of the cilia carpet.
An active droplet, which swims due to solutocapillary Marangoni flow, has recently been realized [1]. Water droplets with a diameter of 50 − 150µm are placed into a surfactant-rich oil phase. The surfactants migrate to the droplet interface where they form a dense monolayer. Bromine dissolved in the water droplets reacts with the surfactants at the interface. It saturates the double bond in the surfactant molecule and the surfactant becomes weaker than the original one. Hence, the "bromination" reaction locally increases the interfacial surface tension. This induces Marangoni flow, which advects surfactants and thereby further enhances the gradients in surface tension. If the advective current exceeds the smoothing diffusion current, the surfactant mixture phase-separates. The droplet develops a polar symmetry and starts to move in a random direction, which fluctuates around such that the droplet performs a persistent random walk. While the droplet swims with a typical swimming speed of 15µm/s, brominated surfactants are constantly replaced by nonbrominated surfactants from the oil phase by means of desorption and adsorption. Finally, the swimming motion comes to an end when the fueling bromine is exhausted.
In ref. [38] we developed a diffusion-advection-reaction equation for the surfactant mixture at the droplet interface and coupled it to the axisymmetric flow field initiated by the Marangoni effect. In a parameter study we could then map out a state diagram including the transition from the resting to the swimming state and an oscillating droplet motion. In this paper we combine our theory with the full three-dimensional solution for the Marangoni flow, which we derived for an arbitrary surface tension field at the droplet interface in ref. [43]. Omitting the constraint of axisymmetry and adding thermal noise to the dynamic equation of the surfactant mixture, we will focus on two new aspects of droplet dynamics that we could not address in ref. [38]. First, while reaching the stationary uniaxial swimming state, the surfactant mixture phase-separates into the two surfactant types. We illustrate the coarsening dynamics and demonstrate that it proceeds in two steps. An initially slow growth of domain size is followed by a nearly ballistic regime. This is reminiscent to coarsening in the dynamic model H [62]. Second, even in the stationary swimming state the surfactant composition fluctuates thermally and thereby initiates random changes in the swimming direction, which diffuses on the unit sphere. As a result the droplet performs a persistent random walk, as observed in experiments [1], which we will characterize in detail.
The article is organized as follows. In sect. 2 we recapitulate our model of the active emulsion droplet from ref. [38] and generalize it to a droplet without the constraint of axisymmetry. While sect. 3 explains the numerical method to solve the diffusion-advection-reaction equation on the droplet surface, the following two sections contain the results of this article. Section 4 describes the coarsening dynamics of the surfactant mixture before reaching the steady swimming state and sect. 5 characterizes the per-sistent random walk of the droplet in the swimming state. The article concludes in sect. 6.

Model of an active droplet
In order to model the dynamics of the active droplet, we follow our earlier work [38]. We use a dynamic equation for the surfactant mixture at the droplet interface that includes all the relevant processes. We assume that the surfactant completely covers the droplet interface without any intervening solvent. We also assume that the head area of both types of surfactant molecules (brominated and non-brominated) is the same. Denoting the brominated surfactant density by c 1 and the non-brominated density by c 2 , we can therefore set c 1 + c 2 = 1. We then take the concentration difference between brominated and non-brominated surfactants as an order parameter φ = c 1 − c 2 . In other words φ = 1 corresponds to fully brominated and φ = −1 to fully non-brominated surfactants and c 1 = (1 + φ)/2 and c 2 = (1 − φ)/2. Finally, we choose a constant droplet radius R.

Diffusion-advection-reaction equation
The dynamics of the order parameter φ at the droplet interface can be expressed as [38]: which we formulate in the form of a continuity equation with an additional source and thermal noise (ζ) term. ∇ s = (1 − n ⊗ n)∇ stands for the directional gradient on a sphere with radius R, where ∇ is the nabla operator and n the surface normal. The current is split up into a diffusive part j D and an advective part j A , which arises due to the Marangoni effect. We summarize them below and in sect. 2.2. The source term describes the bromination reaction as well as desorption of brominated and adsorption of nonbrominated surfactants to and from the outer fluid. Both processes tend to establish an equilibrium mixture with order parameter φ eq during the characteristic relaxation time τ R . Ad-and desorption dominate for φ eq < 0 while bromination dominates for φ eq > 0. The source term is a simplified phenomenological description for the ad-and desorption of surfactants. A more detailed model would include fluxes from and to the bulk fluid [63]. We will explain the thermal noise term further below. The general mechanism of eq. (1) to initiate steady Marangoni flow is as follows. The diffusive current j D smoothes out gradients in φ, while the advective Marangoni current j A amplifies gradients in φ. Hence, j D and j A are competing and as soon as j A dominates over j D , φ experiences phase separation. As a result, the resting state becomes unstable and the droplet starts to swim.
We now summarize features of the diffusive current j D , more details can be found in ref. [38]. We formulate a Flory-Huggins free energy density in terms of the order parameter of the surfactant mixture, which includes entropic terms and interactions between the different types of surfactants: Here, ℓ 2 is the head area of a surfactant at the interface. We introduce dimensionless parameters b 1 (b 2 ) to characterize the interaction between brominated (non-brominated) surfactants and b 12 describes the interaction between the two types of surfactants. The diffusive current is now driven by a gradient in the chemical potential derived from the total free energy functional F [φ] = f (φ) dA: where the Einstein relation D = λk B T /ℓ 2 relates the interfacial diffusion constant D to the mobility λ. To rule out a double well form of f (φ), which would generate phase separation already in thermal equilibrium, we only con- This also means that the diffusive current j D ∝ −∇ s φ is for all φ indeed directed against ∇ s φ. In the following we assume b 12 = (b 1 + b 2 )/2 and therefore require b 1 + b 2 < 4.
We formulate the thermal noise term in eq. (1) as Gaussian white noise with zero mean following ref. [64]: Here, the strength of the noise correlations is connected to the mobility λ of the diffusive current via the fluctuationdissipation theorem. In order to close eq. (1), we now discuss the advective Marangoni current j A .

Marangoni flow
The advective current for the order parameter φ is given by where u| R is the flow field at the droplet interface. It is driven by a non-uniform surface tension σ and therefore called Marangoni flow [65,63]. In our case, we have a nonzero surface divergence ∇ s · u| R = 0. In fact, it can be shown that an incompressible surface flow cannot lead to propulsion of microswimmers [66]. In order to evaluate u| R , one has to solve the Stokes equation for the flow field u(r) surrounding the spherical droplet (r > R) as well as for the flow fieldû(r) inside the droplet (r < R). Both solutions are matched at the droplet interface by the condition [63], where P s = 1 − e r ⊗ e r is the surface projector. Equation (5) means that a gradient in surface tension σ is compensated by a jump in viscous shear stress. Here, T = η[∇ ⊗ u + (∇ ⊗ u) T ] is the viscous shear stress tensor of a Newtonian fluid with viscosity η outside of the droplet and the same relation holds forT of the fluid with viscosityη inside the droplet. We have performed this evaluation in ref. [43] for a given surface tension field and only summarize here the results relevant for the following. Alternative derivations are found in ref. [67,68,69,70].
In spherical coordinates the Marangoni flow field u| R at the interface reads [43,67,68,69] with spherical harmonics Y m l (θ, ϕ) given in appendix A. Here, are the expansion coefficients of the surface tension, where Y m l means complex conjugate of Y m l , and [43,68,70] v is the droplet velocity vector. It is solely given by the dipolar coefficients (l = 1) of the surface tension and determines propulsion speed v D ≥ 0 as well as the swimming direction e with |e| = 1. Note that by setting m = 0, eqs. (6)-(8) reduce to the case of an axisymmetric droplet swimming along the z-direction, as studied in ref. [38]. In ref. [43] we give several examples of flow fields u| R . In general, Marangoni flow is directed along gradients in surface tension, i.e. u| R ∇ s σ. This is confirmed by eq. (6) and also clear from fig. 2 (b), which we discuss later. However, according to eq. (6) higher modes of surface tension contribute with a decreasing coefficient [43]. Note the velocity field in eq. (6) is given in a frame of reference that moves with the droplet's center of mass but the directions of its axis are fixed in space and do not rotate with the droplet. Finally, the velocity fields inside (û) and outside (u) of the droplet in both the droplet and the lab frame can be found in the appendix of ref. [43].
The surface tension necessary to calculate v D and u| R is connected to the order parameter φ by the equation of state, σ = f − ∂f ∂c1 c 1 − ∂f ∂c2 c 2 , which gives [38] (9) This implies that for b 1 > b 2 > 0, ∇ s φ points along ∇ s σ. Moreover, since the Marangoni flow u| R is oriented along ∇ s σ, as noted above, we conclude that for φ > 0 the advective current j A = φu| R points "uphill", i.e., in the direction of ∇ s φ, in contrast to j D [38].
This completes the derivation of the surface flow field u| R as a function of the expansion coefficients s m l of the surface tension. Together with the equation of state σ(φ) the advective current j A in eq. (4) is specified. Finally, using the diffusion current j D from eq. (2), the diffusionadvection-reaction equation (1) becomes a closed equation in φ.
The swimming emulsion droplet is an example of a spherical microswimmer, a so-called squirmer [39,40,41,42,43]. Squirmers are often classified by means of the socalled squirmer parameter β [5]. When β < 0, the surface flow dominates at the back of the squirmer, similar to the flow field of the bacterium E. coli. Since such a swimmer pushes fluid outward along its major axis, it is called a 'pusher'. Accordingly, a swimmer with β > 0 is called a 'puller'. The algae Chlamydomonas is a biological example of a puller. Swimmers with β = 0 are called 'neutral'.
For an axisymmetric emulsion droplet swimming along the z-direction, the squirmer parameter is given by with coefficients s m l from the multipole expansion (7) of the surface tension σ [43]. A generalization of this formula to droplets without axisymmetry and swimming in arbitrary directions is derived in ref. [43]. The relevant expressions are presented in appendix B.

Reduced dynamic equations and system parameters
In order to write eq. (1) in reduced units, we rescale time by the characteristic diffusion time τ D = R 2 /D and lengths by droplet radius R, and arrive at (11) where the Gaussian noise variable fulfills The dimensionless velocity field at the interface and the droplet velocity vector read, respectively, All quantities in eqs. (11) and (13), including j D , u| R , t, ∇ s , ζ, and v D , are from now on dimensionless, although we use the same symbols as before. Writing the dynamics equations in reduced units, introduces the relevant system parameters M, ν, κ, φ eq , and ξ, which we discuss now. The Marangoni number M quantifies the strength of the advective current in eq. (11) and is given by M = (b1−b2)R λ(η+η) . It is the most important parameter of our model, as it determines whether the droplet swims. In eq. (13a) we introduced the ratio of shear viscosities, ν =η/η, for the fluids inside and outside of the droplet, respectively. In our study we consider a water droplet suspended in oil and set ν ≈ 1/36 [1]. The interaction parameters b 1 and b 2 not only appear in M but also as b 1 + b 2 in the diffusive current in eq. (2) and in the equation of state σ(φ) in eq. (9). Therefore, they need to be set individually. Assuming the head area of a surfactant ℓ 2 to be on the order of nm 2 , we can fit eq. (9) to the experimental values σ(φ = 1) ≈ 2.7mN/m and σ(φ = −1) ≈ 1.3mN/m[1] to find b 1 ≈ 0.6 and b 2 ≈ 0.3. We keep these values fixed throughout the article.
Parameter κ = τ D /τ R tunes the ratio between diffusion and relaxation time and the equilibrium order parameter φ eq measures whether ad-and desorption of surfactants (φ eq < 0) or bromination (φ eq > 0) dominates. In this study we set κ = 0.1 and φ eq = 0.5. A parameter study for these parameters can be found in [38]. Finally, the reduced noise strength ξ = ℓ/R ∝ 1/ √ N , where N is the total number of surfactants at the droplet interface, connects the the droplet size R to the molecular length scale ℓ.
The following sect. 3 describes, how we solved the dynamic equation (11) numerically. Readers not interested in the details can proceed immediately to sect. 4, where we present our first results.

Finite volume method on a sphere
To numerically solve the rescaled dynamic equation (11) for the order parameter field φ, we had to decide on an appropriate method. The most widely used numerical methods for solving partial differential equations are the finite difference method (FDM), the finite element method (FEM), and the finite volume method (FVM) [71,72]. We ruled out FDM due to numerical complications of its algorithm with spherical coordinates. They are most appropriate for the spherical droplet surface but one needs to define an axis within the droplet. The FEM is also very delicate when writing a numerically stable code for our model. This is mainly due to the advective term in eq. (11), which commonly causes difficulties in FEM routines [71]. In contrast, the FVM is especially suited for solving continuity equations. Therefore, it is much more robust for field equations that incorporate advection and we chose it for solving eq. (11) on the droplet surface.
In order to generate a two-dimensional FVM mesh that is as uniform as possible and quasi-isotropic on a sphere, we chose a geodesic grid based on a refined icosahedron [73]. An icosahedron has f 0 = 20 equilateral triangles as faces and v 0 = 12 vertices. In each refinement step, each triangle is partitioned into four equilateral triangles and the three new vertices are projected onto the unit sphere enclosing the icosahedron. Hence, after the n-th refinement step, the resulting mesh has f n = 4 n f 0 triangular faces and v n = v n−1 + 3 8 4 n f 0 grid points. 1 The "finite volume" then refers to a small volume (in this case an area) surrounding each grid point of the mesh. Thus we have to construct the Voronoi diagram of the triangular mesh. The Voronoi diagram consists of v n elements, 12 of which are pentagons associated with the vertices of the original icosahedron while the rest are hexagons. Unless otherwise noted we use a Voronoi mesh with v 3 = 642 FVM elements. The geodesic icosahedral grid is a standard grid in geophysical fluid dynamics. A comprehensive review on numerical methods in geophysical fluid dynamics can be found in [74].
In the following we will outline how we convert the diffusion-advection-reaction equation (11) to a set of ordinary differential equations for a vector φ comprising the values φ i of the order parameter field at the center points of all FVM elements. FVM was developed for treating current densities in a continuity equation and we illustrate the procedure for the diffusion term of eq. (11). We start by integrating over element i with area A i and use the divergence theorem, where n i is the outward normal at the element boundary: In the last term of eq. (14a), the line integral is converted into a sum over the N straight element boundaries of length l ij and n ij is the normal vector at the corresponding boundary. Figure 1 illustrates the relevant quantities.
In the second line the directional derivative n ij · ∇ s φ resulting from j D in eq. (2) is approximated by a difference quotient. The prefactor in j D , which we abbreviated by D(φ i , φ j ) in eq. (14b), also contains φ. It is interpolated at the boundary between elements i and j by means of the central differencing scheme as (φ i + φ j )/2. Finally, we write the whole term as the product of local diffusion matrix D i and vector φ. After applying this technique to all elements, the matrices D i are combined into one matrix D for the whole mesh. The same procedure is carried out for the advective term in eq. (11) but discretizing j A = M φu| R needs more care. While u| R is directly calculated at the boundary between elements i and j, the order parameter φ is treated differently. If the local Peclet number Pe = h ij M |u| R |/D(φ i , φ j ) is larger than 2, the central differencing scheme fails to converge. Instead a so-called upwind scheme is used, which takes into account the direction of flow [71]. For outward oriented flow, i.e. u| R · n ij > 0, one uses the element order parameter φ i , while for inward flow, i.e. u| R · n ij < 0, one uses the order parameter of the neighboring element φ j . In the case Pe < 2, φ is interpolated by the central Finally, the linear terms in φ and its time derivative are simply approximated by φ i andφ i . In the end, we are able to write the discretized eq. (11) as a matrix equation for the vector φ: (15) where the diagonal matrix M carries the areas of the elements, and with diffusion matrix D, advection matrix A, and element noise vector z, which describes typical Gaussian white noise with zero mean and variance one, In appendix C we derive eq. (16b) by integrating eq. (12) over two FVM elements i and j. Finally, the set of stochastic differential equations are integrated in time by a standard Runge-Kutta scheme.
In the following we present results obtained with the described numerical scheme.

Dynamics towards the swimming state
This section focuses on the dynamics of the active emulsion droplet from an initial resting state with swimming speed v D = 0 to a stable swimming state with swimming speed v D > 0. After a comparison with the axisymmetric model of the droplet from our previous work [38], where we also did not include thermal fluctuations, we investigate the coarsening dynamics of the order parameter φ at the droplet interface while reaching the swimming state.

Swimming speed v D
In order to test the simulation method, we start our analysis with a set of parameters, for which we found a swimming state in the inherent axisymmetric model [38]. They are given by Marangoni number M = 3, reduced reaction rate κ = 0.1, and equilibrium order parameter value φ eq = 0.5. We keep these values fixed throughout the following unless otherwise noted. The initial condition for solving eq. (15) is an order parameter field that fluctuates around φ eq : φ(θ, ϕ) = φ eq + δφ(θ, ϕ). The small fluctuations δφ(θ, ϕ) ≪ 1 are realized by random numbers drawn from the normal distribution N (φ eq , α 2 ) with mean φ eq   [38]. Note that the Marangoni flow u|R is directed along the gradients of φ and surface tension σ. and variance α 2 = 10 −5 and added at the grid points of the simulation mesh. Furthermore we set the noise strength to ξ = 10 −3 . Figure 2(a) shows the droplet swimming speed v D as a function of elapsed time.
First of all, we notice the good agreement with the corresponding graph of v axi D of the axisymmetric system of ref. [38], which we also plot in fig. 2(a). The same applies to the order parameter profile φ and the surface velocity field u| R of the swimming state, when averaged about the swimming axis e, see fig. 2(b). Thus, the full three- dimensional description presented in this work is consistent with the axisymmetric model of ref. [38]. The same is true for the squirmer parameter β from eq. (28), for which we find β ≈ −1.2 for M = 3. This is fairly close to the value of the axisymmetric model (β ≈ −0.8) and confirms that the swimming active droplet is a pusher. We stress that the Marangoni number M is the crucial parameter in our model, as it determines whether the droplet rests or swims. For small M , the homogeneous state φ = φ eq is stable, i.e., any disturbance δφ of the initially uniform φ is damped by the diffusion and reaction terms of eq. (11). As a result, the droplet rests. The transition to the swimming state occurs at increasing Marangoni number M via a subcritial bifurcation as illustrated in fig.  3, which shows swimming speed v D and squirmer parameter β plotted versus M . We use here a system without thermal noise, i.e., ξ = 0, in order to monitor the complete transition region of the subcritical bifurcation. At a transition value M tr the advective term of eq. (11) overcomes the damping terms. The homogeneous state becomes unstable and the droplet starts to swim with a finite swimming speed v D . As usual for a subcritical bifurcation, the transition to the swimming state takes place in a finite interval of M . There, the transition Marangoni number M tr depends on the initial disturbance strength α 2 of the uniform order parameter profile. The inset of fig. 3 confirms this statement. Next, we will discuss the biaxial evolution and the coarsening dynamics of the order parameter field, which we could not study in the axisymmetric description.

Transient biaxial dynamics
The good agreement of the rotationally averaged order parameter profile φ ϕ and the axisymmetric φ axi from our earlier work, both plotted in fig. 2 (b), suggests that in the steady swimming state, the full three-dimensional solution is also nearly axisymmetric about the swimming axis e. However, in non-steady state we expect φ to deviate from axisymmetry, which we quantify by introducing an appropriate measure for the biaxiality of the order parameter field φ. In analogy to characterizing the orientational order of liquid crystals, we define for the order parameter profile the traceless quadrupolar tensor [75] with surface normal n, unit tensor ½, and the surface integral is performed over the whole droplet interface. Just as in the case of the moment of inertia tensor, the eigenvalues and eigenvectors of Q characterize the symmetries of the order parameter field φ. If two eigenvalues of Q are equal, φ is said to be uniaxial. On the other hand, if all eigenvalues of Q are distinct, φ is biaxial. Finally, the case of three vanishing eigenvalues, i.e., Q = 0, describes an isotropic or uniform order parameter field φ or at least with tetrahedral or cubic symmetry. A measure for the degree of biaxiality, which incorporates the three mentioned cases, is given by the biaxiality parameter [76,77] If the order parameter field φ is axisymmetric or isotropic, ∆ = 0, while with increasing biaxiality ∆ approaches 1. In fig. 2 (a), we plot ∆ as a function of time. At the initial time t = 0, the order parameter profile is roughly uniform with ∆ ≈ 0 (not visible). As the droplet speeds up, the biaxiality parameter ∆ fluctuates strongly between 0 and 1. Starting at t ≈ 3, ∆ sharply decreases towards zero before the swimming speed becomes maximal. Finally, in the steady swimming state, ∆ is nearly zero but still fluctuates due to the thermal noise in the order parameter profile φ, which we indicate by the error bars in fig. 2 (b). Hence, during the speed up of the droplet, the order parameter field φ clearly is not axisymmetric.

Coarsening dynamics
The period of strong biaxiality goes in hand with the coarsening dynamics of the order parameter profile towards steady state. Figure 4(a) shows the order parameter profile φ(θ, ϕ) at various time steps for the same simulation run as in fig. 2. Shortly after the simulation starts with the nearly uniform initial condition, small islands or domains with φ > φ eq and φ < φ eq emerge, which rapidly grow until t ≈ 1, where the droplet hardly moves, see fig.  2 (a). Then the coarsening or demixing process is slowed down. The domains coalesce on larger scales and the droplet speeds up significantly. Since the droplet interface area is finite, the domains coalesce at some point to one large region which covers about half of the interface. From then on the droplet interface is covered by only two regions with φ < φ eq and φ > φ eq . The domain wall between the two regions is close to the equator and the droplet has reached its top speed [compare v D (t ≈ 5) in fig. 2(a)]. Then, the domain wall moves towards the southern pole to its final position. Since its circumference shrinks, the droplet speed v D slows down to its stationary value, which it reaches at t ≈ 9. Thus, overshoot of the swimming velocity in fig. 2(a) is the result of two processes taking place on different time scales: the coarsening process and the final positioning of the domain wall at a somewhat larger time scale, which depends on the parameters.
Note that depending on the final position of the domain wall separating the two regions, the droplet is either a pusher or a puller. If the domain wall with increasing φ is situated in the southern hemisphere (π/2 < θ < π), the droplet is a pusher. If it is located in the northern hemisphere (0 < θ < π/2), a puller is realized. However, in our simulations the swimming droplet is always a pusher irrespective of the system parameters. This is due to the fact that the advective Marangoni current j A at the interface of the swimming droplet is always directed towards the southern hemisphere. This flow also moves the domain wall away from the equator and towards the posterior end of the droplet. The squirmer parameter β varies in the range −2 < β < 0 depending on Marangoni number M (see fig. 3) and equilibrium order parameter φ eq . This is in agreement with earlier observations in ref. [38]. The timeframe t > 9, where the swimming speed fluctuates around its steady-state value, will be covered in sect. 5.
To quantify further the spatial structure of the order parameter profile during coarsening, we examine the angular power spectrum |s m l | 2 of the surface tension. It is related to φ in eq. (9). Using the orthonormality relation of spherical harmonics Y m l (θ, ϕ), given in appendix A, one can compute the total power P of the surface tension σ: Here, the polar power spectrum g l characterizes the variation of the surface tension and thus the order parameter field φ along the polar angle θ. In particular, g l for small l quantifies the large-angle variations of σ. Note that g 1 is directly related to the swimming speed v D calculated from eq. (8) in the polar coefficients s m 1 . Using s −1 Figure 4(b) depicts the polar power spectrum g l normalized by the total power P at the same time steps of the coarsening dynamics discussed before in Fig. 4(a). We also show an ensemble average of g l /P . At the initial time t = 0, the spectrum of g l is solely characterized by frequencies or polar contributions of the noisy initial condition φ(t = 0) = φ eq + δφ. Thus, the maximum frequency or polar number l of the spectrum at t = 0 is set by the level of refinement of the simulation mesh. During the initial period of fast coarsening until t = 1, the polar power spectrum shifts from high to low frequencies indicating the increase of domain sizes. Then the higher frequencies vanish more and more from the spectrum, as the phases associated with φ < φ eq and φ > φ eq separate. Eventually, the spectrum g l strongly peaks at l = 1 while the remaining coefficients become insignificant in comparison. Finally, from t = 5 to t = 9, the first coefficient g 1 of the angular power spectrum decreases again while the second and third coefficients g 2 and g 3 rise. This confirms that in the final stage the droplet slows down its velocity v D and tunes its squirmer parameter β by shifting the domain wall further away from the equator.
In order to quantify further the temporal evolution of the coarsening dynamics, we will now investigate the average domain size as a function of time. We define the mean linear size of a phase domain by Here, v + n denotes the averaged number of grid points in a connected region, where φ is larger than φ eq , and v n is the total number of grid points. Thus, the domain size lies within the range 1/v n ≤ L ≤ 1, and L(t) should increase during the coarsening dynamics towards the steady swimming state. The fluctuations δφ of the initial profile are normal distributed with zero mean such that at t = 0 half of the grid points have φ > φ eq . They cannot all be isolated but rather belong to small connected regions with L ≈ 5/v n , where we extracted the factor √ 5 from our simulations at t = 0. Furthermore, we expect the maximum length to be around L ≈ 1/2. So, in our simulations L(t) lies in the interval 5/v n ≤ L ≤ 1/2. Figure 5 shows L(t) averaged over 200 simulation runs for different noise strengths ξ. The other parameters are the same as before. We clearly see a separation of time scales of the coarsening dynamics for both cases, with and without noise. At early times, we find in both cases a power law behavior L(t) ∝ t 0.1 . Without noise, coarsening quickly speeds up at a rate L(t) ∝ t 1/2 and then slows down again to L(t) ∝ t 0.1 . In contrast, thermal fluctuations in the order parameter profile hinder early coarsening and the mean domain size continues to grow slowly with L(t) ∝ t 0.1 over several decades and then crosses over to a fast final coarsening with rate L(t) ∝ t 0.8 . The crossover time is only determined by the diffusion time τ D and does not depend on noise strength ξ. Interestingly, a similar observation to the second case has been made for coarsening in the dynamical model H, where the Cahn-Hilliard equation couples to fluid flow at low-Reynolds number via an advection term. A slow coarsening rate L(t) ∝ t 1/3 in a diffusive regime at short times is followed by an advection driven regime with L(t) ∝ t at later times [62,78,79,80]. Although we cannot simply reformulate our model as an advective Cahn-Hilliard equation, since the phase separation in our case is driven by the interfacial flow u| R itself, we observe similar coarsening regimes as in model H, when we include some noise.

Dynamics of the swimming state
We now consider the time regime t > 9, where the droplet moves in its steady swimming state. However, as can be observed in fig. 2 (a), the droplet speed v D (t > 9) in the swimming state strongly fluctuates since we have added a thermal noise term to the diffusion-advection-reaction equation (11) for the order parameter field φ. These fluctuations also randomly change the swimming direction e as the inset of fig. 6 illustrates, where we show an exemplary swimming trajectory r(t) = r(0) . Therefore, we expect the droplet to perform active Brownian motion or a persistent random walk. In a droplet with axisymmetric profile the swimming direction is perpendicular to the domain wall separating both phases. When the order-parameter profile fluctuates, we also expect the domain wall to fluctuate and thereby the swimming direction e. There are no other reasons to change the orientation of e. In ref. [43] we showed that a spherical and isotropic emulsion droplet, with Marangoni flow at its surface, does not experience a frictional torque, which could also change the swimming direction. Thus, for an arbitrary surface tension profile σ(θ, ϕ) a spinning motion of the droplet does not occur. But this also means that fluctuating flow fields in the surrounding fluid, which have to fulfill boundary condition (5) at the droplet surface, cannot generate a stochastic torque acting on the droplet. Therefore, in contrast to a rigid colloid, spherical emulsion droplets do not exhibit conventional thermal rotational diffusion.

Active Brownian motion of the droplet
To characterize the active Brownian motion of the droplet, we first discuss the mean squared displacement (MSD) ∆r 2 = [r(t) − r(0)] 2 , where we average over an ensemble of trajectories. Here, the droplet is already in the swimming state at t = 0, thus the MSD does not include the droplet's acceleration towards the steady swimming state as discussed in sect. 4. Figure 6 shows the MSD for a droplet with noise strength ξ = 5 · 10 −3 . At early times, the droplet moves ballistically since the MSD grows as ∆r 2 ∝ t 2 , while between t = 10 and t = 100 it crosses The trajectory is reminiscent of an active particle with constant speed and rotationally diffusing orientation vector e(t).
over to diffusive motion with ∆r 2 ∝ t. This motion persists as t → ∞. As expected, in the absence of noise, ξ = 0, we always observe ballistic motion ∆r 2 ∝ t 2 (not shown). The MSD for ξ = 10 −3 in fig. 6 does not cross over to diffusion in the plotted time range. In the following, we will discuss the influence of the noise strength ξ on the Brownian motion in more detail but we will first introduce what has become the standard model of an active Brownian particle [15,81,11,12]. If we assume the droplet speed v D and orientation vector e to be independent random variables, we can factorize the MSD as For active Brownian particles without any aligning field the swimming direction diffuses freely on the unit sphere, which one describes by the rotational diffusion equation ∂ t p(e, t) = D r ∇ s 2 p(e, t). Thus the orientational correlation function decays as [82,83] e(0) · e(t) = e −t/τr .
Here, the rotational correlation time τ r = 1/(2D r ) is the characteristic time it takes the droplet to "forget" about the initial orientation e(0). Hence, for times t < τ r the droplet swims roughly in the direction of e(0), while at later times t > τ r the orientation becomes randomized. Under the assumption of a constant swimming speed, i.e. v D (t ′ )v D (t ′′ ) = (v D ) 2 , one finds for the MSD Expression (20)   diffusive motion with for t ≫ τ r . Here, D eff is the effective translational diffusion constant. It neglects any contribution from thermal translational motion, which is o.k. for sufficiently large v D . Indeed, for the active droplet the rotational correlation function e(0) · e(t) decays exponentially as demonstrated in fig. 7 for different noise strengths and by fits to eq. (19). The rotational correlation time τ r , which acts as fitting parameter, is shown in the inset for various values of noise strength ξ. For ξ = 5 · 10 −3 , we find τ r ≈ 30, which is in agreement with the cross-over region from ballistic to diffusive motion in the MSD curve of fig. 6. Furthermore, from the asymptotic behavior at t ≪ τ r and t ≫ τ r of the MSD in fig. 6, we find v D ≈ 0.3 and D eff ≈ 1, respectively. This gives the rotational correlation time τ r = 3D eff /(v D ) 2 ≈ 33, which is close to the value determined from the orientational correlations. Thus D eff and τ r comprise the same information about the droplet trajectory r(t). However, the measurement of τ r in experiments or simulations can be done on much shorter time scales than D eff . Note that the relative fluctuations of the swimming speed about its mean value are small, as fig.  2(a) demonstrates. Therefore, they do not have a strong effect on D eff and we can safely use the mean value v D in eq. (21).
We do not know published experimental data for trajectories of active droplets in an unbounded fluid. However, fig. 1 of ref. [1] shows a trajectory of an active droplet confined between two glass plates. One can estimate the rotational correlation time τ r to be on the order of 100s. To compare this value with our model, we recapitulate the noise strength ξ = ℓ/R, which connects surfactant head size ℓ with droplet radius R, see sect. 2.3. If we assume, ξ ≈ 10 −4 . . . 10 −3 , we find from fig. 7 a rotational correlation time τ r ≈ 10 4 given in units of diffusion time Typical values for D are on the order of 10 −5 cm 2 /s [84]. Thus, for a droplet with R on the order of 10µm, one finds τ D ≈ 0.1s and the rotational correlation time τ r ≈ 10 3 s. This is only a factor 10 larger than the estimated value of 100s from ref. [1]. Given some uncertainties in our estimate such a difference can be expected. Nevertheless, two causes for the discrepancy are thinkable. First and foremost, our model droplet is allowed to move freely in the bulk fluid, while the real droplet of ref. [1] is confined between two plates, which limits the degrees of freedom and thus alters τ r . Secondly, active emulsion droplets are usually immersed in a surfactant laden fluid well above the critical micelle concentration. Hence, the surfactants from the bulk adsorb in form of micelles. This leads to local disturbances in the surfactant mixture at the front of the swimming droplet, and hence to an additional randomization of the droplet trajectory. We recently modeled the adsorption of micelles explicitly in a different system [43].

How fluctuations randomize the droplet direction
Now, we develop a theory how the noise strength ξ influences the rotational diffusion of the droplet direction. By increasing ξ in the diffusion-advection-reaction equation (11), the order parameter profile φ is subject to stronger fluctuations. In particular, these fluctuations affect shape and orientation of the domain wall separating the two regions with φ < φ eq and φ > φ eq from each other. The surface flow field is largest in this domain wall and thereby the orientation of the wall on the droplet interface determines the droplet swimming vector e. Thus, increasing noise strength ξ results in stronger fluctuations of e and ultimately a more pronounced rotational diffusion. The inset of fig. 7 confirms this scenario for the rotational correlation time τ r . Interestingly, for noise strengths up to ξ ≈ 3 · 10 −3 , one fits the data quite well by τ r ∝ 1/ξ 2 . Since the noise strength ξ was defined as ξ = ℓ/R in sect. 2.3, the rotational diffusion constant D r = 1/(2τ r ) of the active droplet behaves as D r ∝ 1/R 2 , which is in contrast to the scaling D r ∝ 1/R 3 of a passive colloid. In addition, one finds for the total number of surfactants N = 4πR 2 /ℓ 2 that D r ∝ 1/N . This can be understood from a simple hand-waving argument.
Fluctuations in the order parameter profile described in eq. (11) correspond to exchanging surfactant molecules (brominated against non-brominated and vice versa). A single event initiates an angular displacement ∆ϕ ≈ 2π/ √ N of the droplet direction e. These fluctuations take place on the diffusive time scale τ . Furthermore, from eq. (19) one finds a diffusive mean squared angular displacement (∆ϕ) 2 = 4D r t for times t ≪ τ r . Thus, for t on the time scale τ , one finds D r ∝ 1/N . In what follows, we want to explain the scaling τ r ∝ 1/ξ 2 more rigorously by applying perturbation theory to the thermal fluctuations of the order parameter profile around its steady profile.
As mentioned before, small fluctuations of the domain wall result in random changes of the droplet direction. We now apply perturbation theory to the fluctuating order parameter profile, which determines the surface tension profile and thereby the swimming direction according to eq. (8). We consider a droplet, which initially swims in z-direction and changes its direction in x and/or y-direction, hence e = e z + δe. We write down a perturbation ansatz for the surface tension profile, σ = σ 0 +δσ, with the unperturbed axisymmetric part σ 0 = ∞ l=1 s 0 l Y 0 l and the perturbation δσ = s 1 1 Y 1 1 + s −1 1 Y −1 1 , where we only include the coefficients s ±1 1 , which are responsible for changes δe, as one recognizes from eq. (8). By linearizing the equation of state (9) around φ eq , one can connect the coefficients s m l of σ directly to the expansion coefficients of the order parameter field φ. Writing φ = φ 0 +δφ, where φ 0 describes the unperturbed steady-state field and δφ its fluctuations, we find φ 0 = aσ 0 and δφ = aδσ, where the factor a is given in appendix E. Similarly, one decomposes j D and u| R into their steady-state fields and a fluctuating small perturbation (see appendix E). This allows us to derive from the field equation (11) of the order parameter, the dynamic equation linear in the fluctuating perturbations: From our study of the coarsening dynamics we know that the first and second term on the right-hand side describe a relaxation towards steady state on times t < 10. The rotational diffusion of the droplet direction occurs on time scales much larger and can only be due to the noise term. Extracting from Eq. (22) the coefficients s ±1 1 relevant for δe, we obtain A more thorough derivation of Eq. (23) is presented in appendix F. We have decomposed noise ζ into its multipole moments, ζ = l,m ζ m l Y m l . Projecting the variance of eq. (12) onto the relevant spherical harmonics, we obtain the fluctuation-dissipation theorem Assuming a constant speed v D during the reorientation of the droplet, we use eq. (23) in eq. (13b) for the droplet velocity vector to formulate the stochastic equation for rotations of the direction vector e: where we introduced the rotational noise vector By comparing eq. (25) with the Langevin equation for the Brownian motion of a particle's orientation e due to rotational noise η r : ∂ t e = √ 2D r η r × e [85], we identify δζ = η r × e and ξ √ 6πv D (2 + 3ν)a = 2D r .
Hence, the rotational correlation time τ r = 1/(2D r ) scales as τ r ∝ 1/ξ 2 with noise strength ξ. This confirms the fit in the inset of fig. 7 for noise strengths up to ξ ≈ 10 −3 . For larger ξ, the fluctuations start to very strongly disturb the domain wall. The illustration of fig. 8 is no longer valid and with it the the perturbation theory breaks down. Instead, the droplet loses its persistent swimming axis and the motion becomes purely erratic, which manifests itself in a rapidly decreasing τ r . Thus, beyond the time scale, the order parameter profile needs to reach its steady state, and for ξ < 10 −3 , the dynamics of the swimming active emulsion droplet is equivalent to the dynamics of an active Brownian particle with constant swimming velocity and rotationally diffusing orientation vector e.

Conclusions
In this paper we considered an active emulsion droplet, which is driven by solutocapillary Marangoni flow at its interface [1]. A diffusion-advection-reaction equation for the surfactant mixture at the droplet interface, which we formulated in ref. [38], is used together with the analytic solution of the Stokes equation [43]. By omitting the axisymmetric constraint and including thermal noise into the description of the surfactant mixture, we generalized the model of ref. [38] to a full three-dimensional system and thereby were able to focus on new aspects.
First, we explored the dynamics from a uniform, but slightly perturbed surfactant mixture to the uniaxial steady swimming state, where the two surfactant types are phaseseparated. In between the initial and the swimming state, the surfactant mixture is not axisymmetric, which we verified by introducing and evaluating a biaxiality measure. We then investigated in detail the coarsening dynamics towards the swimming state by means of the polar power spectrum of the surface tension σ as well as the average domain size of the surfactant mixture. The coarsening proceeds in two steps. An initially slow growth of domain size is followed by a nearly ballistic regime, which is reminiscent to coarsening in the dynamic model H [62].
Second, we studied the dynamics of the squirming droplet. Due to the included thermal noise, the surfactant composition fluctuates and thereby the droplet constantly changes its swimming direction performing a persistent random walk. Thus, the swimming dynamics of the squirming droplet is a typical example of an active Brownian particle. The persistence of the droplet trajectory depends on the noise strength ξ. It is characterized by the rotational correlation time, for which we find the scaling law τ r ∝ ξ −2 . In fact, we are able to explain this scaling by applying perturbation theory to the diffusion-advectionreaction equation for the mixture order parameter. Thus we can link the dynamics of the surfactants at the molecular level to the dynamics of the droplet as a whole.
We hope that our work initiates further research in the field of active emulsion droplets. A deeper theoretical understanding of the coarsening due to the Marangoni effect could help to understand the power laws that we found in our simulations. Furthermore, various extensions of this work are possible, e.g., the explicit implementation of micellar adsorption as discussed in ref. [43] or taking into account confining plates below and above the droplet via no-slip boundary conditions. Finally, a numerical study of the collective motion of active droplets, which swarm in experiments [1], is still missing in the literature but has been implemented for pure squirmers [28].
Exploring and understanding the swimming mechanisms of both biological and artificial microswimmers is one of the challenges in the field. Here, we demonstrated that this task involves new and fascinating physics. Having gained deeper insights into these mechanisms can help to further improve the design of artificial microswimmers and tailor them for specific needs such as cargo transport.
In eq. (29b) we used the divergence theorem and in eq. (29c) we converted the line integrals into sums over the element boundaries. Furthermore, we discretized δ(r i −r j ) by partitioning the surface into rhombi of area A ♦ (see fig.  1) and defined δ q,p = 1/A ♦ for q = p , 0 for q = p , where q and p are the indices of the respective boundaries of elements i and j. Three cases have to be considered. First, if the elements i and j are neither identical nor neighbors, δ q,p vanishes in eq. (29c) for all q and p. Second, for i = j, δ q,p = 1/A ♦ and n iq · n jp = 1 for all q and p. Finally, for neighboring elements there is one common boundary, where δ q,p = 1/A ♦ and n iq · n jp = −1. Thus, one finds: where N is the number of element boundaries. Here, Q ij = 1 if elements i and j are neighbors and zero otherwise. Note that in eq. (30), we assumed the same edge length l and number of boundaries N for all elements. This is reasonable for a refined icosahedron with 642 FVM elements, as discussed in sect. 3. The form of eq. (30) acknowledges the conservation law for the noise [86]. However, in simulations we did not observe any effect of the next-neighbor correlations and therefore simplified the noise to the expression (16b) in the main text. Furthermore, we take N = 6 and A ♦ = 3/4l 2 , since our grid is mostly hexagonal, which explains the prefactor 2N l 2 /A ♦ = 2 · 12 1/4 in eq. (15), when we redefine the noise vector by the following replacement, z → 2 · 12 1/4 z.

D Average over droplet interface
The average is taken over the azimuthal angle ϕ in the coordinate frame whose z-axis is directed along the swimming direction e.
Here, the front of the moving droplet is at θ = 0.

E Perturbation ansatz
The zero and first-order contributions of φ = φ 0 + δφ, j D = j D,0 + δj D , and u| R = u 0 + δu are given by: with parameters Here we used the values of sect. 2.3 for b 1 , b 2 , b 12 , φ eq and ν.
F Dynamic equation for s ±1

1
To derive a dynamic equation for the expansion coefficients s ±1 1 , we project the dynamic equation (22) for the perturbation δφ onto the spherical harmonics Y ±1 1 [see also eq. (31b)]. Employing the orthonormality relation of the spherical harmonics and using eqs. (27), we ultimately obtain with noise components ζ m l defined in eq. (24). Due to the nonlinear advection term M φu| R in eq. (11), the coefficients s ±1 1 couple to s 0 2 . The term in square brackets on the right-hand side describes a relaxational dynamics for s ±1 1 . In particular, for the parameters chosen we find the swimming droplet to be a pusher. Thus, according to eq. (10) the coefficient s 0 2 > 0 and the term in square brackets is always negative. On time scales larger than the relaxation time, we can ignore the relaxational dynamics and the time dependence of the order parameter perturbation is solely determined by the thermal noise term, which confirms relation (23).
Note that in the dynamic equation for s 0 l equivalent to eq. (33), the advective term ∝ M is always positive and triggers for l = 1 and for sufficiently large M the onset of forward propulsion of the droplet (see fig. 3 and ref. [38]).