Inhomogeneous surface tension of chemically active fluid interfaces

We study the dependence of the surface tension of a fluid interface on the density profile of a third suspended phase. By means of an approximated model for the binary mixture and of a perturbative approach we derive close formulas for the free energy of the system and for the surface tension of the interface. Our results show a remarkable non-monotonous dependence of the surface tension on the peak of the density of the suspended phase. Our results also predict the local value of the surface tension in the case in which the density of the suspended phase is not homogeneous along the interface.


I. INTRODUCTION
The physics of liquid interfaces is crucial in several scenarios 1 including the stability of biofilms in healthcare devices 2 , shipping 3 , self-cleaning surfaces for automotive industry 4 . Clearly, the mechanical and thermodynamic properties of the interface are ultimately determined by the concentrations and interactions among its molecular constituents. It is well known that adding additional suspended particles to a phase-separated binary mixture can alter the interface and in particular it surface tension. Interestingly, in many scenarios of actual interest, fluid interfaces are either formed by complex fluids whose structure may change in time, as it happens for antagonistic salts 5,6 or they are in the vicinity of "devices" that continuously transform some molecules into others. This is the case of biofilm formation 2,7 , where the cellular activity keeps on pumping molecules in and out the cell membrane. Similarly, phoretic colloids attains motion by catalyzing chemical reactions on their surfaces. When such system are in the vicinity of a fluid interface 8,9 , their chemical activity affects the local surface tension possibly leading to a variety of phenomena including the onset of Marangoni flows 10 .
In this contribution, we determine the dependence of the surface tension of a phase-separated binary mixture on the density of a third suspended phase. We show that the presence of this additional suspended phase affects the surface tension. In particular, when the concentration of the suspended phase is not homogeneous, the surface tension varies along the interface. We show that when the density of the suspended phase is not at equilibrium (i.e. its chemical potential is not homogeneous) the imbalance in the surface tension is not compensated for by any conservative force and hence it leads to the onset of Marangoni flows. At variance, when the suspended phase experience an external conservative potential, the system still reaches an equilibrium steady state characterized by a homogeneous chemical potential and an inhomogeneous density profile. In such a scenario, surface tension is still inhomogenous but the force induced by this inhomogeneity will be compensated for by the conservative potential hence leada) Electronic mail: squarcio@is.mpg.de b) Electronic mail: malgaretti@is.mpg.de ing to no Marangoni flows. The structure of the manuscript is as follows. In Sec. II we describe the coarse-grained model of the binary mixture and its interaction with the suspended phase. In Sec. III we compute the density profile of the suspended phase, in Sec. IV we derive the protocol to calculate the surface tension and we discuss some interesting scenario. Finally in Sec. V we derive the expression for the local surface tension and in Sec. VI we draw our conclusions.
A j e 8 s u r p H 1 R 9 2 w y d 2 6 t c V 3 E U S Y n 5 J S c E 4 9 c k g a 5 J U 3 S I p x o 8 k x e y Z v z 5 L w 4 7 8 7 H o r X k F D P H 5 A + c z x 8 x F J L n < / l a t e x i t >< l a t e x i t s h a 1 _ b a s e 6 4 = " r Q P j N j t K 1 Y c + + b j / F e 8 m V F N P 5 u I      p J x h f 5 X 7 n j i r N Y n l r J g n 1 B R 5 J F j K C T S 7 1 k 4 g N y l W 3 5 s 6 A V o m 3 I N V G 5 f u 9 / t F 5 a w 7 K n / 1 h T F J B p S E c a 9 3 z 3 M T 4 G V a G E U 6 n p X 6 q a Y L J G I 9 o z 1 K J B d V + N r t 1 i k 6 t M k R h r G x J g 2 b q 7 4 k M C 6 0 n I r C d A p t I L 3 u 5 + J / X S 0 1 4 4 W d M J q m h k s w X h S l H J k b 5 4 2 j I F C W G T y z B R D F 7 K y I R V p g Y G 0 / J h u A t v 7 x K 2 v W a 5 9 a 8 G 5 v G J c x R h G M 4 g T P w 4 B w a c A 1 N a A G B C O 7 h E Z 4 c 4 T w 4 z 8 7 L v L X g L G Y q 8 A f O 6 w 9 0 V Z I 1 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " S X t G 4 g G C q n u Y e F O C 4 n W X 5 f w 0 l Q U = " > A A A B 6 3 i c b V B N S w M x E J 3 U r 1 q / q h 6 9 B I v g q e x 6 U W 8 F L x 4 r 2 A 9 o l 5 J N s 9 3 Q J L s k W a E s / Q t e P C j i 1 T / k z X 9 j t t 2 D t j 4 Y e L w 3 w 8 y 8 M B X c W M / 7 R p W N z a 3 t n e p u b W / / 4 P C o f n z S N U m m K e v Q R C S 6 H x L D B F e s Y 7 k V r J 9 q R m Q o W C + c 3 h V + 7 4 l p w x P 1 a G c p C y S Z K B 5 x S m w h D d O Y j + o N r + k t g N e J X 5 I G l G i P 6 l / D c U I z y Z S l g h g z 8 L 3 U B j n R l l P

