Non-reciprocal interactions between condensates in chemically active mixtures

Abstract

We study the behaviour of catalytically active droplets in multi-component conserved mixtures affected by noise. Working in the thin interface limit, we analytically determine the state diagram of the system, characterized by multiple dynamical regimes, and verify our findings using numerical simulations. In particular, we show the emergence of a non-reciprocal, chemically-mediated interaction between the droplets, which leads to the formation of (meta-)stable clusters of droplets of different species. We find that the clusters can display self-propulsion in a large part of the parameter space, including regions where the non-reciprocal interactions between the droplets are purely attractive. This surprising feature arises from the non-local nature of the chemical interactions, and points to locality violations as a general mechanism for energy dissipation and emergence of out-of-equilibrium steady states in active matter.

The distinguishing feature of an active system is the continuous energy consumption of its microscopic constituents, which drives the system out of equilibrium leading to complex and novel interactions. In multiple experimental systems—both biological and artificial—such interactions are due to the production and consumption of chemicals, which in turn affect the dynamics of the species constituting the system [1, 2, 3, 4, 5]. Examples include large collections of catalytically active colloids [6, 1, 7, 8], bacterial colonies and collections of cells [9, 10, 11, 12, 13], as well as metabolically active condensates [14, 15, 16, 17, 18, 19, 20], and suspensions of enzymes [21, 22, 23]. Chemical activity of enzymes manifests in a variety of non-equilibrium behaviour including growth and division of droplets [24, 25, 26, 27]. Such mixtures are naturally present in the cytoplasm of cells and it is believed that their activity is crucial in allowing them to perform their biological functions [14, 28].

When multiple chemically active species are present in a mixture, the effective interaction mediated by the chemical often break action-reaction symmetry [29, 30, 31, 32, 33, 34]. This has been studied in droplets [35], mixed bacterial colonies [36, 37, 38] and catalytic colloids [39, 40, 41, 42], showing the formation of self-propelling clusters. The formation of multi-component, multi-phase mixtures has been observed also in vivo in embryo cells [43]. The physical mechanism via which the chemical affects the species is not fully understood and likely system-specific, especially for phase separating biopolymers [27]. It has been argued that the enzymes in these mixtures display an effective chemotactic behaviour [44], leading to a mass current of the species proportional to the chemical gradients [19, 45, 46, 47]. Moreover, it has been shown in multiple experiments on phase separating droplets [15, 16, 17, 48, 49, 50] the chemical affects the species by altering the the inter-species interactions determining the phase separation, effectively changing the critical point of the mixture. These observations cannot be fully accounted for by chemotactic currents (see however [45]), and rather suggest that the free energy of the mixtures depend on the local chemical concentration, as illustrated in Figs. 1(a) and 1(b).

Figure 1: (a) The species gather into droplets due to their interaction. Spatial variations in the densities of the species ($\phi_{1}$ and $\phi_{2}$) create a spatially varying chemical field ($c$) as they produce and consume the chemical. (b) The chemical affects the intra-species interaction and mediates effective inter-species interactions. (c) Configurations of droplets: (I) disjoint, (II) partially overlapping and (III) fully overlapping. (d) Attractive interactions give rise to chasing-behaviour. The centre of the droplet is an emitter (E) of a chemical signal and the interfaces are sensors (S). An example where purely attractive interactions lead to chasing (left) and its implementation with two droplets (right).

In this Letter, we explore the dynamics of a multi-component mixture in the presence of a chemical which affects its material properties and is in turn produced and consumed by the components of the mixture. We find the emergence of (meta-)stable clusters of droplets of different species, displaying different mechanisms of self-propulsion. The first mechanism is chasing of droplets due to non-reciprocal, attractive-repulsive interactions, similarly to what is observed in [29, 40, 36, 38]. The second mechanism is non-locality-induced chasing of mutually attractive droplets, defying the expectation that self-propulsion requires attractive-repulsive interactions. Finally, the third mechanism is due to the interplay between non-reciprocal couplings and fluctuations in deterministically static configurations.

Theoretical framework—In order to describe the dynamics of the system, we start by modelling the evolution of the density fields for $N$ phase-separating species in $d$ dimensions, $\phi_{a}(\bm{r},t)$, $a\in\{1,2,\cdots,N\}$, which are governed by a free energy $\mathcal{F}[\{\phi_{a}\},c]$ that depends on a chemical field $c$ that the two species produce or consume, e.g. through dependence of the material properties on the chemical. We thus consider a generalization of the model B dynamics [51], as follows

$$\partial_{t}\phi_{a}=\bm{\nabla}\cdot\left[M_{a}\bm{\nabla}\dfrac{\delta\mathcal{F}[\{\phi_{b}\},c]}{\delta\phi_{a}}\right]+\sqrt{2M_{a}T}\bm{\nabla}\cdot\bm{\Lambda}_{a},$$ (1a)

$$-D_{c}{\bm{\nabla}}^{2}c+\gamma c=\sum_{a}r_{a}\phi_{a},$$ (1b)

where $M_{a}$’s are the corresponding mobilities, $\bm{\Lambda}_{a}$’s are Gaussian stochastic currents with zero mean and correlations $\langle\Lambda_{a,i}(\bm{r},t)\Lambda_{b,j}(\bm{r}^{\prime},t^{\prime})\rangle=\delta_{ij}\delta_{ab}\delta^{d}(\bm{r}-\bm{r}^{\prime})\delta(t-t^{\prime})$ [52], and $T$ is the temperature. In Eq. (1b), $r_{a}$’s are the production (consumption) rates of the chemical by the two species, $D_{c}$ is the diffusion coefficient of the chemical, and $\gamma\equiv D_{c}\kappa^{2}$ is its decay rate, which introduces a decay length scale $\kappa^{-1}$. The dynamics of $c$ is assumed to be significantly faster than the dynamics of the two species, leading to a separation of time scales in Eq. (1b).

We consider the case where the free energy of the system can be decomposed into independent parts, namely, $\mathcal{F}[\{\phi_{b}\},c]=\sum_{a}\mathcal{F}_{a}[\phi_{a},c]$. While this choice helps to simplify the presentation, we show that our results are not fundamentally modified by the presence of (weak) inter-species interactions (see Appendix A). The free energies are defined as $\mathcal{F}_{a}=\int d^{d}\bm{r}\left[\frac{1}{2}K_{a}(\phi_{a},c)({\bm{\nabla}}\phi_{a})^{2}+U_{a}(\phi_{a},c)\right]$ where the $U_{a}$ terms are double-well potentials with $\mathbb{Z}_{2}$ symmetry (see Appendix B), which indicate the possibility of phase separation in both species, associated with the formation of droplets of the corresponding dense phase in coexistence with the corresponding dilute background. The characteristic length-scales that control the interface thicknesses for the two species are significantly smaller than other relevant length scales when the system is far from criticality, and can thus be set to zero [53].

In the sharp interface limit, the densities in and out of the droplets are given by the two minima of the field potentials $U_{a}(\phi_{a})$, which we denote by $\bar{\phi}_{a,+}$ for the dense phases and $\bar{\phi}_{a,-}$ for the dilute phases, leading to density contrasts of $\Delta\bar{\phi}_{a}=\bar{\phi}_{a,+}-\bar{\phi}_{a,-}$ across the droplet interfaces. Since we are interested in the linear response regime, we can treat these densities as independent of $c$. For the thin interface limit, the surface tensions $\sigma_{a}$ of the interfaces of the droplets can be calculated in terms of $U_{a}$ and $K_{a}$, which depend on $c$ [53, 54]. For relatively small values of $r_{a}$, we can consider a linear response regime where $\sigma_{a}(c)=\sigma_{a}(0)+\frac{\partial\sigma_{a}}{\partial c}c$, where the derivative is evaluated at vanishing concentration (i.e. $\left.\frac{\partial\sigma_{a}}{\partial c}\right|_{c=0}$).

