Electrostatic pair-interaction of nearby metal or metal-coated colloids at fluid interfaces

In this paper, we theoretically study the electrostatic interaction between a pair of identical colloids with constant surface potentials sitting in close vicinity of each other at a fluid interface. By employing a simplified yet reasonable model system, the problem is solved within the framework of classical density functional theory and linearized as well as nonlinear Poisson-Boltzmann (PB) theory. Apart from providing a sound theoretical framework generally applicable to any such problem, our novel findings, all of which contradict common beliefs, include the following: first, quantitative as well as qualitative differences between the interactions obtained within the linear and the nonlinear PB theories; second, the importance of the electrostatic interaction between the omnipresent three-phase contact lines in interfacial systems; and third, the occurrence of an attractive electrostatic interaction between a pair of identical metal colloids. The unusual attraction we report on largely stems from an attractive line interaction which although scales linearly with the size of the particle, can compete with the surface interactions and can be strong enough to alter the nature of the total electrostatic interaction. Our results should find applications in metal or metal-coated particle-stabilized emulsions where densely packed particle arrays are not only frequently observed but are sometimes required.


I. INTRODUCTION
Since its discovery in the early twentieth century by Ramsden [1] and Pickering [2], stabilizing emulsions using colloidal particles instead of surfactants has become a standard technique. This new class of emulsions are popularly known as Pickering emulsions and finds application in diverse areas such as food [3][4][5], cosmetic [3,6,7], petroleum [8], pharmaceutical or biomedical [9][10][11][12], and optical industries [13], among others. Major advantages that make Pickering emulsions in general preferable over conventional surfactant-stabilized emulsions are its enhanced stability and reduced toxicity [14,15]. Another important aspect is the availability of a wide range of particles differing in their composition and shape [16][17][18]. Among them, metal colloids are gaining increasing attention owing to their numerous applications facilitated by the advancement of opto-electronics and nanotechnology in the twenty-first century [19].
Self-assembly of metal particles at fluid interfaces together with distinct optical, electrical, and catalytic properties of metal particles are exploited in liquid like mirrors [20,21], sensors [22][23][24], detectors [25], filters [26], antennas [23], controllable and targeted drug delivery equipment [9][10][11][12], plasmonic rulers [27], or in purification processes [28]. Due to their antimicrobial effect, silver particle stabilized emulsions are suitable for biomedical and textile applications [29]. For many of these applications, particularly for optical applications, a regular array of particles is needed which can be easily achieved at fluid interfaces. Contrary to the topological defects commonly observed for solid-liquid interfaces, the fluidic nature of these interfaces allows them to self-heal without any external influence, leading to strikingly uniform films spanning over large areas [30,31]. Moreover, due to their defect-free nature such structures are easy to reproduce and being fluidic, they are easily deformable [31]. * majee@is.mpg.de Even a nanoparticle, while at fluid interfaces, typically feels a strong trapping potential several orders of magnitude larger than the thermal energy, originating from the reduction of the fluid-fluid interfacial area [32][33][34]. This restricts the motion of the particle to be only along the interfacial plane. Consequently, the arrangement of particles within a monolayer and the stability of the overall structure are largely dictated by their lateral interaction. For charged particles, one of the dominant contributions to the total lateral interaction usually comes from the electrostatic origin [35][36][37][38], which is the focus of the present study. If the inter-particle separations are large, then, from the electrostatic point of view, the system behaves like a set of interacting point dipoles as the counter-ion cloud surrounding each particle differs in the two adjacent fluid media and leads to an effective dipole normal to the interfacial plane [39]. The large separation distance also allows for a point-particle assumption and for an analytical solution using the linearized PB theory or the Debye-Hückel (DH) theory [40]. Whereas this simple yet useful picture has sparked the interest of a bunch of subsequent studies [37,[41][42][43][44][45][46][47][48][49], little is known about the opposite limit, i.e., when the inter-particle separation is small compared to the size of the particles. Clearly, the dipolar assumption fails in this situation.
Most of the few recent efforts [50][51][52][53][54] targeted toward addressing the short inter-particle separations have considered particles with constant surface charge densities, which, at a simplistic level, describes dielectric or insulating particles the best [55,56]. Conductive metal particles, on the other hand, are characterized by constant surface potentials [56]. In a very recent publication [57], the interaction between particles at an air-water interface with constant surface potentials only in the portions immersed in water is addressed using the linearized PB theory. However, equipotential metal particle surfaces carry the same constant potential irrespective of the adjacent fluid phase. Not only that, metal particle stabilized emulsions often feature significantly polarizable oils [19,31,58] which behave quite differently compared to air. Moreover, the use of the linearized PB theory at short separations may not be accurate as well [51]. Therefore, a proper description of the electrostatic interaction between a pair of metal particles situated very close to each other at a fluid interface is still lacking and we aim toward filling this gap here.
The short separation situation we consider is frequently encountered experimentally for metal particle-stabilized emulsions [13,19,24,31,[58][59][60][61]. In fact, for nanoplasmonic mirrors or detectors, in order to get substantial signal, one needs densely packed array of relatively big particles ( 40 nm in diameter) with average surface to surface distance being smaller than the particle radii [13,19,25,62]. Not only that, our study could be equally useful for interfacial Janus colloids [63][64][65] when the metal caps face each other or for core-shell particles consisting of metal shells with metallic or cost-effective non-metallic cores [24,66]. One should also note that, strictly speaking, all particles including the metal ones in electrolyte solutions are charge-regulated [67,68]. We do not account for such complexities here. However, the constant potential case should still remain insightful as the charge regulation solution is bounded by constant charge and constant potential limits [56].
In order to tackle the problem efficiently, we use some justified approximations (see Fig. 1). First, the short inter-particle separation is exploited to treat the particles as flat plates by ignoring their curvatures. This assumption, motivated by the Derjaguin approximation, has been used before in this context as well [50,51,53]. Second, we assume the fluid-fluid interface to be flat as interfacial deformations are usually negligible for smooth particles up to a few micrometers in size [35,[69][70][71]. Moreover, we consider a 90 • liquid-particle contact angle which corresponds to the maximum reduction of interfacial area for any spherical particle and has been particularly shown to be pivotal for the entrapment of nanoparticles [72,73]. Several other existing theoretical models also rely on these latter two assumptions [45,57,[74][75][76].
The electrostatic problem for the resulting model system, as depicted in Fig. 1(b), is studied within the framework of classical density functional theory (DFT) which requires the free energy of the system as the only input. A subsequent minimization of this free energy then leads to the governing equation for the electrostatic potential along with the boundary conditions and the ensuing effective interaction, which automatically ensures self-consistency, can also be attained easily. We solve the problem separately using the linearized PB theory and the nonlinear PB theory. In each case, the effective interaction is decomposed into parts proportional to the surface areas of the plates and to the lengths of the three-phase contact lines. We also provide results for the separation-independent interactions, i.e., surface, interface, and line tensions present in the system. Apart from offering a comparison, the analytically obtained results within the linear theory also serve as checks for those obtained numerically within the nonlinear PB theory in proper limits. Our results feature several novel aspects that contradict common beliefs. First, we show that the interaction of the three-phase contact lines plays a crucial role in determining the nature of the total electrostatic interaction. This is not intuitively obvious as the line part to the total electrostatic Two identical particles (represented by the yellow circles) with constant surface potentials Ψ P are sitting next to each other at a flat interface (indicated by the green horizontal line) formed by two immiscible electrolytes denoted as medium "1" (for x < 0) and as medium "2" (for x > 0). The surface-to-surface separation between the particles is small compared to their size and each particle is submerged equally in the two fluid phases corresponding to a liquid-particle contact angle of 90 • . (b) Enlarged view of a model system that approximates the boxed region in part (a). Owing to their short separation, the particles are assumed to be parallel flat plates situated at z = 0 and z = L and aligned vertically with respect to the fluid interface. Fluid phase "1" ("2") has dielectric constant ε 1 (ε 2 ). The ionic strength of the added binary monovalent salt and the corresponding inverse Debye length in medium "1" ("2") are given by I 1 (I 2 ), and κ 1 (κ 2 ), respectively. interaction scales linearly with the size of the particle whereas the surface parts scale with the square of the size of the particle. Second, our findings suggest that, depending upon the specific system under consideration, the interaction obtained within the linearized PB theory can differ quantitatively as well as qualitatively from that obtained within the nonlinear PB theory. While the quantitative differences are not unexpected at small inter-particle separations, the qualitative one is indeed surprising. Third and as a most striking result, it is shown that the electrostatic interaction of identical metal particles at fluid interface is not necessarily always repulsive. The unexpected attraction stems from the line interaction energy which can be repulsive as well as attractive and, depending upon the system, can be strong enough to alter the nature of the total electrostatic interaction.

II. MODEL AND FORMALISM
In a Cartesian coordinate system, the region of interest is the space x ∈ (−∞, ∞), z ∈ [0, L] bounded by two plates positioned at z = 0 and z = L and filled with two immiscible fluid phases creating an interface at x = 0 (see Fig. 1(b)). As the plates mimic metal surfaces, they are equipotential. Moreover, they are considered to be identical and are modeled as constant potential surfaces, i.e., they carry the same surface potential Ψ P irrespective of the separation distance L between them. The fluid phase occupying the space x < 0 (x > 0) is denoted as medium "1" ("2"). The dissolved salt is a simple binary compound composed of oppositely charged monovalent ingredients only. Consequently, the ionic strength equals the salt concentration and its bulk value in medium i ∈ {1, 2} is given by I i . The background solvents, i.e., the fluids are treated as structureless incompressible linear dielectrics with dielectric constants ε i = ε r,i ε 0 , i ∈ {1, 2}, where ε r,i is the relative permittivity of medium i and ε 0 is the vacuum permittivity. Please note that the length scale of our interest is the Debye screening length κ −1 i = ε r,i / (8π B I i ), where B = e 2 / (4πε 0 k B T ) denotes the vacuum Bjerrum length with the elementary positive charge e, the Boltzmann constant k B , and the temperature T . Therefore, any phenomena occurring on a smaller length scale such as the scale of bulk correlation length or molecular length are discarded. This includes the structures formed in the liquids due to the presence of the salt ions or of the surfaces and the interface, and associated changes in the ion number density profiles n ± (r) [77]. Thus, both the ionic strength profile I(r) and the dielectric constant profile ε(r) vary steplike at the interface. However, we consider the variation in the local charge density e [n + (r) − n − (r)] as it varies on the scale of the Debye length. On this length scale, the salt ions can be considered as point like objects and, following the standard practice within a mean-field theory, we ignore any ion-ion correlation. With all these, the grand canonical density functional in the units of the thermal energy k B T = 1/β for our system, which is in equilibrium with the ion reservoirs provided by the bulk of the two fluid media, is given by In this expression, within the curly brackets, the first two terms correspond to the entropic ideal gas contribution of the ions with fugacities ζ ± , the third term describes the contribution due to ion solvation expressed via an external potential V ± (r) acting on the ions. The term quadratic in the electric displacement vector D includes all the Coulomb electrostatic interactions in the system due to the presence of surface charges and ions, whereas the last term stands for the work done by the system to maintain the plates at a constant surfaces potential Ψ P . The integration volume V is the space (x, z) ∈ (−∞, ∞) × [0, L] enclosed by the bounding surfaces ∂V, i.e., the two plates at z = 0 and L withν (r) denoting the unit outward normal to these surfaces. Although Eq. (1) is central to both the linearized and the nonlinear PB theory and what we do next is in principle the same within both the theories, the ways in which we proceed from here differ slightly as the former is analytically tackleable whereas one needs to resort to numerical techniques for the latter. Hence, below we discuss them separately.

A. Linear theory
Within the linearized PB theory, one assumes that the deviations in the ion number densities compared to their bulk values are small. Consequently, Ω [n ± ] is expanded in terms of these small deviations to obtain Ω [n ± ] by retaining terms up to quadratic order in the expansion. A minimization of Ω [n ± ] with respect to n ± then provides the equilibrium profiles n eq ± which, together with the Gauss's law, leads to the linearized PB or DH equation to be solved for the electrostatic potential Ψ i (r) in medium i ∈ {1, 2} (see Chapter 2 of [78] for details). In conjunction with this, the following boundary conditions, which also come out of the minimization process, must be satisfied. First, both the electrostatic potential and the x-component of the electric displacement vector should be continuous at the fluid-fluid interface, i.e., at any given z-value, Second, the electrostatic potential must match the surface potential of the two plates at z = 0 and L: Moreover, the electrostatic potential needs to be finite while approaching x → ±∞ which, in fact, is a prerequisite for using a Derjaguinlike approximation and is typically satisfied due to electrostatic screening. Ψ b,i in Eq. (2) refers to the bulk electrostatic potential in medium i ∈ {1, 2}. By construction, the bulk electrostatic potential profile Ψ b (r) reads as This difference (Ψ D ) arises from contrasting solvation energies of the ions in the two media and is known as the Donnan potential or Galvani potential difference [79].
As shown in Sec. III A, Eq. (2) is analytically solvable for our set-up. Once the electrostatic potentials Ψ i (r) are known, they are used to calculate n eq ± [Ψ] by considering the ion density profiles as functionals of Ψ. Finally, inserting these n eq ± [Ψ] back into the expression for Ω [n ± ], one obtains the grand potential Ω (L) = Ω [n ± [Ψ]] of our system which combines the following contributions distinctly different from each other according to their origin: Here, Ω b,i = −2I i /β is the osmotic or entropic energy contribution due to the ideal gas formed by the ions expressed per volume V i of medium i ∈ {1, 2}, γ i is the surface tension acting between each plate and medium i, ω γ,i (L) is the surfaces interaction energy density between the portions of the plates immersed in and acting through medium i, A i is the total surface area of the two plates immersed in medium i, γ 1,2 is the interfacial tension between medium "1" and medium "2", A 1,2 is the total area of the fluid-fluid interface; τ is the line tension present at the three-phase contact lines at (0, 0) and (0, L), and ω τ (L) is the contribution due to interaction between the contact lines expressed per total length of the two contact lines. Clearly, γ i and τ are L-independent as they result from the interaction of a single plate with its adjacent fluid(s). For the purpose of calculating the effective interaction between two colloids as shown in Fig. 1(a), the terms involving Ω b,i and γ 1,2 in Eq. (3) also become irrelevant as the total volume of the fluids and the total area of the fluid-fluid interface, which not only include the boxed region in Fig. 1(a) but also the outer space, do not change while changing the separation distance L between the colloids. Therefore, the effective inter-particle interaction, which is what we are interested in, is exclusively dictated by the contributions involving the surface interaction energy density ω γ,i (L) and the line interaction energy density ω τ (L). Please note that they are constructed such that ω γ,i (L → ∞) → 0 and ω τ (L → ∞) → 0.

B. Nonlinear theory
Within the nonlinear PB theory, the density functional in Eq. (1) is directly minimized with respect to n ± which provides the equilibrium profiles n eq ± for the ion number densities. Using these profiles, one finally ends up getting the nonlinear PB equation along with the same boundary conditions as listed below Eq. (2) to be satisfied by the electrostatic potential Ψ i (r) in medium i ∈ {1, 2}. In principle, after solving Eq. (4), the grand potential Ω (L) can also be obtained from Eq. (1) by inserting the equilibrium ions profiles n eq ± [Ψ]. However, in order to facilitate the whole process, we use the density functional instead and obtain the electrostatic potential Ψ (r) (defined as Ψ (r) = Ψ 1 (r) for x < 0 and Ψ (r) = Ψ 2 (r) for x > 0) by minimizing it numerically using Rayleigh-Ritz-like f inite element method [80]. Please note that the minimization of this functional with respect to Ψ (r) also results in the PB equation [Eq. (4)] together with all the boundary conditions that need to be satisfied. Once the potential Ψ (r), which minimizes Eq. (5) is obtained, it can be inserted back in Eq. (5) to obtain Ω min fem (L) which, at equilibrium, is related to the grand potential Ω (L) by Ω (L) = −Ω min fem (L). Without loss of generality, we first rewrite Eq. (5) by expanding the function cosh {βe (Ψ(r) − Ψ b (r))} in a Taylor series, where m can be interpreted as the degree of nonlinearity. For m = 1, one recovers the linearized PB problem whereas m → ∞ corresponds to the full nonlinear problem. We increase m in Eq.  3), with all the quantities in the right hand side carrying the same meaning as defined for Eq. (3). Numerical extraction of these contributions requires solving some subproblems with the following set-ups: (i) single fluid medium "1" ("2") spanning the lower (upper) half-space in Fig. 1(b) in the absence of any plates; the corresponding grand potential density Ω b,i = −2I i /β is easily obtained by simply setting Ψ(r) = Ψ b (r) and Ψ P = 0 in Eq. (6), (ii) only the two fluid media present in the absence of any plates, (iii) a single plate present in contact with a single semi-infinite fluid medium, (iv) two plates interacting across a single fluid medium, and (v) a single plate touching two semiinfinite fluids in its one side. A sequential subtraction of the grand potentials obtained for these subproblems then enables one to separate all the contributions in Eq. (7). For m = 1, the resulting interaction parameters are checked against those obtained analytically within the linear theory. With increasing m, for all the different system parameters considered here, we observe that Ω (m) (L) = −Ω (m)min fem (L) saturates for m ≥ 4. Therefore, in what follows, we call the results for m = 4 as the solutions under the nonlinear theory. This is in accordance with the conclusion of Ref. [51] and is not unexpected since we deal with the potentials of similar order of magnitude in both studies.

III. RESULTS AND DISCUSSION
A. Linear theory

Electrostatic potential
Our goal is to find the electrostatic potential which satisfies the DH equation (Eq. (2)) and the associated boundary conditions as mentioned below Eq. (2) for the system depicted in Fig. 1(b). To achieve this, we first dissect the original problem into the following subproblems: (i) two plates with constant surface potentials Ψ P interacting across fluid "1" or fluid "2" with bulk potential Ψ b,1 or Ψ b,2 , respectively; the resulting potential in both cases depend on z-coordinate only, (ii) two fluids phases occupying the two half spaces (x ≷ 0) forming an interface at x = 0 in the absence of any plates; the corresponding electrostatic potential, which is solely x-dependent, and the electric displacement vector are continuous at the interface. Exploiting the linear nature of the governing equation, i.e., the DH equation, if we simply add the solution for medium i in subproblem (i) with the corresponding ones obtained in subproblem (ii), the sums are also solutions of the DH equation and they satisfy most of the boundary conditions except for the continuity of the total electrostatic potentials at the interface. In order to get rid of this inconsistency we then seek for a correction function which also satisfies the DH equation and which, when added to the previously obtained sums, meets the remaining continuity condition without affecting any of the conditions that are already fulfilled. Such a function can easily be constructed by means of Fourier series expansion (for details, please see Chapter 3 of Ref. [78]). With all these, the final expressions for the electrostatic potentials in medium i ∈ {1, 2} read as where n 2 π 2 +κ 2 2 L 2 ε 1 p n,1 (L) + ε 2 p n,2 (L) and p n,i (L) = n 2 π 2 L 2 + κ 2 i . The first two terms in Eq. (8) together is the solution of the subproblem (i), i.e., potential distribution for two plates with constant potentials Ψ P interacting across medium i. The remainder represents the correction function. The solution of the subproblem (ii) is canceled by the first term (n = 0) of the series representing the correction function. From Eq. (8), one can trivially check that Ψ i (z = 0) = Ψ i (z = L) = Ψ P . In the limit x → ±∞, the series term vanishes because of the exponential function and one is left with the first two terms describing the potential due to two plates with constant surface potential Ψ P interacting across medium i with bulk potential Ψ b,i . This is exactly the situation far away from the interface. Upon approaching simultaneously the limits x → ±∞ and L → ∞, the series term vanishes and the second term reduces to an exponentially varying function characteristic of a single plate placed in contact with an electrolyte solution. In addition, if we set z → ∞ in Eq. (8), we recover the bulk potential Ψ b,i . As a side remark, we note that in the limit L → ∞ only, the second term reduces to an exponentially decaying potential and the sum can be converted to an integral over q = nπ/L. The resulting expression corresponds to the potential due to a single plate with surface potential Ψ P placed in contact with two immiscible fluids spanning the regions x ≷ 0 with z > 0. This can also be derived independently with the help of Fourier transforms (for details, please refer to the derivation of Eqs. (3.115) and (3.116) in Ref. [78]).

Interaction energies
With the electrostatic potential at hand, we can use it to calculate the grand potential Ω (L) = Ω [n ± [Ψ]] and to derive the expressions for the interaction parameters of our system. As defined in Eq. (3), this is achieved by distinguishing the terms proportional to V i , A i , A 1,2 , and . Moreover, the L-independent interactions are recognized in the terms proportional to A i and by taking the limit L → ∞ (for details, please see Chapter 4 of Ref. [78]).
The surface tensions, as defined in Eq. (3), acting between each surface and its adjacent fluid medium i ∈ {1, 2} is given by and the surface interaction energy between two surfaces with total area A i in medium i ∈ {1, 2} is given by Whereas the amplitudes of both γ i and ω γ,i (L) depend on ε i , κ i , Ψ P , and Ψ b,i , the decay rate of ω γ,i (L) is uniquely governed by κ i . In the limit of vanishing separation between the plates, Eq. (10) predicts a finite, repulsive interaction: Contrary to a constant charge boundary condition, this non-divergent behavior in the L → 0 limit is characteristic of a constant potential boundary condition [56]. In the large asymptotic limit, Eq. (10) provides . Therefore, as dictated by the leading order term, the surface interaction ω γ,i (L) decays monotonically ∼ e −κ i L for large separations between the plates.
The interfacial tension per total interfacial area A 1,2 acting between the two fluid media reads as As expected, the interfacial tension depends only on the properties of the two fluid media and not on any plate properties. Consequently, the expression for γ 1,2 does not differ from what one obtains for plates with constant surface charge densities [53]. The line tension acting at both the three-phase contact lines expressed per total length of the two contact lines is given by and p i (q) = q 2 + κ 2 i , and the line interaction energy (defined in Eq. (3)) acting between the two three-phase contact lines is given by with ξ n (L) and p n,i (L) as defined below Eq. (8) and the expression for τ given in Eq. (12). In the limit of vanishing separation (L → 0) between the plates, all the terms in Eq. (13) vanish except for the last term, i.e., τ. Therefore, in this limit ω τ (L) stays finite unlike what one observes for constant charge boundary condition [50,53]. Using the relation ∞ n=1,3,...
it is trivial to see that the terms ∝ L, i.e., the third and the fourth terms in Eq. (13) cancel each other in the large asymptotic limit (L → ∞). The rest of the terms together lead to an overall exponential decay of the line interaction ω τ (L) in this limit. A striking observation regarding the line interaction for constant surface charge density boundary condition is that it is independent of the Donnan potential Ψ D [50,53]. On the contrary, as Eq. (13) suggests, ω τ (L) does depend on Ψ D while using a constant surface potential boundary condition. This is indeed reasonable since the line part essentially captures the effects of the interface to the total inter-plate interaction and the contrast in the bulk potential should get reflected in it.

B. Nonlinear theory
In this subsection, we discuss the results within the nonlinear theory and compare them with those obtained using the linear theory in Sec. III A. In order to present our numerically obtained results and to gain insight into the variations of different parameters efficiently, we use dimensionless quantities. Accordingly, ω γ,i (L) and ω τ (L) are expressed as βω γ,i (L)/κ 2 1 and βω τ (L)/κ 1 , respectively and their behavior is studied as functions of the scaled separation κ 1 L. In doing so, the number of dimensionless free parameters reduces to only four: βeΨ P , βeΨ D , I = I 2 /I 1 , and ε = ε 2 /ε 1 . As the permittivities of oils used in typical experimental systems widely vary [19,31,58], we consider two experimental systems differing significantly in their fluid contents for our discussion. The variations of different system parameters with respect to one of these systems are also presented which allow one to infer what could be expected to happen for an arbitrary general system.

Water-lutidine system
As our first example, we take a system with particles at a water-lutidine (2,6-dimethylpyridine) interface at temperature T = 313 K above the critical temperature [81]. The waterrich phase is denoted as the medium "1" and the lutidine-rich phase as the medium "2". The dissolved salt is NaI with bulk ionic concentrations I 1 = 1 mM and I 2 = 0.85 mM which leads to a Donnan potential estimated to be 1 k B T/e [77,82]. The relative permittivities of the two fluids are ε r,1 = 72 and ε r,2 = 62 [83,84]. Typically metal surfaces are negatively charged in an electrolyte solution [85]. Therefore, the particles are considered to be carrying a constant potential Ψ P = −3 k B T/e ≈ −77 mV. Please note that the surface potential is expected to be slightly higher than the zeta potentials usually measured in experiments. Therefore, the dimensionless parameters of our system, which we call as the standard system, are I = 0.85, ε = 62/72, βeΨ D = 1, and βeΨ P = −3.
The ensuing interaction energies between the two plates are presented in Fig. 2. For all the plots, the data points correspond to the numerically obtained results and the solid lines represent the analytically obtained expressions in Subsection III A 2. Figure 2(a) shows the variation of the scaled surface interaction energy density ω γ,1 (L) acting between the portions of the plates in contact with medium "1" as function of the varying scaled separation L. As one can see, both the linear and the nonlinear theory predict a monotonically decaying repulsive interaction between the plates, and as expected, the numerically obtained results corresponding to m = 1 in Eq. (7) match perfectly with the analytically obtained solution within the linear theory throughout the range of separations considered here. While the results within the two theories qualitatively agree with each Variations of the surface interaction energy densities ω γ,1 (L) and ω γ,2 (L), both expressed in units of κ 2 1 /β and (c) variation of the line interaction energy density ω τ (L) expressed in units of κ 1 /β as functions of the scaled separation κ 1 L between the two surfaces for standard set of parameters (βeΨ P = −3, βeΨ D = 1, ε = 62/72, and I = 0.85) corresponding to particles trapped at a water-lutidine interface. As shown by all the plots, the numerically obtained results for m = 1 in Eq. (7) match perfectly with the analytically obtained results within the linearized PB theory. Whereas both ω γ,1 (L) and ω γ,2 (L) are repulsive and decay monotonically with increasing separation distance κ 1 L under linearized as well as nonlinear PB theory, the linear theory underestimates them at very short separations and overestimates at larger separations. On the other hand, the linear theory fails to qualitatively predict the correct variation for ω τ (L). Whereas ω τ (L) is attractive everywhere in the linear case, it can be attractive as well as repulsive within the nonlinear theory.
other, the linear theory underestimates the interaction at very short separations and overestimates at relatively larger separations. This discrepancy gradually drops down as the separation distance is increased. Similar conclusions can be drawn from Fig. 2(b) concerning the variation of the scaled surface interac- The parameters (βeΨ P = −3, βeΨ D = 3.8, ε = 10.3/80, and I = 2.9 × 10 −4 ) used for the plots correspond to a system consisting of particles trapped at a water-octanol interface. As shown by all the plots, whereas both ω γ,1 (L) and ω γ,2 (L) are repulsive and decay monotonically with increasing separation distance κ 1 L under linearized as well as nonlinear PB theory, the linear theory underestimates them at very short separations and overestimates at larger separations. Regarding the variation of the line interaction ω τ (L), whereas the linear theory qualitatively predicts the correct behavior, it is inaccurate in predicting the magnitude of ω τ (L) and the location of the minimum. tion energy density ω γ,2 (L) acting between the portions of the plates in contact with medium "2" as function of the varying scaled separation L within the two theories. The magnitudes of ω γ,1 (L) and ω γ,2 (L) are of the same order as the two fluids do not significantly differ in their properties and the surface potentials of the plates are also the same in both media. However, the surface interaction in medium "2" is slightly stronger com-pared to medium "1" due to less screening. Figure 2(c) shows the variation of the scaled line interaction energy ω τ (L) as function of the scaled separation L within both the theories. As one can infer, within the range of separations considered, the linear theory predicts a monotonically decaying attractive interaction with increasing separation distance between the plates. However, within the nonlinear theory the line interaction shows a qualitatively different behavior. It changes non-monotonically upon increasing the separation distance, shows a maximum at κ 1 L ≈ 1 followed by a minimum at around κ 1 L ≈ 2 and then decays monotonically. Clearly, the interaction is repulsive in between the maximum and the minimum. Please note that both the maximum and the minimum occur at distances much larger than the molecular length scale (typically a few angstroms). Not only that, up to around κ 1 L ≈ 2, the linear theory significantly overestimates the strength of the line interaction as well.

Water-octanol system
Water-lutidine falls in a class of systems where the bulk properties, i.e., the relative permittivities and the bulk ionic strengths, and consequently the Debye screening lengths vary little in the two fluid phases. Although such combination of immiscible fluids are used in experiments [19,85], there is another frequently used category of systems where moderately polarizable oils are used [19]. Therefore, as a second example, we consider a system consisting of particles trapped at a water-octanol interface. At room temperature T = 300 K these two fluids are characterized by significantly different relative permittivities with ε r,1 = 80 (for water as medium "1") and ε r,2 = 10.3 (for octanol as medium "2"). The ion-partitioning at this interface results in highly contrasting bulk ionic strengths as well. For I 1 = 10 mM in the water phase, I 2 = 2.9 × 10 −3 mM in the octanol phase, leading to a Donnan potential Ψ D ≈ 3.8 k B T/e. The associated inverse Debye lengths in the two phases are given by κ 1 ≈ 0.324 nm −1 and κ 2 ≈ 0.015 nm −1 . The particles are considered to be carrying the same surface potentials as before, i.e., Ψ P = −3 k B T/e ≈ −77 mV. Hence, the dimensionless parameters of this system are given by I = 2.9 × 10 −4 , ε = 10.3/80, βeΨ D = 3.8, and βeΨ P = −3.
The resulting interactions for this system are shown in Fig. 3. Both ω γ,1 (L) and ω γ,2 (L), as shown in Figs. 3(a) and 3(b), respectively predict monotonically decaying repulsive interactions under the linear as well as the nonlinear theory. But, as before, the linear theory underestimates the interactions at very short separations whereas overestimates them at relatively larger separations. Depending upon the fluid medium and the separation distance L, this mismatch can be significant (see, for example, the data presented in Fig. 3(b)). Except for very short separations, the interaction in medium "2", i.e., ω γ,2 (L) is stronger compared to that in medium "1" (presented by ω γ,1 (L)). This is due to a weaker screening in the medium "2" compared with medium "1" which is evident from the rapid decay of ω γ,1 (L) in Fig. 3(a). At very short separations, however, one enters into the electric double layer and with decreasing As one can see, the interaction energy decreases initially, shows a minimum, and finally increases monotonically within both the linear (red crosses) and the nonlinear PB theories (blue open circles). The initial decreases in the interaction energy corresponds to a positive effective force at short separations L, implying that the interaction is repulsive there. Beyond the position of the minimum, it turns into an attractive one. Whereas this feature remains the same within the two theories, the location of the minimum (as indicated by the arrows in the inset) shifts from L ≈ 48 nm to a significantly shorter separation L ≈ 38 nm while calculated using the nonlinear PB theory. Therefore, the particles will come much closer than what is predicted by the linear theory. separations, the ions, which are strongly attracted to the surfaces, need to be removed. As the ionic strength is higher in medium "1", the system needs to remove more ions leading to a stronger repulsion compared to medium "2". Concerning the line interaction energy ω τ (L), as shown in Fig. 3(c), it decays non-monotonically with increasing separation L between the two plates. Initially its magnitude increases implying a repulsive interaction, which, after the occurrence of a minimum, turns into an attractive interaction. Although the linear theory predicts the qualitative behavior correctly, it significantly overestimates the strength of the interaction over the entire range of separations considered here and also fails to accurately locate the position of the minimum. It is worth mentioning that the line interaction for water-octanol system is much stronger compared to that for the water-lutidine system owing to greater contrast between the combination of the fluids used. While comparing the magnitudes in Figs. 2 and 3, please note that κ 1 ≈ 0.324 nm −1 for the water-octanol system whereas κ 1 ≈ 0.106 nm −1 for the water-lutidine system. So far we have presented results concerning all the different parts, i.e., ω γ,1 (L), ω γ,2 (L) and ω τ (L) to the total electrostatic interaction. This indeed allows one to investigate the system in a detailed way and gain better insight from a the- Variations of the surface interaction energy density ω γ,2 (L) expressed in units of κ 2 1 /β (left panels) and variations of the line interaction energy density ω τ (L) expressed in units of κ 1 /β (right panels) as functions of the scaled separation κ 1 L between the two plates for various sets of free dimensionless parameters βeΨ P (panels (a) and (b)), βeΨ D (panels (c) and (d)), ε (panels (e) and (f)), and I (panels (g) and (h)). In all cases, only a single parameter has been changed at a time from the standard water-lutidine system and the interactions are compared. The parameters that are changed are mentioned in each plot. As shown by all the plots in the left panel, ω γ,2 (L) is always repulsive and it increases with increasing difference βe |Ψ P − Ψ D | or βeΨ D or ε. On the other hand, whereas an increase in I strengthens ω γ,2 (L) at very short separations, it weakens ω γ,2 (L) at relatively larger separations. As shown by panel (b), ω τ (L) shows a maximum followed by a minimum for strongly or moderately negative surface potentials. On the other hand, for positively charged or weakly negatively charged particles ω τ (L) remains attractive everywhere. With increasing values of βeΨ D , ω τ (L) changes from overall attractive to a non-monotonically varying behavior showing a maximum at very short separations followed by a minimum at larger separations. For βeΨ D = 0, ω τ (L) almost vanishes as the contrast between the two fluid media reduces in this case. With decreasing ε, the height of the maximum of ω τ (L) increases and it shifts to the left whereas the minimum becomes shallow and shifts to larger distances. On the other hand, with decreasing I, whereas the position of the maximum does not change, the minimum shifts to larger distances and the magnitude of ω τ (L) increases. Insets of panels (g) and (h) show the variations at larger separation distances as the Debye screening length is larger for the smaller I values.
oretical point of view. We are also not required to restrict ourself to fixed particle sizes. However, experimentally it might be challenging to disentangle all these individual interaction parameters. Also, one might reasonably wonder about the relative importance of the line interaction energy ω τ (L) in comparison with the surface interaction energies ω γ,1 (L) and ω γ,2 (L). As Eq. (3) or (7) suggests, the surface contributions to the total grand potential are proportional to the surface areas of the plates whereas the line contribution scales with the length of the three-phase contact lines. Therefore, for large surfaces, the surface contributions dominate anyway in due course. But what happens for typical system parmeters? To address these concerns, we present in Fig. 4 the total interaction energy ∆Ω(L) = ω γ,1 (L)A 1 + ω γ,2 (L)A 2 + ω τ (L) (i.e., the total L-dependent part in Eq. (3) or (7)) for a wateroctanol system with βeΨ P = −3, βeΨ D = 3.8, ε = 10.3/80, I 2 = 2.9 × 10 −2 mM, I 1 = 0.1 M, A 2 = A 1 = 2573 nm 2 , and = 114 nm. The effective areas A i and and effective length of the three-phase contact lines are rough estimates for 60 nm particles. The quantity, i.e., ∆Ω(L) we plot in Fig. 4 is related to the effective electrostatic force F(L) = −∂ (∆Ω) /∂L which can be considered as experimentally more relevant quantity. As the plot suggests, the total interaction energy ∆Ω(L) initially decreases with increasing separation, reaches a minimum (marked with the arrows in the inset of Fig 4) and then increases within both the linear and the nonlinear PB theories. This implies a repulsive interaction till the position of the minimum and an attractive interaction beyond that. As both the surface parts are repulsive in Fig. 3 throughout the range of separations κ 1 L considered here, the attractions clearly come from the line parts. Although both the linear and the nonlinear theory predict an attractive interaction beyond certain separations, the equilibrium separation (given by the position of the minimum) decreases significantly from L ≈ 48 nm to L ≈ 38 nm while calculated using the nonlinear PB thoery. It is worth mentioning that the line interaction starts to be comparable in magnitude with the surface parts at much shorter distances. Therefore, depending on the system, the line contribution can indeed dominate over the surface contributions even for typical system parameters. We also note that, as discussed in Subsection III B 3 and as can be seen in Figs. 5(a) and 5(b), the minimum can be moved to smaller separations and its depth can be increased by reducing the surface potential βeΨ P as the magnitudes of the surface parts reduce whereas that of the line part increases with decreasing surface potential.

Variation of system parameters
Now we discuss the effects of changing the free parameters, i.e., βeΨ P , βeΨ D , ε, and I, of our system within the nonlinear PB theory. The comparisons are done with respect to the standard (water-lutidine) system. Please note that, the permittivity ratio ε = ε 2 /ε 1 and the ionic concentration ratio I = I 2 /I 1 can vary due to changes in the respective quantities in either or both fluids. However, the resulting total interaction is sensitive only to the ratios ε and I. For our analysis, we have chosen to vary ε and I by changing ε 2 and I 2 , respectively while keeping ε 1 and I 1 fixed. Therefore, only the surface interaction ω γ,2 (L) in medium "2" and ω τ (L) changes due to such variations but ω γ,1 (L) remains unaltered. ω γ,1 (L) is also independent of the Donnan potential βeΨ D . However, it changes due to variation in βeΨ P but this happens in a fashion similar to what one observes for ω γ,2 (L). Therefore, here we show only the variations of ω γ,2 (L) and ω τ (L) upon changing the system parameters.
As one can see from Fig. 5(a), the surface interaction ω γ,2 (L) in medium "2" increases with absolute value of the increasing contrast between surface potential Ψ P and the bulk potential Ψ b,2 . Consequently, no change is observed while varying βeΨ P from 3 to −1 (please note that βeΨ b,2 = 1 here). Similarly, ω γ,2 (L) increases with increasing Donnan potential βeΨ D and the permittivity ratio ε; see Figs. 5(c) and 5(e). Concerning the variation of ω γ,2 (L) upon increasing the ionic strength ratio I, ω γ,2 (L) increases at very short separations whereas it decreases at relatively larger separations; see Fig. 5(g). At these larger separations, with increasing ionic concentration, electrostatic screening increases. Therefore, the effective electrostatic interaction diminishes. However, at very short separations, within the electric double layer, the number of ions strongly attracted to the surfaces is more for higher ionic strengths. Consequently, more work needs to be done in order to decrease the separation distance as it requires removing these ions. This results in higher values for ω γ,2 (L). Please note that irrespective of the specific values of the system parameters ω γ,2 (L) is always repulsive.
The variations of ω τ (L) with different system parameters are shown in the right panel of Fig. 5. As one can see from Fig. 5(b), ω τ (L) shows a maximum followed by a minimum for strongly or moderately negative surface potentials. Therefore, the line interaction can be repulsive as well as attractive depending upon the separation distance L. On the other hand, for positively charged or weakly negatively charged particles the line interaction ω τ (L) remains attractive everywhere and decays monotonically. With increasing values of βeΨ D , ω τ (L) changes from overall attractive monotonically decreasing behavior to a non-monotonically varying behavior showing a maximum at very short separations followed by a minimum at larger separations; see Fig. 5(d). A special situation is presented by βeΨ D = 0 which leads to an almost vanishing ω τ (L) as the contrast between the two fluid media reduces in this case. With decreasing ε, the height of the maximum of ω τ (L) increases and it shifts to the left whereas the minimum becomes shallow and shifts to larger distances; see It is important to note that, the equilibrium separation between the particles is determined by the total interaction energy ∆Ω(L) = ω γ,1 (L)A 1 + ω γ,2 (L)A 2 + ω τ (L) ; the curves in Fig. 5 do not directly dictate it. However, they provide important informations regarding what one should expect for a given experimental set-up. For example, the surface interactions being always repulsive, inter-particle attraction can only come from the line part. As the plots suggest, the line part can either be monotonically attractive or it can change from repulsive to attractive beyond the position of a minimum at short separations. I. Values of the surface tensions γ 1 (third column) and γ 2 (fourth column) resulting from interactions of the plates across medium "1" and medium "2", respectively, the interfacial tension γ 1,2 (fifth column) acting between the two fluids, and the line tension τ (last column) acting at the three-phase contact lines for all the systems considered in this study (listed in the first column). The first block represents the results for our standard water-lutidine system and the last block provides results corresponding to the water-octanol system. The intermediate blocks provide results for system where a single parameter has been changed at a time from the standard system. These systems are denoted by the parameter that is changed. The surface and the interfacial tension values are given in the units of k B T/nm 2 while the line tension values are provided in the units of k B T/nm. For each system, the results obtained analytically within linearized PB theory and numerically within nonlinear as well as linear PB theory are presented. As expected, the numerically obtained values for m = 1 in Eq. (7) (denoted as 'Linear') agree with those obtained analytically. As it is evident from the presented values, γ 1 , γ 2 , and γ 1,2 are always negative for the systems considered here and the linear theory always underestimates their absolute values. However, the line tension τ appears to be the most sensitive to the variations in the system parameters and to the type of theory used to calculate its value. It can be positive as well as negative, decrease as well as increase and even change sign while calculating within the linear and the nonlinear theories. Even if it is monotonically attractive, at very short separations (below L ≈ 20 nm) the surface parts usually dominate. Only in the case of small ε values, the surface part is strongly suppressed. However, the line part itself becomes repulsive in this case (see Fig. 5(f)). The curves plotted in Figs. 2, 3 and 5 also suggest that if a minimum exists at a non-vanishing separation, it typically occurs between L ≈ 20 nm−40 nm (please note that κ 1 for Fig. 3 is different than that in Figs. 2 and 5). Therefore, in view of these results, we conclude that overall electrostatic attraction can start to kick in at or above L ≈ 20 nm and the particles will arrange themselves with an inter-particle separa-tion L 20 nm. It is worth mentioning that the relevance of the line part and the precise location of the minimum not only depend on the system parameters considered in Fig. 5 but also on the surface areas exposed to the fluids and on the length of the three-phase contact line. However, from Fig. 5, one can easily say that the magnitude of the line part can be strongly enhanced by increasing the contrast in the ionic strengths in the two media (or equivalently, decreasing the ionic strength ratio I = I 2 /I 1 ). On the other hand, for a given separation, the surface parts can be suppressed by decreasing the relative permittivity value, by reducing the surface potential, or by in-creasing the electrostatic screening. In fact, this is also evident if one compares the two experimental systems (water-lutidine and water-octanol) discussed previously. What remain to be discussed are the L-independent interactions in Eqs. (3) or (7). Their values are listed in Table I for all the different system parameters considered here. As one can see, in all cases the numerically obtained values for m = 1 in Eq. (7) (denoted as 'Linear') agrees perfectly with the analytically obtained results (denoted as 'Analytical') using Eqs. (9), (11), and (12). As one can see from Table I, the surface tension γ 1 acting between the plates and medium "1" is always negative. However, its absolute value is higher within the nonlinear theory compared to the linear theory. Similar observations hold true for the surface tension γ 2 acting between the plates and medium "2". Both of them also increases with the absolute value of the difference between the surface potential βeΨ P and the bulk potential βeΨ b,i . As one can infer from the fourth column, the absolute value of γ 2 decreases with decreasing βeΨ D , ε, and I within both the nonlinear and the linear theories. As expected, the interfacial tension γ 1,2 (given in the fifth column) remains unaffected due to changes in the surface potential Ψ P and is insensitive to the sign of the Donnan potential βeΨ D . With decreasing ε and I, the magnitude of γ 1,2 also decreases. For all the situations considered here, γ 1,2 is always negative except for βeΨ D = 0 when it vanishes completely. The quantity that is most sensitive to changes in the system parameters is the line tension τ presented in the last column. It can be positive as well as negative. Not only that, it can decrease as well as increase and even change sign while calculating within the linear and the nonlinear theories. The values for the line tensions that we obtain are within the range 0.1 pN − 1 pN, which is in accordance with the values reported in the literature [86].
We note that, although not shown here, the agreement of the results corresponding to m = 1 in Eq. (7) with those obtained analytically within the linear theory have been checked for all the systems considered here. As another check, we also observe that with decreasing βe Ψ P − Ψ b,i values, the difference between the results within the linear and the nonlinear theory diminishes.

IV. CONCLUSIONS
To summarise, by using the classical DFT formalism, we have addressed the problem of the electrostatic interaction between a pair of identical particles with constant surface potentials sitting in close proximity of each other at a fluid-fluid interface. The particles could be either monometallic, or metal coated such as Janus particles with metallic caps facing each other or core-shell particles with metallic shells. Considering a simple yet reasonable model system, which exploits the short separation between the particles, we have solved the problem within both the linear and the nonlinear PB theory. Within the linear theory, closed-form analytical solutions are obtained for all the interaction parameters present in the system. For the nonlinear PB theory, the problem is solved numerically. While the results within the linear theory are mostly quantitatively inaccurate, they can be qualitatively reliable depending upon the system. The surface interaction energies between the two particles are always found to be repulsive and monotonically decaying. However, the line interaction energy varies nonmonotonically with increasing separation between the particles and it can be attractive as well as repulsive. Not only that, contrary to common expectation, it is shown that the attractive line interaction can be strong enough to dominate over the sum of the surface parts and to dictate the nature of the overall electrostatic interaction for typical experimental systems. A comparison between two typical experimental setups illustrates that the importance of the line interaction enhances with increasing contrast in the bulk properties of the two fluids forming the interface. Therefore, special attention should be given to line interaction in the vast majority of experiments performed using liquids varying starkly in their bulk properties. We also provide results concerning the interactions such as surface, line, and interfacial tensions that do not depend on the separation distance between the particles. Finally, we note that the presented theoretical framework is based only on a simple minimization of the system free energy and can, in principle, be easily be applied to related interfacial problems dealing not necessarily with planar surfaces or interfaces. Therefore, on one hand, we expect our study to directly contribute towards improving the understanding and modeling of interfacial monolayers of colloidal particles with constant surface potentials. On the other hand, we expect it to serve as a useful reference for future theoretical investigations of related problems.