II. THE MODEL
In this section we outline the construction of the model. The approach we follow consists in two steps: in Sec. II A we recall some basic notions from solution theory, hence we consider a model of a homogeneous binary fluid. In Sec. II B we extend the model to the inhomogeneous setting and, in parallel, we also model the interaction between the binary fluid and an inhomogeneous solute.

A. Homogeneous mixture
As a warm up we consider a binary liquid mixture and denote its components with the labels a and b. Let N a be the number of molecules of species a and analogously for N b . arXiv:2008.13004v1 [cond-mat.stat-mech] 29 Aug 2020 We are assuming incompressibility, hence N a + N b = N, with N a constant.
At constant pressure and temperature, the system is described by the Gibbs free energy G = G 0 + ∆G, where is the Gibbs free energy of the unmixed system and µ (0) a , and µ (0) b are the corresponding chemical potentials for pure components a and b. The free energy of mixing takes the form 11 where β −1 = k B T and the (dimensionless) parameter χ measures the strength of the interaction between the binary mixture components. The chemical potentials µ j = ∂G/∂N j , with j = a, b, are given by where n a = N a /N, n b = N B /N; and evidently, n a + n b = 1.
Parenthetically, we observe that, thanks to Eq. (3), the free energy in Eq. (1) can be rewritten in the compact form G = N a µ a + N b µ b . Notice also that the mixing free energy is invariant under the exchange of N a with N b . In order to exploit the aforementioned symmetry, we parametrize the concentrations in terms of the order parameter φ ∈ (−1/2, 1/2), the latter is defined by At phase coexistence, the two chemical potentials are identical, i.e., µ a = µ b , hence by defining from Eq. (3) we find the following relationship 12 between the order parameter and ∆µ The above can be written in the equivalent form the latter is precisely the mean-field equation in the Curie-Weiss model of ferromagnetism, provided −∆µ is identified with the external magnetic field and φ with the bulk magnetization 13 .
The equilibrium composition of the binary mixture follows from the solutions of Eq. (7) which corresponds to minima of the Gibbs free energy. On expressing N a and N b in terms of φ, and by expanding in powers of φ, we find In order to get a quantitative picture, we expand Eq. (8) in powers of φ up to O(φ 4 ), thus Eq. (8) becomes thus, we have obtained the standard φ 4 double well potential for the bulk free energy. By combining Eq. (9) and Eq. (1), the total Gibbs free energy per particle becomes withn defined in Eq. (4). Let us consider the phase coexistence, thus we set ∆µ = 0. If χ < 2, then the potential Eq. (10) exhibits a (stable) global minimum only for φ = 0, while if χ > 2 a new pair of stable solutions of the form φ = ±φ 0 appears. In the unmixed regime, the equilibrium concentrations are given byn ± φ 0 with The two-phase coexistence line spanned by χ > 2 terminates in a critical point, the latter is located in χ = 2. The existence of the aforementioned solutions can also be inferred by inspecting the mean-field equation Eq. (7) with ∆µ = 0.