We assume that the long-time behaviour of the system is characterized by the presence of spherical droplets of radii $R_{a}$ (one for each species), with volume $\pazocal{V}_{a}=V_{d}R_{a}^{d}$ and surface area $\pazocal{S}_{a}=dV_{d}R_{a}^{d-1}=d\pazocal{V}_{a}/R_{a}$, where $V_{d}=\pi^{d/2}/\Gamma(\frac{d}{2}+1)$ is the volume of unit sphere in $d$ dimensions. We can use the dynamical equation for the concentration fields (1b) to extract effective stochastic equations of motion for the centre of mass of the droplet of species $a$, which we denote by $\bm{x}_{a}$ (see Fig. 1), as follows [54]

$$\dot{\bm{x}}_{a}=-\mathcal{M}_{a}{\cal S}_{a}\frac{\partial\sigma_{a}}{\partial c}{\bm{\nabla}}{\bar{c}}_{{\rm s},a}(\bm{x}_{a})+\sqrt{2\mathcal{M}_{a}T}\,\bm{\eta}_{a}(t),$$ (2)

where ${\bar{c}}_{{\rm s},a}(\bm{x}_{a})\equiv\langle c\rangle_{{S}_{a}}=\frac{1}{d{V}_{d}}\int d\hat{\bm{n}}_{a}c(\bm{x}_{a}+R_{a}\hat{\bm{n}}_{a})$ with $\hat{\bm{n}}_{a}$ representing the radial unit vector describing the points on the boundary $S_{a}$ of droplet $a$, $\mathcal{M}_{a}=M_{a}d/(\pazocal{V}_{a}\Delta\bar{\phi}_{a}^{2})$ is an effective mobility coefficient, and $\bm{\eta}_{a}$’s are Gaussian white noise terms with zero mean and correlations $\langle\eta_{a,i}(t)\eta_{b,j}(t^{\prime})\rangle=\delta_{ij}\delta_{ab}\delta(t-t^{\prime})$. We note that Eq. (2) is independent of the details of the free energy, and can be regarded as a generalized (diffusio-)phoretic transport mechanism resulting from interfacial stresses [55]. It is interesting to note that Eq. (2) describes the dynamics of droplets even when they overlap. Within the thin interface approximation, it is possible to derive exact expressions for ${\bar{c}}_{{\rm s},a}(\bm{x}_{a})$ in different dimensions (see Appendix C for details).

When the droplets do not overlap, we can obtain a closed-form expression for the surface-average concentration and use it to express the drift velocity as follows (see Appendix C for the derivation)

$$\dot{\bm{x}}_{a}=-\frac{d^{2}M_{a}}{R_{a}\Delta\bar{\phi}_{a}^{2}}\frac{\partial\sigma_{a}}{\partial c}\,{\cal I}_{d}\left(\kappa R_{a}\right){\bm{\nabla}}c(\bm{x}_{a}).$$ (3)

where ${\cal I}_{d}(b)=2^{\frac{d}{2}-1}\Gamma\left(\frac{d}{2}\right)b^{-\frac{d}{2}}\left[dI_{\frac{d}{2}}({b})+{b}I_{\frac{d}{2}+1}({b})\right]$, with $I_{\nu}$ being the modified Bessel function of the first kind. Note that the effective phoretic mobility arising from this mechanism shown in Eq. (3) exhibits a leading order $1/R$ size dependence in any dimension $d$, which is in sharp contrast with the case of standard diffusiophoresis where the mobility is independent of size [1]. Moreover, in the derivation of this result no thin slip-layer assumption is used [55, 44], and near-field effects are taken into consideration [56, 57], as represented by the function ${\cal I}_{d}\left(\kappa R_{a}\right)$, which introduces additional $d$-dependent size dependence.

Our theoretical framework leads to the emergence of non-reciprocal interactions between different droplets. To demonstrate this effect and the phenomenology that emerges from it, we focus on the two-body dynamics of droplets comprising two different species. We define the relative position of the two droplets as $\bm{x}=\bm{x}_{2}-\bm{x}_{1}$ and their midpoint as $\bm{X}=\frac{1}{2}(\bm{x}_{1}+\bm{x}_{2})$. The dynamics of these coordinates will take the general form

$$\dot{\bm{x}}=-\mathcal{M}\bm{\nabla}_{\bm{x}}W_{-}(\bm{x})+\sqrt{2\mathcal{M}T}\,\bm{\eta}_{x}(t),$$ (4a)

$$\dot{\bm{X}}=\bm{v}(\bm{x})+\sqrt{\mathcal{M}T/2}\,\bm{\eta}_{X}(t),$$ (4b)

where $\bm{\eta}_{x}(t)$ and $\bm{\eta}_{X}(t)$ are two independent Gaussian white noise terms with zero mean and unit variance, $\mathcal{M}=\mathcal{M}_{1}+\mathcal{M}_{2}$, and $\bm{v}(\bm{x})=-\frac{1}{2}\mathcal{M}\bm{\nabla}_{\bm{x}}W_{+}(\bm{x})$, with the potentials $W_{\pm}(\bm{x})=\frac{\mathcal{M}_{2}}{\mathcal{M}}{\cal S}_{2}\frac{\partial\sigma_{2}}{\partial c}{\bar{c}}_{{\rm s},2}(\bm{x})\pm\frac{\mathcal{M}_{1}}{\mathcal{M}}{\cal S}_{1}\frac{\partial\sigma_{1}}{\partial c}{\bar{c}}_{{\rm s},1}(\bm{x})$ depending only on $\bm{x}$ due to translational invariance of the system.

To demonstrate the rich behaviour of the system, we start with $d=1$, where a droplet of species $a$ is characterized by the positions $x_{a,l}=x_{a}-R_{a}$ and $x_{a,r}=x_{a}+R_{a}$ of its left and right interfaces, respectively. Consequently, we can obtain a simple expression for the surface-averaged concentration as ${\bar{c}}_{{\rm s},a}(x_{a})=[c(x_{a}-R_{a})+c(x_{a}+R_{a})]/2$. Note, however, that this choice is made for simplicity of the presentation, as analytical results can be obtain in all dimensions, as shown in Appendix C. From here onward, We set $\bar{\phi}_{1,+}=\bar{\phi}_{2,+}=1$, $\bar{\phi}_{1,-}=\bar{\phi}_{2,-}=0$, $D_{c}=1$, and $\gamma=1$, and hence express all length scales in terms of the decay length of the chemical $\kappa^{-1}=\sqrt{{D_{c}}/{\gamma}}=1$ [58]. Moreover, we assume $R_{2}>R_{1}$ without loss of generality.

Most of our results hold for any profile of the chemical concentration. To close the equations, we need to solve for the concentration field using Eq. (1b), which yields

$$\partial_{x}c(x)=\frac{1}{2}\sum_{a}r_{a}\left(e^{-|x-x_{a}+R_{a}|}-e^{-|x-x_{a}-R_{a}|}\right),$$ (5)

in 1D. Using Eq. (5), we obtain explicit expressions in the form of Eq. (4b), with the potentials $W_{\pm}(x)$ given in Appendix D. We observe that the expressions for $W_{\pm}(x)$ depend on the coupling constants only through the combinations

$$\mathcal{C}_{a\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}b}\equiv-r_{b}M_{a}\frac{\partial\sigma_{a}}{\partial c},$$ (6)

for $a\neq b$. This is because a droplet $a$ climbs (or descends) chemical gradients at a rate proportional to $-M_{a}\frac{\partial\sigma_{a}}{\partial c}$, while due to Eq. (5) it senses a chemical gradients proportional to the production rate $r_{b}$ of the other droplet. The $\mathcal{C}_{a\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}b}$ coupling constant then specifies the strength of the interaction of experienced by $a$ toward $b$. We will refer to them as chasing rates. Note that for positive $\mathcal{C}_{a\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}b}$ species $a$ is attracted to species $b$, and, in general, $\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}\neq\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}$, as expected in active systems [29, 40].