B. Inhomogeneous mixture
In this section we extend the treatment of the binary liquid outlined in Sec. II A to instances in which inhomogeneous densities are allowed. The typical situation is the one in which the two coexisting phases are separated by an interface. On top of that, we also allow for the presence of third substance, which will be termed the solute.
In order to describe the inhomogeneous system, we promote the number densities to be spatially varying fields, therefore we promote N j to N j (r), with j = a, b. Then we introduce the densities ρ a (r) = N a (r)/V and ρ b (r) = N b (r)/V, where ρ = N/V is a constant reference density. The total number of molecules for each species is found by integrating over the volume, i.e., N a = ∫ dr ρ a (r), and analogously for b. Having in mind a situation in which the solute is dilute, we model it through an ideal gas with weak interactions with the binary liquid. The Gibbs free energy is now given by with (r) the solute density and the ideal gas contribution. Some comments are in order: i) although we have introduced two density fields, ρ a and ρ b , the functional in Eq. (12) depends on a single density field. In fact, ρ a and ρ b are linearly dependent because ρ a + ρ b = ρ. The functional given by Eq. (12) depends on either ρ a , or ρ b , or any linear combination of them with the exception of ρ a + ρ b . The natural choice is to parametrize the densities by following Eq. (4), thus we introduce the order parameter φ(r) defined by ii) the coupling K in front of the square gradient terms in Eq. (12) quantifies the energetic cost for density inhomogeneities. Within van der Waals theory of liquid-vapor interfaces, the coefficient in front of the square gradient is proportional to the second moment of the inter-particle interaction potential; taking the opposite sign, hence the coefficient is positive 14 . Since incompressibility imposes ρ a (r) + ρ b (r) = ρ, the two square gradients in Eq. (12) actually coincide. iii) in principle, an external potential V(r) can be coupled to the solute density. The weak interaction between the solute and the binary liquid mixture is taken into account by the coefficients ω a and ω b . Notice that, if ω a = ω b , then the solutebinary interaction does not depend on the concentration φ of the binary liquid, meaning that the solute and the binary liquid are decoupled. On the other hand, if ω a ω b , then the solutebinary interaction becomes proportional to φ(r) (r). We note that the presence of the potential V(r) in Eq. (12) makes it different from those models that accounts for a inhomogeneous temperature profile 15 .
iv) the homogeneous model described in Sec. II A is retrieved in the double limit ρ a , ρ b → constant and ω a , ω b → 0. On the other hand, by turning off the interaction with the solute, and keeping in mind that ρ a and ρ b are not independent, Eq. (12) reduces to the Cahn-Hilliard effective Hamiltonian 16 .
v) Since we are interested in inhomogeneous systems, we consider the generalized chemical potentials which are given by 17 .
Again, by taking the limit of homogeneous system with ω a = ω b = 0, we retrieve the results given in Eq. (3). By following the same guidelines of Sec. II A, in order to find a coarse-grained model, we expand the free energy functional in powers of the density fields φ(r) and (r) by keeping the lowest powers. The result of this procedure leads to the free energy functional where V (r) is the effective external potential As a consistency check, we observe that for a homogeneous system in absence of the solute, the functional in Eq. (17) reduces to Eq. (9), provided N is identified with ∫ dr ρ. The difference of generalized chemical potentials ∆µ φ (r) = µ a (r) − µ b (r) follows straightforwardly from Eq. (16) (19) Alternatively, Eq. (19) follows by taking a functional derivative of the free energy with respect to the order parameter, i.e., It is convenient to factorize the coefficient K ρ 2 in Eq. (17), and to introduce the rescaled bulk free energy and the rescaled coupling parameter Thanks to the above definitions, the free energy functional Eq. (17) can be written as follows with the effective Hamiltonian Notice that Eq. (24) is formally identical to the Landau-Ginzburg free energy for a uniaxial ferromagnet with order parameter field φ(r) in the presence of an inhomogeneous bulk field h(r) = −ω (r) 18,19 .

C. Equilibrium
At thermodynamic equilibrium, the chemical potentials of both components are equal, therefore their difference given by Eq. (20) vanishes. Moreover, the chemical potential of the solute is constant. A simple calculation gives and therefore the solute density is given by a Boltzmann distribution, which we write in the compact form where is a reference density set by the chemical potential µ , while is the total external potential acting on the solute. Notice that V(r) receives a contribution also from the binary liquid provided ω a,b 0.
Concerning the density of the binary mixture, the latter turns out to be found as the solution of which follows by equating the two generalized chemical potentials in Eq. (16). Thanks to Eq. (19), the above gives We observe that a nonzero solute density (r) acts as an external source for the binary liquid density φ(r). As a result, the binary liquid density changes accordingly due to the presence of the solute. The solution of Eq. (31) will be discussed in Sec. III. We conclude by noting that in the absence of the solute-binary interaction, i.e., ω a = ω b , the right hand side of Eq. (31) vanishes and the density profile of the binary liquid is given by the classical van der Waals theory 14 .

D. Out of equilibrium
Let us consider now a non-equilibrium regime in which the solute is produced within the binary mixture, while the latter is assumed to be at local equilibrium for the assigned solute density (r). Within this scenario, the production rate of the solute is considered to be rather slow, thus the binary liquid adjusts its composition by following the solute density according to Eq. (31). Accordingly, the binary mixture is at equilibrium with the specified solute density profile and hence there are no density currents within the binary mixture, i.e., J φ ∝ ∇∆µ φ = 0, the latter follows from ∆µ φ = 0. The density φ(r) satisfies Eq. (31), however the source term ∝ (r) has to be identified, as we are going to show.
Contrary to the equilibrium setting discussed in Sec. II C, the chemical potential of the solute is no longer constant, therefore it exists the net current of solute J (r) = −L ∇βµ (r) ; (32) within linear response theory, the Onsager coefficient reads L = D where D is a diffusion coefficient, thus Eq. (32) becomes The time evolution of the solute density then follows accordingly with the continuity equation if no reactions are taking place, and if the solute is not produced/destroyed within the system, i.e., ∂ t = −∇ · J . More generally, however, the solute is produced/destroyed due to certain chemical reactions whose occurrence is confined within small region of the system; hence sources and sinks have to be included in the description. We thus consider the following equation for the time evolution of the solute density The terms on the right hand side of Eq. (34) correspond to the density of sources, S(r), and sinks, (r)/τ, where we have assumed a uniform decay of the solute concentration with characteristic time τ. It is understood that S(r) may take both signs, thus a local sink would correspond to S(r) < 0 in certain regions of space. We will see an example of this sort in Sec. III B. By combining Eq. (33) and Eq. (34), the non-equilibrium stationary state is given by the solution of D∇ 2 (r) + βD∇ · (r)∇V(r) = (r)/τ − S(r), with asymptotic boundary condition → 0 at |r| → ∞ (away from sources/sinks) and with a specific boundary condition on the sources/sinks. We emphasize that Eq. (35) and Eq. (31) form a coupled system of differential equations for the densities φ(r) and (r). In particular, from Eq. (31) we see how the density φ(r) decouples from (r) provided ω a = ω b . Analogously, if ω a = ω b , then the total potential V(r) does not depend on φ(r) and, correspondingly, also Eq. (35) decouples from φ(r).
The above considerations suggest to seek for a solution for the density fields φ(r) and (r) by means of a perturbative analysis based on the parameterω, which is assumed to be small. In order to illustrate the approach, we consider the absence of an external potential, i.e., V(r) = 0, and write the solution of Eq. (35) as follows in which k (r) = O(ω k ). At leading-order in powers ofω, Eq. (35) reads The above will be solved in Sec. III B for a specific choice of the solute source S(r).