Deterministic dynamics—We can now study the behaviour of the system as a function of chasing rates and droplet sizes. In Fig. 2(a), we present a classification of the different states of the dynamics by using the dimensionless parameters $\mathcal{R}\equiv R_{1}/R_{2}$ as radius and $\vartheta\equiv\arg(\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}+\mathrm{i}\,\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1})$ as angle (see Appendix B for the details of the rules used in the classification).

Figure 2: (a) State diagram obtained by simulating Eq. (1b), presented a radial plot with radius $\mathcal{R}=R_{1}/R_{2}$ from 0 to 1 and angle $\vartheta=\arg(\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}+\mathrm{i}\,\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1})$. The black [Eq. (7)] and orange [Eq. (8)] lines delimit the boundaries between the different configurations predicted by theory. (b) Effective potentials for $W_{-}(\bm{x})$ in the various state. (c) Steady-state velocity and relative position from Eq. (9) and (10) (dashed black), and comparison with the result from simulations, where colour indicaters time, with the steady-state shown in orange. The simulations are obtained with $T=10^{-5}$.

Starting our analysis with no noise, we note that Eq. (4b) always admits $x=0$ as a stationary point, corresponding to the symmetric configuration with the small droplet encapsulated inside the large one, i.e. configuration (III) of Fig. 1(c). This fixed point is stable when

$$\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}<\frac{R_{2}}{R_{1}}\frac{1}{\tanh R_{1}}\,\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2},$$ (7)

which means that if Eq. (7) is satisfied, droplets starting will reach a symmetric configuration of complete engulfment when starting from asymmetric initial conditions [see Fig. 2(a)]. For sufficiently large initial distances [droplets in configuration (I) of Fig. 1(c)], the potential decreases with $|x|$ when

$$\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}>\frac{R_{2}}{R_{1}}\dfrac{\tanh R_{2}}{\tanh R_{1}}\,\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}.$$ (8)

In this case, two well-separated droplets experience an effective repulsion (and their distance increases indefinitely). Note that Eq. (7) and (8) can be satisfied at the same time. In this case, $W_{-}(x)$ has a metastable minimum at $x=0$, while disjoint droplets are repulsive [see Fig. 2(b)]. Finally, when both Eq. (7) and (8) are violated, $W_{-}(x)$ becomes a double-well potential with two minima [see Fig. 2(b)], representing configuration (II) of Fig. 1(c). In this case, the fixed point at $x=0$ is unstable, while two stable fixed points appear at $\pm x^{*}$, where

$$x^{*}=\dfrac{1}{2}\ln\left[1+\left(e^{2R_{2}}-e^{2R_{1}}\right)\dfrac{\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}R_{1}-\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}R_{2}}{\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}R_{1}+\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}R_{2}}\right].$$ (9)

This corresponds to the partially overlapping configuration (II) of Fig. 1(c), with the smaller droplet being on the left and right of the larger one, respectively. The steady-state velocity $v^{*}=v(x^{*})$ of the midpoint of the cluster is then found [see Eq. (4b)]

$$v^{*}=\frac{\sinh\left(R_{2}-R_{1}\right)\,\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}\,\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}}{(\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}R_{1}+\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}R_{2})}e^{-x^{*}}.$$ (10)

These predictions are verified by simulating Eq. (1b) with a small value of noise ($T=10^{-5}$) to check if the final state is in a disjoined, overlapping, or engulfed configuration [respectively, in configurations (I), (II), or (III) of Fig. 1(c)], and it is indeed observed that the shape of the effective potential $W_{-}$ determines the behaviour; see Figs. 2(a),  2(b), and 2(c).

It is remarkable that self-propelling clusters are observed in the region of fully attractive, albeit non-reciprocal, droplet interactions where $\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}>0$ and $\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}>0$. This may be surprising, as self-propulsion normally emerges when the chasing rates have opposite signs. However, since droplets are extended objects, additional opportunities to realize conditions for effective non-reciprocal interactions may arise from the fact that the sources and sinks are chemicals are in the bulks of the droplets while the force-sensing operates at their interfaces. The chasing rates determine how the interface of each droplet interacts with the centre of the other. This implies, for example, that configurations in which the leftmost interface of droplet $2$ lies between the leftmost interface and centre of mass of droplet $1$ can lead to a net repulsive effect on droplet $1$, despite all the forces being attractive from the point of view of the respective interfaces; see Fig. 1(d). Note that this effect is only possible in active systems.

As a consequence of the change in self-propulsion mechanism when $\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}$ changes sign (corresponding to $\vartheta=\frac{\pi}{2}$) the cluster inverts its velocity. While for $\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}<0$ droplet $1$ on the left is repelled by droplet $2$ that is chasing it, for $\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}>0$ droplet $1$ chases droplet $2$ due to the effective non-reciprocity described above. This is shown in Fig. 2(c), where we present the velocities and relative positions for a set of simulations, within, and just outside, the overlapping regime. In both cases, the solutions relax over time towards the analytical predictions.

Figure 3: (a) Histogram of the jumping times for the relative position in the overlapping regime at $T=2\times 10^{-2}$, following a Poisson distribution with measured jumping time $\langle\tau\rangle=1.29\times 10^{3}$, compared with $\tau_{K}=1.44\times 10^{3}$. Inset: Plot of centre positions show a clear distinction between running and tumbling. (b) MSD of the cluster from simulations of Eq. (1b) in the overlapping (red) and engulfed (blue) regimes. The latter follows an AOU process with coefficients obtained from Eq. (4b) (dashed orange). (c) Steady-state probability distribution for $x$. Points show simulations of overlapping (red) and attractive (blue) regimes, and compared with the predictions from Eq. (4b) (dashed lines).

Stochastic dynamics—We now turn our attention to the finite temperature case, focusing on the attractive and overlapping regimes. In the overlapping regime of Fig. 2, Eq. (4a) describes a Langevin dynamics in a double-well potential. At sufficiently low temperatures, the dynamics of the relative position is confined in the close neighbourhood of the two minima $\pm x^{*}$, while rare trajectories exhibit barrier-crossing, going from one minimum to the other. The mean first-passage time can be calculated from Kramers’ law as $\tau_{K}=\frac{2\pi}{\sqrt{|W_{-}^{\prime\prime}(0)W_{-}^{\prime\prime}(x^{*})|}}\exp\left\{\left[W_{-}(0)-W_{-}(x^{*})\right]/T\right\}$ [59, 60]. Each time the relative position passes from one minimum to the other, the midpoint $X$ of the cluster inverts its direction of motion. This results in an emergent run-and-tumble dynamics [61, 39], which can be modelled as $\dot{X}=v^{*}e(t)+\sqrt{\mathcal{M}T/2}\,\eta_{X}(t)$, where $e(t)=\pm 1$ is a random Poisson process with rate $1/\tau_{K}$ characterizing the tumbling motion [62]. The probability distributions for the stochastic jump times and relative positions obtained from our simulations of Eq. (1b) show good agreement with the prediction from Kramers theory [see Fig. 3(a)] and the Boltzmann distribution associated with the potential $W_{-}$ [see Fig. 3(c)]. The run-and-tumble dynamics of the midpoint is then characterized by a transition of the mean-squared-displacement (MSD) $\langle X(t)^{2}\rangle$ from ballistic to diffusive behaviour at $\tau_{K}$, which is is consistent with the behaviour observed from Eq. (1b), as shown in Fig. 3(b).