A. General formalism
In this section we analyze the binary liquid mixture in the phase separated regime in the presence of an external field which, within our formalism, is represented by the density of solute. In particular, we show how to compute the density φ across a domain wall separating two coexisting phases. We employ the so-called double parabola approximation, which consists in replacing the φ 4 double well with two parabolic branches, i.e., The latter has been successfully employed in the study of short range critical wetting in three dimensions 18,20,21 , for which a wealth of formal results can be applied to the model studied in this paper.
To be definite, we consider a phase-separated binary mixture with the component of higher density standing below the interface, as depicted in Fig.1. The mean-field equation governing the density profile in the presence of the external field can be found by imposing the stationarity of the functional given by Eq. (24) in which φ Ξ (r) denotes the solution of Eq. (39) with the appropriate boundary conditions. The advantage of the doubleparabola approximation is that it allows for an analytical treatment, because Eq. (24) gives a pair of linear partial differential equations, which can be solved by means of Green's functions techniques. Eq. (39) gives where Ω a denotes the region occupied by the a-rich phase, i.e. φ > 0, while Ω b denotes the region occupied by the brich phase, i.e. φ < 0. The above differential equations have to be supplemented with the asymptotic boundary conditions φ(r) → ±φ 0 for z → ∓∞. Another pair of boundary conditions follows by imposing φ(r ) = 0 along the interface, the latter corresponds to the locus of points in which the domains Ω a and Ω b come in touch. According to the crossing criterion mentioned above, the interface is formally defined by We specialize our attention to a flat interface parallel to the (x, y) plane, and let z = be the interface location. The regions Ω a and Ω b are respectively given by the half-spaces z < and z > , while the plane z = defines the locus of the interface. In order to solve the non-homogeneous problem specified by Eq. (40), we first consider the homogeneous one ( = 0). The rescaled Green function K(r, r ) satisfies where κ = ξ −1 b is the inverse bulk correlation length. In three spatial dimensions, the above gives the Ornstein-Zernike kernel 18 By virtue of translational invariance along the (x, y) directions, the three-dimensional problem can be reduced to a onedimensional problem. By integrating Eq. (43) over the longitudinal coordinates (x and y), we obtain the rescaled Green's function which satisfies the one-dimensional version of Eq. (42) The solution of Eq. (40) with asymptotic conditions φ(z → ±∞) = ∓φ 0 enforcing a domain wall, can be split as follows where φ Ξ (z; 0) is the density profile for the homogeneous problem, i.e., when = 0, which reads where defines the interface position. The density profile given by Eq. (47) is plotted in Fig. 2. The interface profile given in Eq. (47) smoothly interpolates between the bulk phases ±φ 0 . The typical variations in density occur on the length scale set by the bulk correlation length. We also point out that is not fixed for an isolated interface standing in the bulk; this is a natural consequence of invariance under translations along the z axis. The archetypal situation in which the translational invariance is broken is provided by the presence of a wall. Such an instance occurs within the study of phase separation in confined geometries and wetting transitions. There the interface height z = with respect to a wall located in z = 0 turns out to be determined by the adsorption properties of the wall; see for instance 22 . In our problem instead, will be determined by a free energy minimization compatible with a canonical criterion in which the mass of the binary liquid is conserved (see Sec. III B).
Coming back to Eq. (46), the term ∆φ Ξ (z; ) gives the devi-ation from the homogeneous solution. The latter satisfies the inhomogeneous problem with boundary condition implied by the crossing criterion (φ Ξ ( ; ) = 0) together with the decomposition Eq. (46). The solution of Eq. (48) subjected to Eq. (49), is given by In order to prove that Eq. (50) is the correct solution, it is enough to apply the operator −∂ 2 z + κ 2 by using Eq. (45) and observe that the term with two K-kernels vanishes because z . For z > the second square bracket vanishes, while the first one gives (z). A similar reasoning applies to the case z < , now with the first square bracket vanishes while the second one gives (z).

B. Density profiles due to a source-sink doublet
We are now in the position to apply the ideas developed so far. The interface hosts chemical reaction with the solute which, in turn, is produced on one side of the interface and then it disappears on the other. This specific instance is modeled by introducing sources of solute on one side of the interface and sinks on the other one; hence a source-sink doublet.
In order to single out the essential features of the model, we consider a flat interface. Consequently, any density inhomogeneity will depend only on z. By virtue of translational invariance along the directions x and y, Eq. (37) gives with the source-sink doublet modeled via The magnitude S 0 > 0 corresponds to a source localized in z = a and a sink localized in z = −a. It is clear that a simultaneous sign reversal of both S 0 and a is a symmetry of Eq. (52). The specific functional form of localized source/sinks given by (52) has been used for mathematical convenience. It is indeed simple to extend the calculation to the case in which Eq. (52) is replaced by a smooth function, for instance given by a Gaussian profile centered in z = a and another one with opposite amplitude centered in z = −a.

By introducing the Fourier transform
with inverse and analogously for S(z), the density of solute in the stationary state is formally given by The Fourier transform of Eq. (52) follows straightforwardly therefore Eq. (55) immediately leads us to with penetration depth and overall amplitude¯ 0 = S 0 λ/(2D). Let us consider the interface separating the components of the binary mixture and the density profile across it. At order O(ω 0 ) the density profile is the one corresponding to a free interface, the latter is given by φ Ξ (z; 0); see Fig. 2. At order O(ω) the density profile is given by notice that in the above we have plugged the leading order solution 0 in the correction term ∆φ(x; ), thus that contribution is of order O(ω). By plugging the solute profile given by Eq. (57) into Eq. (59), we obtain a family of profiles parametrized by the interface location . The position of the interface is obtained by imposing the mass conservation of the two species. In our model this implies that Such a constraint is satisfied in the absence of solute by the free profile φ Ξ (z; 0) by taking = 0. Setting = 0 in the presence of the solute still gives a vanishing total mass because the resulting profile is odd with respect to z = 0. Therefore the interface of the binary liquid exhibits the profiles shown in Fig. 3.  Fig. 2 only within a finite distance within the interfacial region. The above feature is actually a direct consequence of the fast decay which characterizes the solute density (z); see Eq. (57). Notice also how the binary-solute interaction can lead to both monotonic (ω a < ω b ) and non monotonic (ω a > ω b ) density profiles.

IV. FREE ENERGY IN THE PRESENCE OF AN EXTERNAL FIELD
In this section we examine the energetic aspects by computing the free energy for the binary liquid in the presence of the solute. It is instructive to begin by considering the case in which the solute is absent 14 . Then we will introduce a non vanishing solute concentration and we will examine its consequences on the free energy.