In the engulfed state at zero temperature, the steady-state is a static configuration with the smaller droplet symmetrically placed inside of the larger one, leaving no sign of the chemical activity in the droplet dynamics at large times. However, when temperature is increased, the relative position starts to fluctuate around $x=0$, which in turn affects the dynamics of the midpoint. At small temperatures we can expand Eq. (4b) around $x=0$, obtaining the equations of motion of an active Ornstein-Uhlenbeck process (AOU) [63]. We thus retrieve a signature of the chemical activity of the system in the MSD of the midpoint, and quantify it analytically (see [58]) and via numerical simulations of Eq. (1b) [see Figs. 3(b) and 3(c)]. We find that active fluctuations enhance the diffusivity of the engulfed clusters at large times, substantially broadening the range of parameters for which the droplet dimer exhibits self-propulsion. Such active fluctuations can provide the possibility to control the self-organization of these compound active particles [64].

Higher dimensions—We now briefly examine the role of dimensionality to probe the robustness of our findings. Due to rotational symmetry, $W_{-}$ is a radial function, which can be a attractive, repulsive, or in the form of (inverted) Mexican hat potential, depending on the radii and chasing rates. When the potential forms a Mexican hat, its degenerate minimum is a circle or sphere of radius $x^{*}$. When the potential is strong as compared to the noise, the dynamics is confined to this minimum, such that the relative position satisfies $\bm{x}(t)\simeq x^{*}\hat{\bm{e}}(t)$, where $\hat{\bm{e}}(t)$ is a unit vector, which evolves according to $\dot{\hat{\bm{e}}}(t)=\frac{1}{x^{*}}\sqrt{2\mathcal{M}T}\,\hat{\bm{e}}(t)\times[\hat{\bm{e}}(t)\times\bm{\eta}_{x}(t)]$. As before, the dynamics of the midpoint reads $\dot{\bm{X}}=|\bm{v}(x^{*})|\hat{\bm{e}}(t)$ (ignoring its own noise term). We thus find that the self-propelled cluster behaves effectively like an active Brownian particle [65, 7] with angular diffusion coefficient $T\mathcal{M}/(x^{*})^{2}$.

Figure 4: (a) State diagram at finite (large disk) and vanishing (small disk) decay rate $\gamma$ in $d=3$. (b) Results from simulations of (1b) in $d=2$ at $T=10^{-4}$. (b) Droplet trajectories (droplet 1 in blue, droplet 2 in black). Inset: Snapshot of the densities $\phi_{a}$, red and blue, and the chemical field $c$. (c) MSD of the midpoint for droplets in the attractive (blue) and overlapping (red) configurations. (d) The effective potential for $\bm{x}$, inferred from the probability distribution of the relative position, which is consistent with the theoretical analysis. Dashed lines are running averages to guide the eye.

We analytically calculate $W_{-}$ in $d=3$ and use it to build the state diagram (see Appendix C and Ref. [58]), which results qualitatively similar to the 1D case. At finite decay rate $\gamma$ [Fig. 4(a)], non-locality allows for the formation of self-propelling clusters of mutually attractive droplets, highlighting the general nature of this mechanism and of its relevance in the dynamics of multi-component mixtures. Working in 3D allows us to study Eq. (1b) in the case of $\gamma=0$, which is characterized by the lack of overlapping clusters at $\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2},\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}>0$ [Fig. 4(a)]. While the same procedure could be carried over also in $d=2$ (see Appendix C), we choose to take a complementary approach and present a numerical simulation of Eq. (1b) in this case. Again, in 2D we find overlapping, self-propelling clusters of droplets with both positive chasing rates [Fig. 4(b)], as well as attractive clusters with AOU behaviour. The MSD versus time plot [Fig. 4(c)] shows a transition from ballistic to diffusive motion at a time consistent with our analysis. We also measure the distribution of the separation of the droplets, which in 2D is related to the effective potential as $W_{-}(\bm{x})\propto-\ln\frac{1}{x}P(x)$, to estimate the effective potential in the engulfed and overlapping regimes [Fig. 4(d)].

Discussion—We have investigated the dynamics of catalytically active multi-component mixtures. Our study indicates that droplet motion within these mixtures is primarily governed by the properties of their surface tension and reaction kinetics when the interactions causing the phase separation depend on concentrations of chemicals that participate in the catalytic activity. We observe self-propelling clusters in a large part of the state-diagram, due to the emergence of multiple mechanisms for self-propulsion, which depend on the interplay between the non-local nature of the chemical interactions, the sizes of the droplets, and thermal fluctuations.

A prominent emerging feature in this system is that the effective coupling constants controlling the non-reciprocal interactions between two droplets depend inversely on the the droplet sizes, implying that two droplets can cross the threshold of mutually attractive to chasing interaction simply by ripening. The resulting strong and non-trivial dependence of the properties of self-propulsion on the sizes of the droplets could be of significant help in experimental validation of the aforementioned findings. The self-propulsion induced by chemical activity leads to enhanced diffusion of the droplet clusters, providing a method to probe the chemical activity of cellular sub-compartments in vivo. We therefore argue our results can help design and analyze chemically active matter, in addition to shedding light on the interplay between activity and non-locality in extended systems. In the future, it will interesting to extend this study to the case of chemically active suspensions with a large number of components [66], investigate the cases when metabolic fluxes can themselves be the drivers of non-equilibrium phase separation [45], and explore possible roles of the emerging phenomenology in origins of life scenarios [21].

End Matter

Appendix A The role of inter-species interaction

To test the robustness of our results with respect to the introduction of inter-species interactions, in particular, the existence of the purely attractive chasing states, we add a Florey-Huggins interaction term of the form $\chi\phi_{1}\phi_{2}$ to the free energy.

In Fig. 5, we illustrate the effect of $\chi$ on the state diagram, as obtained from simulation of Eq. (1b) in 1D, together with the self-propulsion velocity. We indeed observe the presence of chasing states at non-zero $\chi$ and positive $\mathcal{C}_{a\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}b}$. In particular, when $\chi<0$, the passive interaction between the droplets due to the cross term is also attractive. Nevertheless, simulations show the formation of self-propelling, overlapping states.

Figure 5: Simulations with non-zero $\chi$. To the left: state diagram categorizing the behaviour of the system as in Fig. 2. The green triangle indicates states where all interactions are attractive, but the system still self-propels. To the right: the corresponding steady-state velocity. Parameters: $R_{1}=1$, $R_{2}=1.5$ ($\mathcal{R}=\frac{2}{3}$), and $\partial_{c}\sigma_{2}(0)=0.8$.

Appendix B Details of the simulations

For numerical integration of the equations of motion, Eq. (1b), we have to make a specific choice for free energy. We consider

$$\mathcal{F}_{a} =\int d^{d}\bm{r}\left[\frac{1}{2}K(\bm{\nabla}\phi_{a})^{2}+\frac{1}{4}A_{a}(c)\left(\phi_{a}^{2}-\bar{\phi}_{a}^{2}\right)^{2}\right],$$ (11)

$$A_{a} =\left(\sqrt{A_{0}}+{\frac{3}{2\sqrt{2}}}\alpha_{a}c\right)^{2},$$ (12)

which yields $A_{a}(c=0)=A_{0}$ and $\partial_{c}\sigma_{a}(0)=\alpha_{a}$. This particular choice of the Ginzburg-Landau type potential and function $A_{a}(c)$ is for convenience, but different choices would lead to the same phenomenology (as long as the system is in the linear response of the dependence of the surface tension on the chemical concentration). For the simulations, we set $K=1$, $M_{a}=1$, $D_{c}=1$, $-r_{1}=r_{2}=1$, and $A_{0}=256$. In the simulations for the phase diagrams, velocity, and position, we include a finite but negligible temperature $T=10^{-5}$ for improved numerical stability. We furthermore set $c(\bm{q}=0)=0$, as this improves stability without affecting the results.

To classify the phases, we initialize two simulations in the configurations (I) and (III), and check if the final configurations are (I), (II), or (III). This gives us data-points such as $[(\mathrm{I}),(\mathrm{III})]\rightarrow[(\mathrm{I}),(\mathrm{I})]$, which indicate that both simulations evolve to the (I) configuration. This is thus classified as repulsive. The complete classification rules are as follows: $[(\mathrm{I}),(\mathrm{III})]\xrightarrow[]{\text{repulsive}}[(\mathrm{I}),(\mathrm{I})]$, $[(\mathrm{I}),(\mathrm{III})]\xrightarrow[]{\text{engulfed}}[(\mathrm{III}),(\mathrm{III})]$, $[(\mathrm{I}),(\mathrm{III})]\xrightarrow[]{\text{meta-stable}}[(\mathrm{I}),(\mathrm{III})]$, and $[(\mathrm{I}),(\mathrm{III})]\xrightarrow[]{\text{overlapping}}[(\mathrm{II}),(\mathrm{II})]$.

We simulate Eq. (1b) using pseudo-spectral methods with periodic boundary conditions, The system is stiff due to the fourth-order derivative and sharp gradients, which is dealt with by using exponential time-stepping [67]. For the simulations in $1D$, we use a length $L=16$ and discretize it in $N=512$ grid points, such that the thickness of the interface is $\xi=\sqrt{2/A_{0}}=\sqrt{8}\mathrm{d}x$. In $2D$, we set $N=64$ and $L=4$ for the simulation of the attractive cases, and $N=128$ and $L=8$ for the overlapping cases.

Appendix C Surface-averaged concentration in different dimensions

Here we present the details of the derivation of the expression for effective phoretic velocity of droplets in arbitrary dimensions. The deterministic part of Eq. (2) reads

$$\dot{\bm{x}}_{a}=-\frac{d^{2}M_{a}}{R_{a}\Delta\bar{\phi}_{a}^{2}}\frac{\partial\sigma_{a}}{\partial c}{\bm{\nabla}}{\bar{c}}_{{\rm s},a}(\bm{x}_{a}),$$ (13)

where

$${\bar{c}}_{{\rm s},a}(\bm{x}_{a})=\frac{1}{dV_{d}}\int d\hat{\bm{n}}_{a}c(\bm{x}_{a}+R_{a}\hat{\bm{n}}_{a}).$$ (14)

We can solve for the chemical concentration by using Eq. (1b) in Fourier space, which yields

$$c(\bm{x})=\sum_{b}\left(\frac{r_{b}}{D_{c}}\right)\int\frac{d^{d}\bm{q}}{(2\pi)^{d}}\,\frac{\phi_{b}(\bm{q})}{q^{2}+\kappa^{2}}\,e^{i\bm{q}\cdot(\bm{x}-\bm{x}_{b})},$$ (15)

where

$$\phi_{b}(\bm{q})=\int d^{d}\bm{y}\phi_{b}(\bm{y})e^{-i\bm{q}\cdot\bm{y}}=\int d\hat{\bm{n}}_{b}\int_{0}^{R_{b}}r^{d-1}dr\,e^{-i\bm{q}\cdot\hat{\bm{n}}_{b}r},$$ (16)

since $\phi_{b}(\bm{y})=1$ for $|\bm{y}|\leq R_{b}$ and $\phi_{b}(\bm{y})=0$ for $|\bm{y}|>R_{b}$ in the thin interface limit. Combining (15) and (16), we obtain

$${\bar{c}}_{{\rm s},a}(\bm{x}_{a})=\sum_{b}\left(\frac{r_{b}}{D_{c}}\right)\frac{1}{dV_{d}}\int d\hat{\bm{n}}_{a}\int d\hat{\bm{n}}_{b}\int_{0}^{R_{b}}r^{d-1}dr$$

$$\hskip 36.98866pt\times\int\frac{d^{d}\bm{q}}{(2\pi)^{d}}\,\frac{e^{i\bm{q}\cdot(\bm{x}_{a}-\bm{x}_{b}+\hat{\bm{n}}_{a}R_{a}-\hat{\bm{n}}_{b}r)}}{q^{2}+\kappa^{2}}.$$ (17)

By introducing the integral

$${\cal J}_{d}(p)=\frac{\int_{0}^{\pi}d\theta\,\sin^{d-2}\theta\,e^{ip\cos\theta}}{\int_{0}^{\pi}d\theta\,\sin^{d-2}\theta},$$ (18)

which can be evaluated as

$${\cal J}_{d}(p)=2^{\frac{d}{2}-1}\Gamma\left(\frac{d}{2}\right)p^{-\frac{d}{2}}\left[dJ_{\frac{d}{2}}(p)+pJ_{\frac{d}{2}+1}(p)\right],$$ (19)

in terms of the Bessel function $J_{\nu}$, we can write a closed form expression for the surface averaged concentration around droplet $a$ as follows

$${\bar{c}}_{{\rm s},a}(\bm{x}_{a})=\sum_{b\neq a}\left(\frac{r_{b}}{D_{c}}\right)\frac{(dV_{d})^{2}}{(2\pi)^{d}}\int_{0}^{\infty}\frac{dq}{q}\,\frac{{\cal J}_{d}(qR_{a})}{q^{2}+\kappa^{2}}$$

$$\hskip 36.98866pt\times\,{\cal J}_{d}(q|\bm{x}_{a}-\bm{x}_{b}|)\int_{0}^{qR_{b}}u^{d-1}du{\cal J}_{d}(u),$$ (20)

in any dimension $d$. The specific forms ${\cal J}_{1}(p)=\cos p$, ${\cal J}_{2}(p)=J_{0}(p)$, and ${\cal J}_{3}(p)=\frac{\sin p}{p}$ can be used to write explicit expressions for ${\bar{c}}_{{\rm s},a}(\bm{x}_{a})$ in different dimensions, namely

$${\bar{c}}_{{\rm s},a}^{(d=1)}(x_{a})=\frac{2}{\pi}\sum_{b\neq a}\left(\frac{r_{b}}{D_{c}}\right)\int_{0}^{\infty}\frac{dq}{q}\,\frac{\cos(qR_{a})\sin(qR_{b})}{q^{2}+\kappa^{2}}$$

$$\hskip 56.9055pt\times\,\cos(q|x_{a}-x_{b}|),$$ (21)

for $d=1$,

$${\bar{c}}_{{\rm s},a}^{(d=2)}(\bm{x}_{a})=\sum_{b\neq a}\left(\frac{r_{b}R_{b}}{D_{c}}\right)\int_{0}^{\infty}dq\,\frac{J_{0}(qR_{a})J_{1}(qR_{b})}{q^{2}+\kappa^{2}}$$

$$\hskip 56.9055pt\times\,J_{0}(q|\bm{x}_{a}-\bm{x}_{b}|),$$ (22)

for $d=2$, and

$${\bar{c}}_{{\rm s},a}^{(d=3)}(\bm{x}_{a})=\frac{2}{\pi}\frac{1}{|\bm{x}_{a}-\bm{x}_{b}|}\sum_{b\neq a}\left(\frac{r_{b}}{D_{c}R_{a}}\right)\int_{0}^{\infty}\frac{dq}{q^{3}}\,\frac{\sin(qR_{a})}{q^{2}+\kappa^{2}}$$

$$\qquad\times\,\left[\sin(qR_{b})-(qR_{b})\cos(qR_{b})\right]\sin(q|\bm{x}_{a}-\bm{x}_{b}|),$$ (23)

for $d=3$. These expressions are in agreement with those obtained using a more direct integration method.

In the special case when the droplets do not overlap, we can obtain a closed form expression in the form of generalized phoresis. To this end, we use a formal representation of the Taylor expansion as

$${\bar{c}}_{{\rm s},a}(\bm{x}_{a})=\frac{1}{dV_{d}}\smash{\int}d\hat{\bm{n}}_{a}e^{R_{a}\hat{\bm{n}}_{a}\cdot\bm{\nabla}}c(\bm{x}_{a}),$$ (24)