A. Surface tension in the presence of an inhomogeneous field
In the absence of external fields, the surface tension is defined as the difference -per unit area -between the free energy of the system with an interface and the free energy of the system without the interface. The excess free energy (in k B T units) is thus given by where G[φ 0 ; 0] is the free energy of the system characterized by the uniform density φ 0 . For a flat interface parallel to the xy-plane we can perform the integrations along x and y, this gives the cross sectional area A π = ∫ dxdy, which factorizes. Hence the surface tension is given by The calculation of Eq. (62) is a standard exercise which we briefly recall. Equation (62) is actually the functional in which φ Ξ is the profile Eq. (47). It is well known that for "one-dimensional" interfaces the above equation admits a mechanical interpretation in which the energy minimization is equivalent to Hamilton's principle of least action 14 . From this mechanical analogy it follows that Eq. (63) admits a first integral. As a result, the contribution stemming from the square gradient -the kinetic term in the mechanical analogy -and the contribution due to the potential in Eq. (63), actually coincide. Moreover, one can show that which implies that U(φ) suffices for the determination of the surface tension. By using Eq. (64) it is thus immediate to obtain If we would have used the φ 4 double well we would have obtained the same result up to a factor 2/3. It is worth notice that for an effective Hamiltonian which is invariant under the global transformation φ ↔ −φ, then the excess free energy with coincides with the one appearing in Eq. (61) if we neglect a set of zero measure in which the gradient is ill defined. The term subtracted in Eq. (66) corresponds to the bulk free energy of a uniform system with density φ 0 for z < and density −φ 0 for z > . Since that δF Ξ = δ F Ξ , both Eq. (61) and Eq. (66) lead to the same value of the surface tension. In order to construct the notion (which is not unique) of excess free energy in the presence of an external field, we adapt the definition Eq. (66) as follows thus the excess free energy per unit area (per k B T) reads with We observe that K(z, ) = exp(−|z − |/ξ b ) vanishes exponentially far away from the interface position. It thus follows that both terms in (A8) are always finite since the solute density is expected to be bounded far away from z = .
B. Modified surface tension due to a source-sink doublet As a concrete example, we compute the deviation from the unperturbed surface tension γ 1 [ ] for the source-sink doublet we discussed in Sec. III B. It is convenient to write γ 1 [ ] in units of γ 0 and express the result as follows where g is the dimensionless parameter which measures the strength of the solute-binary liquid coupling, while Σ(κa, κλ) is a dimensionless scaling function of the scaling variables κλ = λ/ξ b and κa = a/ξ b ; we recall that ξ b is the bulk correlation length. A simple calculation gives Notice that Eq. (74) is an odd function of the distance, a, of the source/sink from the interface. The scaling function Eq. (74) is plotted in Fig. 4. For a fixed penetration depth λ and correlation length κ (see Eqs. (57) and (58)), the modified surface tension γ 1 [ ] exhibits a maximum when the separation between sources and sinks attains the value The latter gives rise to the maximum (for g > 0) Regarding the above as a function of κλ, it is simple to show that it reaches the maximum value 8g/e ≈ 2.9g for κλ = 1. The maximum deviation of the surface tension is thus reached when the penetration depth matches the bulk correlation length. The explicit solution shown above indicates that γ 1 [ ] becomes a strong effect only within a certain window of parameters. The latter is presumably a more general feature due to the weakness of the binary-solute interaction. This scenario is in clear contrast with the effects generated by surfactants. We conclude this section by examining the force required to be exerted on the source-sink doublet in order to keep them in equilibrium. The free energy of the system is F = K ρ 2 H [φ Ξ ; ]. By using Eq. (A10) and keeping terms up to O(ω), we have with γ r given by Eq. (A8). It is convenient to isolate those terms which depend on the source, thus we write F = γ 0 A π 1 + gF with F = (γ 1 [ ] + γ r [ ])/γ 0 , hence the interesting piece to analyze is F. We have the scaling form analogous to (74) The function F(κa, κλ) is odd with respect to a and it is characterized by a sigmoidal profile which interpolates between ±4κλ for a → ∓∞. Since F(κa, κλ) is a monotonous function, it reaches its minimum value for κa → +∞. For g > 0 the free energy is minimum when a → +∞, while for g < 0 the free energy is minimum when a → −∞. The above considerations lead us to conclude that the source-sink doublet needs an external force f in order to keep a constant position. The force turns out to be given by hence Since F(κa, κλ) decreases monotonically upon increasing a, the force f (κa, κλ) is always positive; clearly it is an even function of a; see Fig. 5. The force attains a certain nonzero value when a = 0, it decreases for κ|a| 1 and reaches the maximum value for κa ≈ 1. Let us take g > 0, meaning that the solute is produced in the less dense phase of the binary and it is destroyed in the denser phase of the binary. The force turns out to be positive, thus the separation between source and sink tends to increase.