which can be written formally as follows

$${\bar{c}}_{{\rm s},a}(\bm{x}_{a})={\cal I}_{d}\left(R_{a}\nabla\right)c(\bm{x}_{a}),$$ (25)

where ${\cal I}_{d}(b)\equiv{\cal J}_{d}(p=-ib)$ [defined via Eq. (18)] is given as ${\cal I}_{d}(b)=2^{\frac{d}{2}-1}\Gamma\left(\frac{d}{2}\right)b^{-\frac{d}{2}}\left[dI_{\frac{d}{2}}(b)+bI_{\frac{d}{2}+1}(b)\right]$, in terms of the modified Bessel function of the first kind $I_{\nu}$. In particular, we have ${\cal I}_{1}(b)=\cosh b$, ${\cal I}_{2}(b)=I_{0}(b)$, and ${\cal I}_{3}(b)=\frac{\sinh b}{b}$. Since we are interested in contributions to $c$ from $b\neq a$ and because the droplets do not overlap, the right-hand (source) side of Eq. (1b) is zero in the region of interest, and we can write $\nabla^{2n}c=\kappa^{2n}c$ for $n\geq 1$. Combining this with Eqs. (13) and (25), we obtain

$$\dot{\bm{x}}_{a} =-\frac{d^{2}M_{a}}{R_{a}\Delta\bar{\phi}_{a}^{2}}\frac{\partial\sigma_{a}}{\partial c}\,{\cal I}_{d}\left(\kappa R_{a}\right){\bm{\nabla}}c(\bm{x}_{a}).$$ (26)

Note that the effective phoretic mobility arising from this mechanism shown in Eq. (3) exhibits a more complicated size dependence than the leading $1/R$ form, which depends on $d$.

Appendix D Effective potentials in 1D

We report here the expression for $W_{\pm}$ in Eq.(4b) for $d=1$. We consider the three subdomains corresponding to the three possible configurations, as illustrated in Fig. 1 (c): (I) disjoint, (II) partially overlapping, and (III) completely overlapping. These correspond to $|x|>R_{1}+R_{2}$, $R_{2}-R_{1}<|x|<R_{1}+R_{2}$, and $|x|<R_{2}-R_{1}$ respectively. The effective potentials $W_{\pm}(x)$ are given by:

$$W_{\pm}(x)=\frac{1}{\mathcal{M}}\begin{cases}\vskip 2.84544pt-2^{-\frac{(1\pm 1)}{2}}e^{-|x|}\left[\frac{\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}}{R_{2}}\cosh R_{2}\sinh R_{1}\mp\frac{\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}}{R_{1}}\cosh R_{1}\sinh R_{2}\right]+\text{const},&\text{(I)\phantom{II}}\\ \vskip 2.84544pt2^{-\frac{(3\pm 1)}{2}}e^{-(R_{1}+R_{2})}\left[\left(\frac{\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}}{R_{2}}\pm\frac{\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}}{R_{1}}\right)\cosh x+e^{-|x|}\left(\frac{\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}}{2R_{2}}\mp\frac{\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}}{2R_{1}}\right)(e^{2R_{2}}-e^{2R_{1}})\right],&\text{(II)\phantom{I}}\\ -2^{-\frac{(1\pm 1)}{2}}e^{-R_{2}}\left[\frac{\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}}{R_{2}}\sinh R_{1}\pm\frac{\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2}}{R_{1}}\cosh R_{1}\right]\cosh x+\text{const},&\text{(III)}\end{cases}$$ (27)

in different cases.

Non-reciprocal interactions between condensates in chemically active mixtures

Supplemental Material
Jacopo Romano^1,2, Martin Kjøllesdal Johnsrud¹, Benoît Mahault¹ and Ramin Golestanian^1,3

¹ Max Planck Institute for Dynamics and Self-Organization (MPI-DS), 37077 Göttingen, Germany.
² SISSA — International School for Advanced Studies and INFN, via Bonomea 265, 34136 Trieste, Italy
³ Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3PU, United Kingdom.

S..1 Rescaling of Eq. (1b) in the linear chemical response regime

The force acting on a droplet interface in response to chemical modulations depends on how the surface tension $\sigma$ and the binodal densities $\bar{\phi}_{\pm}$ vary with the chemical concentration $c(\bm{x})$. For weak coupling, one can linearize these dependencies as: $\sigma_{a}(c)=\sigma_{a}(0)+\partial_{c}\sigma_{a}(0)c$, $\phi_{a,\pm}(c)=\bar{\phi}_{a,\pm}(0)+\partial_{c}\bar{\phi}_{a,\pm}(0)c$, where the subscript $a$ denote the droplet species. The velocity $v_{a}$ of a droplet of species $a$ is then a function of the response coefficients: $v_{a}=v_{a}(\partial_{c}\sigma_{a}(0),\partial_{c}\bar{\phi}_{a,\pm}(0),c)$, with $v_{a}=0$ for $c={\rm const}$. In the weak-response regime, we retain only terms up to first order in the responses.

We argue now that the dependence of $v_{a}$ on $\partial_{c}\bar{\phi}_{a,\pm}(0)$ is negligible at first order. In fact, using the results of Ref. [54] (see in particular Sec. 4), we note that this dependence arises from two distinct contributions. The first is a correction of order $\mathcal{O}(\partial_{c}\bar{\phi}_{a,\pm}(0))$ to the droplet mobility. Since the force generated by chemical gradients is of order $\mathcal{O}(\partial_{c}\sigma_{a}(0))$, this correction affects the velocity only at the subdominant order $\mathcal{O}(\partial_{c}\bar{\phi}_{a,\pm}(0)\partial_{c}\sigma_{a}(0))$. The second contribution corresponds to an effective interfacial force that originates from mass displacements associated with changes in the binodal densities. This term scales as $\mathcal{O}(\partial_{c}\bar{\phi}_{a,\pm}(0)\,v_{b})$, where $v_{b}$ is the velocity of the droplet of the other species. However, since $v_{b}=\mathcal{O}(\partial_{c}\sigma_{b}(0))$, this contribution is also subdominant. Therefore, in our analysis of the droplet dynamics in one dimension and for weak responses, we can assume $\partial_{c}\bar{\phi}_{a,\pm}(0)=0$ without loss of generality. In higher dimensions, the analysis of systems with space-dependent binodal densities becomes more challenging, since for interfaces with finite curvature the dependence of $v_{a}$ on $\partial_{c}\bar{\phi}_{a,\pm}(0)$ is no longer subdominant. For these systems, we nevertheless set $\partial_{c}\bar{\phi}_{a,\pm}(0)=0$ to simplify the analysis, although the generalization to include this correction is not expected to change the qualitative features of the results.

We now discuss the transformations leading to the reduction of the parameter space presented in the main text. Since $\bar{\phi}_{\pm,a}$ are constants we can rescale the fields as

$$\phi_{a}\rightarrow\bar{\phi}_{a,-}+\Delta\bar{\phi}_{a}\,\phi_{a},\hskip 14.22636ptc\rightarrow\frac{r_{1}\bar{\phi}_{1,-}+r_{2}\bar{\phi}_{2,-}}{\gamma}+c,$$ (S1)

such that the newly defined density fields for the species go from $0$ outside of the droplet to $1$ inside of it, while the chemical field has been rescaled by its background value. The equation for the chemical in (1a) naturally introduces a timescale $t_{\gamma}=1/\gamma$ and a length scale $\kappa^{-1}=\sqrt{D_{c}/\gamma}$. It is then natural to rescale (1a) with respect to these two, so that both length and time become dimensionless. The free energy has a natural unit of energy ${\cal E}$ related to the energy scales of the molecular interactions, which we also use to dimensionalize the temperature. These transformations account for rescaling the coefficients of (1a) to:

$$r\rightarrow\sqrt{\frac{D_{c}}{\gamma}}r,\hskip 14.22636ptt\rightarrow\frac{t}{\gamma},\hskip 14.22636ptr_{a}\rightarrow\frac{\gamma}{\Delta\bar{\phi}_{a}}r_{a},\hskip 14.22636ptM_{a}\rightarrow\frac{M_{a}(\Delta\bar{\phi}_{a})^{2}D_{c}^{1+d/2}}{{\cal E}\gamma^{d/2}},\hskip 14.22636pt\mathcal{F}\rightarrow{\cal E}\mathcal{F},\hskip 14.22636ptT\rightarrow{\cal E}T.$$ (S2)

S..2 Ripening in multi-droplet systems

In the main text, we focused on systems composed of a single isolated droplet per species. However, when multiple droplets of the same species coexist in the mixture, Ostwald ripening can cause their sizes to evolve on timescales comparable to those of droplet motion. To illustrate this effect, we consider the simple but instructive case of two droplets of species $a$, separated by a distance $L$ in $d$ spatial dimensions. We neglect shape deformations and assume that both droplets remain approximately spherical throughout their evolution.

The interface point $\bm{\mathcal{X}}$ on droplet $I=1,2$ move with velocity $\dot{\bm{\mathcal{X}}}=\dot{R}_{I}\,\mathbf{n}(\bm{\mathcal{X}})+\dot{\mathbf{x}}_{I},$ where $R_{I}$ and $\mathbf{x}_{I}$ denote the droplet radius and centre of mass, respectively, and $\mathbf{n}(\bm{\mathcal{X}})$ is the outward normal unit vector at $\bm{\mathcal{X}}$. We denote by $\mathcal{S}_{I}$ and $\mathcal{V}_{I}$ the surface area and volume of droplet $I$. Furthermore, we define the surface-averaged force on droplet $I$ as

$$F_{I}=\frac{1}{\mathcal{S}_{I}}\int_{{S}_{I}}dS_{\bm{\mathcal{X}}}\left[\frac{d-1}{R_{I}}\,\sigma(c(\bm{\mathcal{X}}))+\mathbf{n}(\bm{\mathcal{X}})\cdot\nabla\sigma(c(\bm{\mathcal{X}}))\right].$$ (S3)

Applying $\frac{1}{\mathcal{S}_{I}}\int_{{S}_{I}}dS_{\bm{\mathcal{X}}}$ to the interface equations of Ref. [54] for $I=1,2$, we obtain

$$\begin{cases}\dot{R}_{1}\,\mathcal{S}_{1}\,G_{a}(R_{1})+\dot{R}_{2}\,\mathcal{S}_{2}\,G_{a}(L)+\dot{\mathbf{x}}_{2}\cdot\nabla G_{a}(L)\,\mathcal{V}_{2}=F_{1}-\mu,\\[4.0pt] \dot{R}_{1}\,\mathcal{S}_{1}\,G_{a}(R_{1})+\dot{R}_{1}\,\mathcal{S}_{1}\,G_{a}(L)+\dot{\mathbf{x}}_{1}\cdot\nabla G_{a}(L)\,\mathcal{V}_{1}=F_{2}-\mu,\\[4.0pt] \dot{R}_{1}\,\mathcal{S}_{1}+\dot{R}_{2}\,\mathcal{S}_{2}=0,\end{cases}$$ (S4)

where $G_{a}$ is the Green’s function of the Laplace equation, satisfying $M_{a}\nabla^{2}G_{a}(\mathbf{x})=\delta^{d}(\mathbf{x})$.

We are interested in the regime of well-separated droplets, $L\gg R_{I}$. In this limit, $\dot{\mathbf{x}}_{I}\cdot\nabla G_{a}(L)\,\mathcal{V}_{I}\ll\dot{R}_{I}\,\mathcal{S}_{I}\,G_{a}(L),$ so the centre-of-mass motion has a negligible influence on ripening. Solving Eq. (S4) gives the ripening rate:

$$\frac{d\mathcal{V}_{1}}{dt}=-\frac{d\mathcal{V}_{2}}{dt}=\frac{F_{2}-F_{1}}{\,2G_{a}(L)-G_{a}(R_{2})-G_{a}(R_{1})\,}.$$ (S5)

Equation (S5) shows that the impact of ripening on droplet dynamics depends strongly on dimensionality. For $d\leq 2$, the Green’s function $G_{a}(L)$ diverges in magnitude as $L$ increases, implying that droplets exchange only negligible amounts of mass once they are sufficiently far apart. Consequently, for $d\leq 2$, clusters of droplets can persist and travel distances much larger than their size, provided the droplets of species $a$ and $b$ start at separations comparable to their radii at $t=0$.

In contrast, for $d>2$, one has $G_{a}(L)\to 0$ as $L\to\infty$, and the ripening speed saturates at a value comparable to the centre-of-mass speed. In this regime, ripening takes place on timescales comparable to cluster formation and motion. Therefore, the dynamics described in the main text will predominantly occur for relatively larger droplets interacting with other species, as they grow at the expense of smaller droplets of their own species in their local neighbourhoods.

S..3 Droplet interaction in 3D using direct integration

We provide a direct calculation of the droplet interaction in $d=3$ and use it to determine the phase transition conditions for the diagram in Fig. 4(a). Without loss of generality (for $T=0$), we assume that the droplets are aligned along the vertical axis, such that the motion is restricted to the $z$-direction. Using Eq. (2) we obtain

$$\dot{z}_{a}=-\dfrac{9M_{a}}{4\pi R_{a}}\frac{\partial\sigma_{a}}{\partial c}\int d\Omega\cos\theta\,\partial_{z}c_{b}(\theta,\varphi),$$ (S6)

since only the chemical produced by droplet $b\neq a$ affects the velocity of droplet $a$. Droplet $b$ induces a perturbation $c_{b}$ of the chemical field, with respect to the background value given as

$$c_{b}(\bm{r})=\begin{cases}c^{\text{in}}_{b}(\bm{r})=1-e^{-R_{b}}(1+R_{b})\,\dfrac{\sinh|\bm{r}-\bm{x}_{b}|}{|\bm{r}-\bm{x}_{b}|}&\text{ for }\;\;|\bm{r}-\bm{x}_{b}|<R_{b},\vskip 8.5359pt\\ c^{\text{out}}_{b}(\bm{r})=(R_{b}\cosh R_{b}-\sinh R_{b})\,\dfrac{e^{-|\bm{r}-\bm{x}_{b}|}}{|\bm{r}-\bm{x}_{b}|}&\text{ for }\;\;|\bm{r}-\bm{x}_{b}|>R_{b}.\end{cases}$$ (S7)

Figure S1: Schematic of the three different droplet configurations, (I), (II) and (III), for $a=1$ in blue and $b=2$ in red.

The most difficult sub-case is the one in which the droplets intersect, corresponding to configuration (II). We place droplet $a$ at the origin and droplet $b$ above it. The droplets partially intersect such that a portion of the interface of $a$ lies inside droplet $b$. We denote by $\bar{\Omega}_{a}$ the solid angle subtended by the part of droplet $a$ that is immersed in $b$, and by $\bar{\theta}_{a}$ the corresponding polar angle, as illustrated in Fig. S1(b).

Let $|z|=|z_{b}-z_{a}|$ be the centre-to-centre distance. The integral in Eq. (S6) can be evaluated piecewise, using the “inside” expression for $c_{b}$ (Eq. (S7)) for $0\leq\theta\leq\bar{\theta}_{a}$, and the “outside” expression for $\bar{\theta}_{a}\leq\theta\leq\pi$. The intersection angle is given by the geometric relation

$$\bar{\theta}_{a}=\frac{z^{2}+R_{a}^{2}-R_{b}^{2}}{2|z|R_{a}},$$ (S8)

The equations of motion follow from evaluating the integrals in Eq. (S6) for each droplet and substituting Eq. (S8) to express the forces in terms of the relative distance $z$. We obtain

$$\dot{z}_{a}=-\dfrac{9M_{a}\partial_{c}\sigma_{a}(0)}{2R_{a}}\left[\int_{\bar{\theta}_{a}}^{\pi}d\theta\,\sin\theta\cos\theta\,\partial_{z}c^{\text{in}}_{b}(R_{a},\theta)+\int^{\bar{\theta}_{a}}_{-\pi}d\theta\,\sin\theta\cos\theta\,\partial_{z}c^{\text{out}}_{b}(R_{a},\theta)\right].$$ (S9)

The (I) and (III) subcases can be obtained from (S9) by setting $\bar{\theta}_{a}=-\pi$ for the smaller droplet in configuration (II), and $\bar{\theta}_{a}=\pi$ in all the other cases. We performed integral in (S9) analytically using the symbolic manipulation software Mathematica and use it to study the effective potential in (4a) for the 3D case. In configuration (III), we find that the fixed point at $x=0$ is a minimum if the following condition holds

$$\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}<\dfrac{(R_{2}+1)R_{2}^{2}\sinh R_{1}}{R_{1}^{2}[R_{1}\cosh R_{1}-\sinh R_{1}]}\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2},$$ (S10)

leading to (meta-)stable engulfed droplets. The condition for repulsive droplets in configuration (I) is

$$\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}<-\dfrac{R_{2}^{2}\sinh R_{1}[R_{2}\cosh R_{2}-\sinh R_{2}]}{R_{1}^{2}\sinh R_{2}[R_{1}\cosh R_{1}-\sinh R_{1}]}\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2},$$ (S11)

whereas for

$$\mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}>\dfrac{(R_{2}+1)R_{2}^{2}\sinh R_{1}}{R_{1}^{2}[R_{1}\cosh R_{1}-\sinh R_{1}]}\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2} \mathcal{C}_{2\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}1}>-\dfrac{R_{2}^{2}\sinh R_{1}[R_{2}\cosh R_{2}-\sinh R_{2}]}{R_{1}^{2}\sinh R_{2}[R_{1}\cosh R_{1}-\sinh R_{1}]}\mathcal{C}_{1\mathrel{\hbox{\rule[2.3pt]{4.0pt}{0.4pt}}\mkern-4.0mu\hbox{\char 41\relax}}2},$$ and (S12)

$V$ is a symmetry-breaking Mexican hat potential leading to the formation of self-propelling clusters in configuration (II). We use the conditions (S12), (S11), (S12) to build the phase diagram of the model, leading to the results in Fig. 4 of the main text.

S..4 Stochastic dynamics $d$-dimensional droplets

We now derive the expression for the stochastic force acting on the droplets centres of mass at finite temperature ($T\neq 0$) in the relevant case $\bar{\phi}_{a,+}=1$ and $\bar{\phi}_{a,-}=0$, although the analysis can be readily generalized. The first step is to determine the effective noise acting on a point of the interface of each droplet due to the stochastic current $\bm{J}_{a}^{s}=\sqrt{2M_{a}T}\,\bm{\nabla}\cdot\bm{\Lambda}_{a}$. This problem has been studied in the literature for both passive and active systems [68]. Our derivation is based on a generalization of the results in Ref. [54].

We start by inverting the Laplacian in Eq. (1a) to obtain

$$\int d^{d}{\bm{r}}^{\prime}G_{a}(\bm{r}-\bm{r}^{\prime})\,\partial_{t}\phi_{a}=\dfrac{\delta\mathcal{F}_{a}(\phi_{a},c)}{\delta\phi_{a}}(\bm{r},t)+\bar{\mu}^{s}_{a}(\bm{r},t)-\bar{\mu}_{\infty,a},$$ (S13)

where $G_{a}(\bm{r})$ is the Green’s function of the $M_{a}\nabla^{2}$ operator [i.e., the solution of $M_{a}\nabla^{2}G_{a}(\bm{r})=\delta^{d}(\bm{r})$], and

$$\bar{\mu}^{s}_{a}(\bm{r},t)=\int d^{d}{\bm{r}}^{\prime}G_{a}(\bm{r}-\bm{r}^{\prime})\bm{\nabla}_{\bm{r}^{\prime}}\cdot\left[\sqrt{2M_{a}}\bm{\Lambda}_{a}(\bm{r}^{\prime},t)\right].$$ (S14)

Comparing Eq. (S14) with the results of Ref. [54] (see in particular Sec. 2.2), we find that $\bar{\mu}^{s}_{a}$ enters as a stochastic thermodynamic force (denoted $\bar{\mu}(\bm{r},t)$ in Ref. [54]) in the interface equations. We calculate its variance using Eq. (S14) and find

$$\langle\bar{\mu}^{s}_{a}(\bm{r},t)\bar{\mu}^{s}_{b}(\bm{r}^{\prime},t^{\prime})\rangle=-2T\delta_{ab}G_{a}(\bm{r}-\bm{r}^{\prime})\delta(t-t^{\prime}).$$ (S15)

The noise acting on the centre of a droplet, $\bm{x}_{a}$, can then be derived from the noise on its interface. In 1D, this corresponds to the following system of equations

$$\begin{cases}-\dot{x}_{a}\,G_{a}(2R_{I})={\partial_{x}\sigma}(x_{a,r})+\bar{\mu}^{s}_{a}(x_{a,r},t)-\bar{\mu}_{a,\infty},\\ \dot{x}_{a}\,G_{a}(2R_{a})=-{\partial_{x}\sigma}(x_{a,l})+\mu^{r}_{I}(x_{a,l},t)-\bar{\mu}_{a,\infty}.\end{cases}$$ (S16)

Solving (S16) for $\dot{x}_{a}$ we find the following equation of motion for the droplet at non-zero temperature

$$\frac{2R_{a}}{M_{a}}\,\dot{x}_{a}=-[{\partial_{x}\sigma}(x_{a,r})+{\partial_{x}\sigma}(x_{a,l})]+2\sqrt{\frac{R_{a}}{M_{a}}T}\,\eta_{a}(t),$$ (S17)

where $\eta_{a}(t)$ is a Gaussian white noise as defined in the main text.

In higher dimensions, the expression of for the speed of the droplet’s centre is a generalization of the results obtained for the determinist case in Ref. [54]. In particular, we find

$$\frac{\pazocal{V}_{a}}{dM_{a}}\dot{\bm{x}}_{a}=-\int_{S_{a}}dS_{\bm{\mathcal{X}}}\;\nabla\sigma_{a}(\bm{\mathcal{X}})+\int_{S_{a}}dS_{\bm{\mathcal{X}}}\;\hat{\bm{n}}(\bm{\mathcal{X}})\mu^{s}_{a}(\bm{\mathcal{X}},t),$$ (S18)

where $\int_{S_{a}}dS_{\bm{\mathcal{X}}}$ denotes a surface integral over the droplet interface in which ${\bm{\mathcal{X}}}$ is the integration variable. The first integral in the r.h.s. of Eq. (S18) gives the deterministic force already discussed for the 3D case in Section S..3. The second integral gives the stochastic force on the droplet’s centre. Using (S15) we find that its correlator is given by

$$\int_{S_{a}}dS_{\bm{\mathcal{X}}}\int_{S_{b}}dS_{\bm{\mathcal{X}}^{\prime}}\;\hat{n}_{i}(\bm{\mathcal{X}})\hat{n}_{j}(\bm{\mathcal{X}}^{\prime})\langle\mu^{s}_{a}(\bm{\mathcal{X}},t)\;\mu^{s}_{b}(\bm{\mathcal{X}}^{\prime},t^{\prime})\rangle=2T\frac{\pazocal{V}_{a}}{dM_{a}}\delta_{ij}\delta_{ab}\delta(t-t^{\prime}),$$ (S19)

leading to the noise correlators reported in the main text.