V. LOCAL SURFACE TENSION
So far we have considered an external field whose intensity (z) does not change along the longitudinal directions parallel to the interface. Here we go one step beyond this hypothesis by considering an external field which may weakly change along the directions parallel to the interface, i.e., x and y. In order to get analytical insight, we assume a slow variation of along the parallel direction which occurs on a length scale λ much larger than the penetration length λ, hence λ λ. Accordingly, the Laplacian operator in Eq. (37) can be decomposed as follows ∇ 2 = ∇ 2 + ∂ 2 z . The separation of scales we are assuming implies We seek for a solution of Eq. (37) in the form of an expansion in the small parameter ≡ λ/λ , hence we write Due to the explicit dependence on the coordinates x , we cannot any longer factor out the cross sectional area, therefore the free energy reads where γ[ ] is the free energy per unit area, as defined by Eq. (77), i.e., F /A π for the translationally invariant case. Accordingly, we have where For a mildly undulated interface -thus in absence of overhangs -it is possible to use the Monge parametrization. The latter allows us to describe the interface in terms of the height function z = (x ), thus the area A [ ] of the interface is given by For weak curvatures, the latter can be approximated by means of the square-gradient term where in the last line we exploit the smallness of the external field which allows to expand the square root in Eq. (86) and identify the first term as the area A π for a flat interface. Therefore, up to quadratic corrections in the radii of curvature of the interface, the free energy per unit cross sectional area is now given by the ratio F /A π . The above result may be interpreted as a local surface tension; the latter can be identified with The functional γ 1 [ ] encodes the nontrivial spatial dependence of the local surface tension.
A. At and away from equilibrium The notion of local surface tension provided by Eq. (87) needs to be taken carefully. In order to understand its implications, we will discuss the equilibrium and the non equilibrium settings.
Away from equilibrium. Within the assumption of local equilibrium for the binary mixture, the binary liquid relaxes on time scales which are fast compared to the dynamics of the solute. The equality of generalized chemical potentials, i.e., ∆µ φ (r) = 0, is used in order to find the density of the binary liquid. In particular, if the solute density is invariant under translations along the interface ( = (z)), then the chemical potentials of both components of the binary liquid are equal across the interface, thus µ a (z) = µ b (z). On the other hand, if the solute density (x , z) breaks the invariance under translations along the interface, then the chemical potentials of the a and b components may acquire a dependence on the coordinates x , but still within the assumption of local equilibrium The latter means that for every liquid column with cross sectional area dxdy centered in z, the system is at local phase coexistence. At equilibrium. Being at mechanical equilibrium, the sum of all stresses at the interface must vanish. De facto, Eq. (87) supports the existence of non-vanishing local surface tension gradients at equilibrium with and external field. On such a basis one would conclude that the resulting inhomogeneous local surface tension will produce tangential stresses which will move the fluid, and eventually they lead to an incompatibility with the mechanical equilibrium. Although the local surface tension given by Eq. (87) is inhomogeneous, the corresponding Marangoni stresses have to be inevitably compensated. As we are going to show, such a compensation requires a non-trivial condition for solute density.
At equilibrium we further require the chemical potentials µ a and µ b to be equal, but also their gradient along the interface has to be identical. Therefore Eq. (89) is supplemented with The solute density has to satisfy Eq. (90) in order to have absence of flows at equilibrium. For the sake of completeness quote the expressions for the individual chemical potentials of both the binary liquid components We conclude by recalling that for the φ 4 model, the potential W(φ) is given by Eq. (21) with φ 0 given by Eq. (11) and K ρκ 2 = 4( χ − 2). For the doubleparabola model instead one uses Eq. (38) where κ is the same as in Eq. (92).

VI. CONCLUSIONS
We have analyzed the dependence of the free energy and of the surface tension of a fluid interface on the density of a third suspended phase. By means of a simplified model for the free energy of the unperturbed interface and via a perturbative approach we have derived a close formula for both the free energy and the surface tension.
First, we have specialized to the case in which the density of the suspended phase is homogeneous along the interface. In this scenario our results show that, even at lowest order in the perturbative approach, the surface tension is sensitive to the presence of a third suspended phase. In particular, for the case-model that we have studied, our results predict a nonmonotonous dependence of the surface tension on the distance of the source of the suspended phase.
Second, when the density of the suspended phase is not homogeneous along the interface, our results predict a dependence of the surface tension on the local value of the density of the suspended phase. If the systems is at equilibrium, i.e. such an inhomogeneity is due to some external field, then the local force associated with the inhomogenous surface tension will be balanced such that there will be no flow. In contrast, if the inhomogenous density of the suspended phase occurs in a non-equilibrium scenario, then the gradient of the surface tension will lead to the onset of Marangoni flows 10 .
Finally, our results can open the route to control diverse phenomena. For example, if catalytic colloids accumulate at one edge of the droplet, the imbalance in the density profile they generate will induce a net Marangoni flow that can induce net motion of the droplet [23][24][25] . Similarly, in the case in which the suspended phase attaches to the solid substrate it will induce an inhomogeneous solid-liquid surface tension that can pull droplets uphill 26 . the stability of Pickering emulsions upon tuning the release of an additional suspended phase.