Self-assembly of repulsive interfacial particles via collective sinking

Duck-Gyu Lee a, Pietro Cicuta b and Dominic Vella *a
aMathematical Institute, Radcliffe Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, UK. E-mail:; Tel: +44 (0)1865 615150
bBSS, Cavendish Laboratory, Cambridge, CB3 0HE, UK

Received 15th April 2016 , Accepted 23rd June 2016

First published on 23rd June 2016

Charged colloidal particles trapped at an air–water interface are well known to form an ordered crystal, stabilized by a long ranged repulsion; the details of this repulsion remain something of a mystery, but all experiments performed to date have confirmed a dipolar-repulsion, at least at dilute concentrations. More complex arrangements are often observed, especially at higher concentration, and these seem to be incompatible with a purely repulsive potential. In addition to electrostatic repulsion, interfacial particles may also interact via deformation of the surface: so-called capillary effects. Pair-wise capillary interactions are well understood, and are known to be too small (for these colloidal particles) to overcome thermal effects. Here we show that collective effects may significantly modify the simple pair-wise interactions and become important at higher density, though we remain well below close packing throughout. In particular, we show that the interaction of many interfacial particles can cause much larger interfacial deformations than do isolated particles, and show that the energy of interaction per particle due to this “collective sinking” grows as the number of interacting particles grows. Though some of the parameters in our simple model are unknown, the scaling behaviour is entirely consistent with experimental data, strongly indicating that estimating interaction energy based solely on pair-wise potentials may be too simplistic for surface particle layers.

1 Introduction

Since first being observed by Pieranski in 1980,1 the self-assembly of colloidal particles at a liquid–fluid interface has sparked considerable interest as both a system in which to study the nature of phase transitions in two-dimensions2,3 and as a useful tool for designing regular arrays at small scales.4 A particularly striking feature of these self-assembled colloidal crystals is that the constituent colloids maintain an equilibrium separation that can be several times larger than the particle diameter:5 this feature is particularly useful in particle lithography since it allows the fabrication of patterns without particle–particle contact, which can cause cross-talk in various applications.4 At a fundamental level, this separation clearly indicates that the particles are repulsive; detailed studies show that the charge on the particles is, in fact, separated so that the colloids behave as electrical dipoles6,7 or even quadrupoles.8,9

Pieranski1 believed the interaction between particles to be purely repulsive: he attributed the regular spacing of colloids to be due to the geometrical confinement of the system. However, isolated clusters of particles have been observed experimentally,10–13 suggesting that this geometrical confinement is not necessary for the formation of a well-ordered crystal. To form such a bound state requires a long-ranged attractive force to overcome the electrostatic repulsion. The existence of such an attraction was also suggested from the inversion of the pair-correlation function14 and could be at the origin of stable clusters studied by Nikolaides et al.15 However, the question of what provides the attractive interaction that balances this repulsion remains open.

At first sight, it is natural to assume that this attraction could result from the well-known attraction between floating particles that is mediated by meniscus deformation,16–18 sometimes called the ‘Cheerios effect’.19 However, it is also well-known that for pairwise interactions, the interaction energy of these flotation forces17 becomes on the order of the thermal energy, kBT, for particles of radius smaller than around 10 μm,17 which is precisely the scale at which these colloidal crystals are observed. Nevertheless, an attractive force persists and so the question remains: what is the basic mechanism behind the attractive force?

Several different mechanisms have been proposed as the origin of the attractive force, including undulations due to a rough contact line13 and/or enhanced normal forces of electrical origins9,15,20 – an electro-dipping force. However, none of these explanations are able to satisfactorily explain all experimental observations: particle roughnesses of the size suggested by Stamou et al.13 were not observed experimentally21 while observations with imposed electrical fields to vary the strength of the electro-dipping force did not produce the expected variation in particle spacing.7 We are not able to resolve these disagreements here, or to propose a new electrostatic mechanism. Instead, we revisit the tacit assumption that the interaction energy between a pair of particles is a useful way of estimating the typical interaction energy for a large number of interfacial particles.

The qualitative change to the flotation of particles caused by other nearby objects is now well documented.22–24 At macroscopic scales, rafts of dense objects float significantly deeper in the liquid than they do in isolation. This is because the proximity of other particles in the raft constrains the menisci to be more horizontal than they would be for an isolated particle: the particles thus sink deeper into the liquid so that hydrostatic pressure can make up for the loss of supporting force from surface tension. Indeed, this effect can be so dramatic that dense particles that are able to float in isolation may actually sink in the proximity of enough other floating particles.22,23 While the dramatic loss of floating stability is unlikely to be relevant at the very small scales of colloidal particles, the observation that their vertical force balance may be affected by the presence of other particles is likely to be robust. The question we address in the remainder of this paper is how any alteration to the vertical force balance manifests itself in the horizontal force balance condition – is the effective interaction energy between particles substantially modified from the two-body case? This question has been addressed previously using mean-field, or coarse-grained approaches;25–31 here, we study this problem by considering in detail the meniscus deformations caused by individual particles and the collective effect of this deformation.

We begin by considering a two-dimensional model problem in Section 2, which allows us to obtain numerical results for large assemblies of particles. These results can be understood quantitatively using a scaling analysis, which is then extended to the three-dimensional problem of most interest in Section 4. We then conclude in Section 5 by discussing the possible significance of our theoretical results for the spontaneous formation of colloidal islands.

2 Two-dimensional formulation

Here we consider a two-dimensional configuration in detail. With this simplification, we are then able to consider a point-like particle (corresponding to a line in 3D), facilitating our analysis. We aim to gain physical understanding here, that can then be used to inform an understanding of the fully three-dimensional problem.

2.1 Key ingredients of the model

We are interested in understanding in detail the flotation of a series of objects; to be able to form the crystals that are observed, we expect that these objects should be subject to both attractive and repulsive interactions. The question of interest is then really can an ordinarily small pairwise interaction be amplified in a many-body interaction? The simplest purely repulsive interaction in 2D is that between line charges (rather than dipoles, which can become attractive depending on their orientation). By ‘line charge’ here, we mean the two-dimensional analogue of a point charge: charge exists along a line perpendicular to the plane of Fig. 1. We therefore consider a collection of N identical line charges floating at the interface between a liquid, of density ρl, and a fluid, of density ρf < ρl; the interfacial tension is γ. Each line charge has an electric charge +q per unit length (so that electrostatic interactions are purely repulsive) as well as a weight per unit length mg. A sketch of the setup is shown in Fig. 1.
image file: c6sm00901h-f1.tif
Fig. 1 Geometry of floating line charges at the interface between a liquid and a fluid. (a) The general formulation of the problem for many charges, showing the geometrical and force balance aspects as insets. (b) The specialization to the two-body problem. Note that, in each case, the angles of the interface to the horizontal are defined positive above the horizontal.

In this simple model, the point-like particles deform the interface purely due to their weight: there are no wetting effects to be considered. We envisage that this weight-induced deformation will be small since the particles of interest are themselves small and easily supported by surface tension. As a consequence, the attractive interaction due to the weight-induced meniscus deformation between a pair of these particles should also be small, leading to a relatively large equilibrium separation at which capillary attraction balances electrostatic repulsion. Equivalently, we expect the typical interaction energy between a pair of such particles (N = 2) to be very small (compared to the thermal energy). To understand the balance between deformation-induced attraction and electrostatic repulsion, we first consider the two-body problem in some detail.

2.2 The two-body problem

The interaction between two identical line charges will clarify the various dimensionless parameters that control the problem. We orientate our axes so that the two particles are located at x = ±d/2 with d the distance that separates them (see Fig. 1b). The meniscus y = h(x) must satisfy the Laplace–Young equation
image file: c6sm00901h-t1.tif(1)
with primes denoting differentiation with respect to x. Eqn (1) is to be solved subject to a symmetry condition, h′(0) = 0, and the meniscus decay condition, h(±∞) = 0. However, the meniscus may have a discontinuity in slope at each particle (indeed this discontinuity is what leads to a horizontal force between particles in this line-particle limit).
2.2.1 The linearized problem. Under the assumption that the particles only deform the interface slightly, so that the slope of interface deformations h′ ≪ 1, the Laplace–Young eqn (1) may be linearized. The deformation of the interface caused by a single line charge at x = xi may then be found analytically to be
image file: c6sm00901h-t2.tif(2)
where x is the horizontal coordinate measured from the mass,
image file: c6sm00901h-t3.tif(3)
is the capillary length and
W = mg/γ(4)
is the dimensionless weight per unit length of the mass. Here, the prefactor is determined by the vertical force balance on the mass – the vertical force provided from surface tension must balance the weight of the line mass.

If two identical masses float with separation d then by linear superposition,32 we have

image file: c6sm00901h-t4.tif(5)
To determine the separation distance at equilibrium, d, we use the horizontal balance between the (repulsive) electrostatic force and the (attractive) capillary forces, which reads
image file: c6sm00901h-t5.tif(6)
where ε0 is the permittivity of free space and the angles β± are the interfacial inclinations at the contact point, given in terms of the meniscus profile by tan[thin space (1/6-em)]β± = ±h′(d/2). Using (5) to determine the leading order behaviour of cos[thin space (1/6-em)]β+ − cos[thin space (1/6-em)]β for W ≪ 1 we find that the equilibrium separation, d, is the solution of the equation
image file: c6sm00901h-t6.tif(7)
image file: c6sm00901h-t7.tif(8)
is the dimensionless charge parameter, which measures the strength of electrostatic repulsion at separation d = [small script l]c in comparison to the typical force from surface tension.

A sketch of the RHS in (7) (see Fig. 1b) reveals that it is a non-monotonic function of d/[small script l]c; in particular, an equilibrium is only possible for sufficiently weak repulsion, or sufficiently large weight, that C2/W2 ≤ 0.184. For a given value of C2/W2 ≤ 0.184 there are two equilibria, the smallest of which is stable and the largest of which is unstable. For d/[small script l]c ≪ 1, the position of the stable equilibrium is given by

image file: c6sm00901h-t8.tif(9)
The analytical progress allowed by the assumption of a small, linear, deformation yields key insight. In particular, we see that as C2/W2 decreases, so does the equilibrium separation between them, d. As a result of this, the depth at which each particle floats
image file: c6sm00901h-t9.tif(10)
increases as C2/W2 decreases. This sinking is caused by the presence of nearby objects and so we refer to it as ‘collective sinking’ here. The fact that the presence of nearby objects modifies the vertical force balance and hence can cause objects that would float in isolation to sink, has been observed at macroscopic scales previously.22–24,33 Here, we do not consider this sinking transition, but emphasize the key point that the presence of a second particle nearby, via its interfacial deformation, modifies the behaviour of a first particle. This is similar to the capillary collapse studied recently.27–31 We shall shortly go beyond the mean-field approach adopted in these previous works by explicitly considering ensembles of particles accounting for the interface deformation beyond the linear theory just presented. To see the possible effect of the nonlinearities, we first reconsider the two-body problem, accounting for nonlinear meniscus deformation.

2.2.2 The nonlinear problem. The above analysis hinged on the assumption that the slope of the interface, h′ ≪ 1, so that the Laplace–Young eqn (1) could be linearized and the meniscus deformations caused by each particle superposed. While this is a valid assumption for large particle separations and small particle weights, we now investigate how the results of the previous sub-section change once nonlinear interface deformations are considered.

When the meniscus slope is no longer considered to be small, the Laplace–Young equation must be solved numerically. In fact, all that is required is to find the angles that the menisci make with the horizontal, β± in Fig. 1b. This simplification, and the fact that the external menisci extend to infinity, mean that we may make use of well known34 first integrals of the Laplace–Young eqn (1) which give

image file: c6sm00901h-t10.tif(11)
image file: c6sm00901h-t11.tif(12)
To determine the angle β, however, we must resort to a numerical solution of (1) subject to the boundary conditions
image file: c6sm00901h-t12.tif(13)
(These relations express symmetry and vertical force balance, respectively.) Once β and h* have been determined for a particular configuration, the horizontal force balance (6) gives the dimensionless charge C required for flotation at that equilibrium separation. The results of this numerical calculation, and a comparison with the corresponding result for small deformations (7), are shown in Fig. 2a. We observe that the trend is very similar to that observed in the linear theory, although the nonlinear equilibrium separation at fixed C2/W2 decreases as W increases: the nonlinear effect of nearby particles is to draw those particles closer together than would be supposed from the linear theory.

image file: c6sm00901h-f2.tif
Fig. 2 Comparison between the linearized and fully nonlinear approaches to the two body problem. (a) The equilibrium separation deqm as a function of the charge-to-weight ratio C2/W2. The result of the linear analysis (7) (solid black curve) is shown together with the result of full nonlinear computations for W = 0.2 (red dotted), W = 0.4 (orange dashed), W = 0.6 (green dash-dotted) and W = 0.8 (blue dash-double dotted). The black dashed line represents the asymptotic result (9), valid for C2/W2 ≪ 1. (b) The combined energy of the system (compared to that of a flat interface and infinitely separated charges) at the corresponding equilibrium separation, deqm. Here the result of the linear analysis (15) (solid black curve) is shown together with the result of full nonlinear computations for W = 0.2 (red dotted), W = 0.4 (orange dashed), W = 0.6 (green dash-dotted) and W = 0.8 (blue dash-double dotted). Note that two equilibria exist here, but that one corresponds to a higher energy, and hence is unstable.
2.2.3 The energy of interaction. For another perspective on the problem, we consider the energy of the system, U2, which is given in dimensional terms by
image file: c6sm00901h-t13.tif(14)
For the case of linear deformations, this expression may be evaluated and expressed in dimensionless terms as
image file: c6sm00901h-t14.tif(15)
This expression is compared to the fully nonlinear calculations at the equilibrium separation, deqm, in Fig. 2b. From this plot we observe that as the electrical repulsion parameter C2 increases, the depth of the energy well in which the system sits actually decreases: increasing the strength of repulsion decreases the binding that surface tension and gravity are able to supply until ultimately the particles disperse, separating to d = ∞. (This unbinding happens for C2/W2 ≳ 0.184 according to the linear calculation.) We also note that the small deformation (linear) theory is able to give a very good qualitative account of the results of the nonlinear computations provided that the weight per unit length, W, does not become too large. However, the general trend is that, once nonlinear deformations are accounted for, the binding energy is larger (since, as we already saw, the equilibrium separation is smaller).

3 The 2-D many-body problem

3.1 Governing equations

In the last section, we considered the two-body problem in some detail. This allowed us to identify the important dimensionless parameters as the dimensionless weight per unit length, W, and the dimensionless repulsion strength, C2. Furthermore, we showed that in the limit of light particles, W ≪ 1, there is an equilibrium floating arrangement only if the charge-to-weight ratio C2/W2 ≤ 0.184. Finally, we found that the typical energy scale of the interaction in such an equilibrium floating arrangement is W2: as we expect the energy of interaction is small when the weight of the particles themselves is small.

We now turn to the many-body problem: does the presence of many floating objects cause a raft of charged particles to float deeper in the liquid than would be the case without many-body interactions? If so, how does this ‘collective sinking’ influence the typical energy well in which each particle sits?

We consider the same setup as for the two-body problem but with N line charges, i.e. N line charges, each of mass m and charge q per unit length float at a liquid–fluid interface, as shown schematically in Fig. 1. (For simplicity, we shall consider N = 2n + 1 odd, which facilitates our calculations; we do not expect this restriction to have any material effect, especially for large N.) The position of each particle in this ‘raft’ is determined by the balance of forces in both the vertical and horizontal directions.

In the horizontal direction force balance requires the net horizontal force from surface tension on the ith particle, γ(cos[thin space (1/6-em)]βi+ − cos[thin space (1/6-em)]βi), to balance the horizontal component of the electrical repulsive force arising from every other particle. In dimensionless terms we have

image file: c6sm00901h-t15.tif(16)
where di,j = [(xjxi)2 + (yjyi)2]1/2 is the distance between two particles. (Note that the Coulombic repulsion between line charges ∼1/di,j with the additional factor (xjxi)/dij coming from resolving the force in the horizontal direction.)

In the vertical direction, the restoring vertical force from surface tension on the ith particle, γ(sin[thin space (1/6-em)]βi+ + sin[thin space (1/6-em)]βi), must balance both the weight of the particle, mg, and any vertical component of the repulsion between them. We have in dimensionless terms

image file: c6sm00901h-t16.tif(17)
The equations representing force balance give 2N equations for 4N unknowns (for each particle, we know neither its (xi, yi) coordinates nor the meniscus angles on either side of it, βi±). To determine additional equations, we must also obtain additional relationships for the βi±. In principle, these angles may be determined by solving the Laplace–Young eqn (1) subject to the boundary conditions h(xi) = yi, h(xi+1) = yi+1. In practice, this calculation is made simpler by two observations: firstly, for the menisci that extend to ±∞ we may use the first integrals (11) and (12). Secondly, for particles that are close (in comparison to the capillary length [small script l]c) the meniscus may be approximated by the arc of a circle, with radius of curvature R (see lower right inset of Fig. 1a). The radius of curvature of the meniscus between particles i and i + 1, which we denote Ri,i+1, is then determined by noting that the hydrostatic pressure within the liquid along the interface, which we estimate as −(ρlρf)g(yi + yi+1)/2, must be balanced by the capillary pressure drop, γ/Ri,i+1; in dimensionless notation we therefore have
image file: c6sm00901h-t17.tif(18)
Elementary geometry then leads to expressions for βi+ and βi+1 in terms of the particle positions and Ri,i+1:
image file: c6sm00901h-t18.tif(19)
image file: c6sm00901h-t19.tif(20)
With the geometrical relationships (19) and (20), the pertinent results from the Laplace–Young eqn (11), (12) and (18) and the two force balance eqn (16) and (17), we have a closed problem. We solve these equations numerically using Newton's method: an initial guess for the position of the particles (xi, yi) is supplied, which is refined by the iteration step until a convergence criteria for (xi, yi) is met. More details about the numerical method are discussed in Appendix A.

3.2 Numerical results

For rafts consisting of identical particles, there are three parameters that characterize the raft shape: the dimensionless weight per unit length, W, the dimensionless charge per unit length, C, and the number of particles in the raft, N. For given values of these parameters, the theoretical formulation given in the preceding section allows us to determine numerically the position of each particle and the properties of the raft.

Our main interest lies in the effect of varying the number of particles N in a raft, and in understanding how a large number of particles behave collectively. However, it is also of interest to see how, with a fixed number of particles, a raft behaves as the two physical parameters, namely W and C, are changed. Fig. 3a shows how the raft shape changes as W increases. As should be expected, the particles fall deeper into the supporting liquid as they grow heavier, becoming more closely packed as they do so. However, we emphasize that this process is highly nonlinear: the largest W used in Fig. 3a is within 20% of the smallest value and yet the maximum depth of the raft increases by almost a factor of 2 and the particles come significantly closer together. This nonlinearity is a result of the ‘collective sinking’ of the particles: an increase in W brings them closer together, decreasing the vertical supporting force that surface tension is able to provide, causing the particles to lower themselves further into the liquid to achieve that supporting force and, in the process, increasing the attractive force between them.

image file: c6sm00901h-f3.tif
Fig. 3 The effect of changing the physical parameters on the shape of a raft with a fixed number of point-like particles at an interface (here N = 31). (a) Increasing the weight per unit length of each particle, W, causes the raft to sink deeper into the supporting liquid; points show the position of the particles with C = 0.02 (fixed) and W = 0.0602 (●), W = 0.063 (■), W = 0.066 (▲), W = 0.069 (◀) and W = 0.0714 (×). (b) Increasing the charge per unit length of each particle, C, causes the raft to lift out of the supporting liquid (since the electrical repulsion is stronger, the equilibrium distance is larger and particles can reach equilibrium without sinking so deep into the liquid); points show the position of the particles with W = 0.0602 (fixed) and C = 0.0155 (●), C = 0.0195 (■), C = 0.0214 (▲) and C = 0.0228 (◀). In each case, the interface shape is shown by the solid black curves.

Fig. 3b confirms the important role of this ‘collective sinking’ in determining the raft shape: as the charge carried by each particle increases, the distance between those particles also increases (since the repulsive force increases). This means that the vertical surface tension force required to support the particles can be obtained at a lower depth and so the raft rises out of the lower liquid.

To understand better the role of collective sinking, Fig. 4 shows the effect of changing just the number of particles contained in the raft. For very small rafts, N = 3 for example, the interface is barely deformed, and the equilibrium particle separation is relatively large: this is to be expected since the weight per unit length used here, W = 0.0602 ≪ 1, does not lead to a significant lateral capillary force. However, as more of these lightweight, lightly charged, particles are introduced (i.e. N increases), the particles float significantly lower in the liquid (Fig. 4a) and come much closer together (the mean separation between neighbours, [d with combining macron], decreases, as shown in Fig. 4b). We see then that the attractive capillary interaction between neighbours must be becoming stronger with increasing N, since the repulsive electrostatic interaction between neighbours increases as [d with combining macron] decreases.

image file: c6sm00901h-f4.tif
Fig. 4 The effect of changing just the number of particles at the interface. With W = 0.0602 and C = 0.02 fixed, we see that as the raft size grows, particles are not only more closely packed (on average) but also sink lower into the supporting liquid. (a) Particle positions for rafts with increasing N: N = 3 (●), 9 (■), 15 (▲), 21 (◀), 27 (▼), 33 (▶), and 39 (×). (b) Mean distance between neighbouring particles for rafts with W = 0.0602 and C = 0.02. Vertical error bars indicate the standard deviation.

We note that we are not able to find equilibrium cluster shapes with arbitrary values of the weight per unit length W or number of particles, N, in a cluster. In particular, for large, heavy clusters (N and W both large) our algorithm fails to find equilibrium configurations. We interpret this apparent lack of equilibrium solutions as a transition from floating to sinking, as has been observed at macroscopic scales with sufficiently large, heavy particle rafts.22–24 While this is interesting at a macroscopic scale, we do not study this transition here since this is extremely unlikely to be pertinent at microscopic scales. Instead we focus on how the lowering of the cluster in the liquid (but without becoming immersed in the bulk) modifies the interaction between floating particles. To understand how this ‘collective sinking’ can enhance the strength of lateral capillary interactions, we turn to some scaling considerations.

3.3 Scaling analysis and typical energy of interaction

To understand how collective effects can enhance lateral interactions, we study the total energy of the system in scaling terms. This energy consists of gravitational energy (of displaced liquid and particles), interfacial energy and the electrostatic energy of the particles. In scaling terms, the lateral extent of the raft L[d with combining macron]N; as found in the simulations presented in Fig. 3 and 4 we shall assume that L ≪ 1 (i.e. that the raft is small compared to the capillary length). The gravitational and interfacial energy of the liquid displaced by the raft itself and the outer meniscus ∼H2(1 + N[d with combining macron]), where H is the depth of the edge of the raft and ‘∼’ means “scales as”. The gravitational energy of the particles themselves ∼NWH. Finally, the electrical potential energy ∼−C2N2[thin space (1/6-em)]log[thin space (1/6-em)][d with combining macron] since there are N particles, each of which interacts with N − 1 other particles (and where we neglect terms like N2[thin space (1/6-em)]log[thin space (1/6-em)]N since they do not vary with [d with combining macron]). The total energy is then
UNH2(1 + N[d with combining macron]) + WNHC2N2[thin space (1/6-em)]log[thin space (1/6-em)][d with combining macron],(21)
which may be minimized by varying H and [d with combining macron] simultaneously. This minimization gives that [d with combining macron]NC2/H2 and HNW/(1 + [d with combining macron]N) ∼ NW, since N[d with combining macron] ≪ 1, so that image file: c6sm00901h-t20.tif. Comparing this to the corresponding result from the two-body problem, (9), and assuming that the prefactor might be such that this result holds all the way down to N = 2, we then hypothesize that
image file: c6sm00901h-t21.tif(22)
This prediction is compared with our numerical data in Fig. 5a. The comparison shows that numerical data collapse as C, W and N are varied independently; furthermore the scaling and, to a certain extent, the prefactor predicted in (22) are as predicted. (Note that strictly speaking the above scaling analysis applied only to large numbers of particles, N ≫ 1, and so the prefactor in (22) is meant to be indicative.)

image file: c6sm00901h-f5.tif
Fig. 5 The (a) mean particle separation and (b) energy change per particle as a result of collective sinking. Results are shown for rafts with different numbers of particles, N, and different weights per unit length, W. Colour is used to show the number of particles in the raft from dark red (N = 3) to blue (N = 39) while the symbol indicates different values of W: squares show numerical results with W = 0.0602 and N varying while triangles show individual values of N (coded by colour) with W varying in the range 0.0602 ≤ W ≤ 0.66. Here C = 0.02 throughout and the solid lines represent the predictions (a) (22) and (b) (24), which are based on our scaling analysis and comparison with the two-particle problem.

The total energy of interaction of the system in this equilibrium, UN, is also of interest. Using the equilibrium value [d with combining macron] from (22), we find that

UNN2[W2 + C2(1 + log[thin space (1/6-em)]N)] ∼ N2W2,(23)
for C/W ≪ 1. The key observation about this energy is that it is quadratic in N, which means that the binding energy per particle, UN/NN increases with the number of particles: the collective sinking of particles into the liquid increases their binding energy. Put another way, an estimate of the binding energy that focuses only on the energy of pairwise interactions is qualitatively incorrect as the size of the raft increases.

Comparing the scaling in (23) with the exact result for N = 2, (15), and assuming that the prefactor is such that the former scaling with N = 2 reduces to (15) we find that

image file: c6sm00901h-t22.tif(24)
Fig. 5b shows a plot of the binding energy (per particle) determined numerically as the physical parameters of the system are varied. Again, we see that the data collapse using the scaling suggested by (24), and that the binding energy per particle exhibits a similar scaling to that expected from (24), which is shown as the solid line in Fig. 5b. As expected, therefore, we see that collective sinking can cause an increase in the ‘energy well’ in which floating particles find themselves.

We emphasize that the pair-wise calculation, which led to (15), would predict an energy per particle ∼W2. Collective sinking (and also the linear superposition of capillary collapse27,28) leads to an additional multiplicative factor N, which clearly becomes more important as the size of the cluster, N, increases. In particular, while the scaling in (23) holds, the binding energy per particle can become arbitrarily large as a result of collective sinking.

4 The 3-D case

4.1 Scaling analysis

In the three-dimensional case of interacting dipoles that motivated this study it is difficult to perform full numerical calculations: these would involve determining the three-dimensional meniscus shape surrounding many objects and resolving a contact line that is not in general circular, even for floating spheres.35 While such an investigation remains a possibility for the future, here we focus on using the understanding we have gained from our detailed consideration of the two-dimensional problem to understand the 3D problem using scaling arguments.

We consider N dipoles, each of mass m. Assuming that the dipoles are aligned by an external field so that they are repulsive, not attractive, the pairwise interaction energy may be written UA/d3, where A is a constant that will depend on the nature of the dipolar interaction, e.g. A = μ0|[p with combining right harpoon above (vector)]mag|2 for magnetic dipoles or A = |[p with combining right harpoon above (vector)]elec|2/ε0 for electrical dipoles. The scaling behaviour of the dipolar energy of this assembly deserves careful discussion: in the 2D case the sum of pairwise interaction energies meant that the total energy scaled like N2. For dipoles in 3D, however, the scaling is more subtle since the energy of an individual dipole surrounded by an infinite, planar cloud of dipoles with mean nearest-neighbour spacing [d with combining macron] is found by summing over the interaction energies of an infinite series of rings of radius Ri = i[d with combining macron] (i = 1, 2, 3,…), each containing approximately 2πi other dipoles. We therefore find that image file: c6sm00901h-t23.tif, and that the energy of the system due to these dipolar interactions is UdipoleNA/[d with combining macron]3. The gravitational energy of the particles is UparticlesNmgH, while the gravitational (and interfacial) energy of the liquid due to the deformation is Uliquid ∼ ΔρgH2[N1/2[d with combining macron][small script l]c + N[d with combining macron]2], where we have included the displaced liquid from the aggregate itself as well as the external meniscus around the perimeter and Δρ = ρlρf.

Minimizing over H, we find that

image file: c6sm00901h-t24.tif(25)
while minimizing over [d with combining macron] gives
image file: c6sm00901h-t25.tif(26)
Solving for [d with combining macron] gives
image file: c6sm00901h-t26.tif(27)
assuming that N1/2[d with combining macron][small script l]c and introducing an unknown prefactor α. Note that as in the 2D monopole case, the mean separation decreases as the aggregate grows larger.

We emphasize that this result only holds for large clusters, where each dipole effectively has infinitely many other dipoles with which it could interact; the interaction energy is then cut off by the decay of the dipolar potential, rather than the number of neighbours. With smaller clusters, the dipole–dipole interaction energy is instead limited by the number of available dipoles, which are a typical distance [r with combining macron] away. In this case, UdipoleAN2/[r with combining macron]3. Assuming that [r with combining macron]N1/2[d with combining macron], we have that UdipoleAN1/2/[d with combining macron]3. To progress further, we assume that small clusters are approximately spherical,36 with radius of curvature ∼[r with combining macron] so that the surface energy ∼γN[d with combining macron]2; equating with the dipole–dipole energy, we find that [d with combining macron]N−1/10.

4.2 Macroscopic analogue experiments

We are not aware of experimental data at a microscopic scale that would allow the scaling law in (27) to be tested. However, recent experiments on macroscopic paramagnetic spheres floating at an air–water interface showed that these spheres do form a raft of the form considered in this paper.37 By digitizing the images presented in Fig. 2 of Vandewalle et al.,37 we were able to compute the mean particle separation (measured only between nearest neighbours) from these experiments on aggregates with varying numbers of particles, N. We expect to see the particle separation decreasing with increasing N, and more specifically, according to (27), that [d with combining macron]N−1/4. This scaling is confirmed by the results presented in Fig. 6b with the prefactor for this scaling corresponding to α ≈ 3. We also note that while for small cluster sizes, N ≲ 40, the data presented in Fig. 6b appear to flatten out slightly, this is not as much as might be expected on the basis of the [d with combining macron]N−1/10 scaling discussed above for this limit.
image file: c6sm00901h-f6.tif
Fig. 6 Aggregates of dipolar particles at a liquid–fluid interface. (a) N = 94 paramagnetic spheres of radius a = 200 μm in an externally applied field of 50 G form a closed aggregate at an air–water interface (image taken from Vandewalle et al.37 with permission from Springer). Note in particular that the particles near the edge are more spaced than those at the centre where, since the interfacial deflection is larger, the capillary attraction is larger too. (b) The average distance [d with combining macron] between paramagnetic particles floating at the air–water interface (data taken from Vandewalle et al.37). The dashed line shows the best fit from the prediction (27) – a scaling that is tested further in the inset. (c) PS particles of radius a = 1μm at the interface between decane and a 0.1 M aqueous NaCl solution. Again, we note that near the edge of the cluster the particle spacing becomes larger than it is away from the edge – an observation that is quantified in (d) by plotting the variation of particle density along a normal to the edge of the aggregate (highlighted by the red lines in (c)).

4.3 Relevance to microscopic scenarios

Having seen that the scaling law in (27) is able to predict the increased clustering of macroscopic dipolar particles as the size of the aggregate increases, we now extrapolate this scaling law to the microscopic scale that motivated this work: colloidal particles at a liquid–fluid interface. In this setting, the key question is how the typical energy of interaction per particle,
image file: c6sm00901h-t27.tif(28)
compares with the thermal energy, kBT. Based on this scaling law, we see that for this particle-level interaction energy to be larger than the thermal energy, we must have NNc where
image file: c6sm00901h-t28.tif(29)
To make further progress, we need to estimate the size of the various terms in (29). From detailed studies of the pair-correlation function for colloidal layers of PS particles (of radius a = 1 μm and density ρPS = 1050 kg m−3) it has been suggested5 that A ≈ 3 × 105kBT μm3. It is very difficult to be certain of the value of the constant of proportionality α from the analogue magnetic experiment (which is likely to depend on, among other things, the wetting properties of the particles), and we see that the scalings above vary sensitively with α. Therefore we take 0.1 ≲ α ≲ 1 for now, and note the typical ranges of the parameters based on this. We also note that the real experiments of interest occur at an oil–water interface with ρoil = 730 kg m−3, ρwater = 1000 kg m−3 and γ = 52 mN m−1 (see Zeppieri et al.,38 for example); as such, we expect that the driving mass will be image file: c6sm00901h-t29.tif while the appropriate capillary length [small script l]c = [γρg]1/2 with Δρ = ρwaterρoil ≈ 270 kg m−3.

Using the values above, we find with α = 0.1 that Nc ≈ 3 × 104, i.e. an aggregate around 100 particles in each direction should be stable to thermal noise. While large, this number of particles is not infeasible. If instead α = 1 then the critical number of particles in an aggregate is Nc ≈ 3 × 108, which is so large as to be very difficult to observe.

Another test of the scaling laws is the values predicted for the two spatial scales of the aggregate: the mean inter-particle separation, [d with combining macron], and the typical depth of sinking, H. By construction, the value of [d with combining macron] at N = Nc is [d with combining macron] ∼ (A/kBT)1/3 ∼ 70 μm; this is the distance at which the typical electrostatic interaction becomes on the same order as the thermal energy, and so in real aggregates the particle separation is likely to be considerably smaller. More interesting is the prediction from (25) that around N = Nc ≈ 3 × 104 the depth of the aggregate H ∼ 3 nm (using the prefactor α = 0.1 in (27)); with α = 1 we find H ∼ 300 nm. These depths are significantly smaller than the O(10 μm) depths predicted from a previous mean-field model25 and, as yet, not detected. The cluster depths we predict are too small to be directly imaged in microscopy, but should be amenable to optical interferometry or FreSCa cryo-SEM.39

While the importance of collective sinking in aggregates of interfacial colloids remains purely speculative, we can compare the phenomenology of our own experiments with what would be expected on the basis of the collective sinking hypothesis. In particular, the collective sinking hypothesis suggests that isolated clusters of colloids can form and, further, that when they do the particles near the edge of the cluster/aggregate should be more widely spaced than those near the centre of the aggregate. (This is observed both in our numerical simulations, see for example Fig. 4a, and in experiments on macroscopic paramagnetic particles floating at the interface,37Fig. 6a.) Similarly, we are able to see a similar phenomenology in aggregates of PS particles trapped at the interface between decane and aqueous salt solutions (see Fig. 6c). These experiments follow the methods of Parolini et al.5 for purification of reagents. After long periods of equilibration (sometimes overnight or even after a few days) regions of crystalline arrangements coexist with completely empty regions see Fig. 6c. Here, we highlight the clear edge of the cluster (the red line in Fig. 6c) and plot the variation of density with distance normal to this interface at isolated points (highlighted in Fig. 6c). This analysis reveals (Fig. 6d) that there is a more than two-fold decrease in particle density from the bulk of the crystal to the edge. We are not aware of any other explanation for either the existence of a well-defined edge of a cluster or for this spacing, and will explore this phenomenon in detail in a separate work.

5 Conclusions

In this paper, we have presented a toy model of the interaction between repulsive particles at an interface. This model allowed us to consider the interaction of large numbers of particles at an interface and to show that as the number of particles increases the particles actually become more closely bound together. This effect is due to the collective sinking of the particles into the liquid: the proximity of other interfacial particles means that the interface is less curved locally than it would otherwise be and so particles sink lower into the liquid. This in turn increases the magnitude of the attractive interaction between them; while this is qualitatively similar to the capillary collapse studied previously, our detailed calculations with two particles showed that this collective sinking provides a binding energy that is quantitatively stronger than would be predicted by using a linear superposition argument.27,28 Crucially, we expect that the importance of this collective sinking, and the additional binding provided by it, increases as the size of the cluster increases.

We presented detailed numerical results for the flotation of line charges. This allowed us to readily perform detailed numerical simulations of the problem, and to gain understanding that could be translated into a scaling argument and thence into scaling arguments for the problem of several dipolar spheres interacting. To make our models more realistic would require detailed simulations of the meniscus around an array of spherical particles. While this would be an involved procedure, we believe that it may soon be feasible computationally35,40 and, further, may yield new insight beyond existing mean-field theories.25,28 In particular, these mean-field theories use a linear superposition of the far-field, small deflection meniscus around an axisymmetric object, h(r) ∼ K0(r/[small script l]c), even though close to small axisymmetric objects a subtly different meniscus form is more appropriate.41 This subtlety arises from a balance between the nonlinear curvature terms and suggests further that the linear superposition approach may not be valid here, particularly when the particles approach one another on a scale comparable to their radius. At still closer approach, the effects of the particle roughness may become important;13 we do not expect roughness to play a major role, however, since clustering has been observed with particles that are well-separated compared to the particle roughness scale.21

For simplicity, our model did not include the electro-dipping force that is believed to be important in at least some observations of colloidal self-assembly. As a result, our theoretical predictions are unlikely to be directly applicable to colloidal self-assembly. Nevertheless, the mechanism that we have investigated here should be important regardless of what causes the force normal to the interface. In particular, while the gravitational contribution discussed here may not be dominant in all situations of interest, a similar effect will exist with an electro-dipping force. We hope that our calculations and scaling arguments will motivate further detailed study of this possibility.

Appendix A: numerical method

To solve the equations of vertical and horizontal force balance, (16) and (17), we used Newton's method. We firstly arrange the equations to take the form: image file: c6sm00901h-t30.tif where image file: c6sm00901h-t31.tif is the vector of particle positions for half of the raft (using symmetry). The set of particle positions, X*, that solves the vector function F(X*) = 0 is obtained by starting from an initial guess, X(0) and repeating the iteration scheme
X(n+1) = X(n)J−1F(X(n))(30)
where J is the Jacobian matrix of F(X), i.e. Jij = ∂Fi/∂xj. This iteration continues until the maximal element of F(X(n)) (in absolute terms) is below some residual, which we set to be ε = 10−13 here.


We acknowledge many discussions with L. Parolini and useful comments on an earlier draft of this paper from Martin Oettel. This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-2012R1A6A3A03039558) and the Korea Institute of Machinery and Materials (KIMM) under Grant NK196D.


  1. P. Pieranski, Phys. Rev. Lett., 1980, 45, 569–572 Search PubMed.
  2. R. K. Kalia and P. Vashishta, J. Phys. C: Solid State Phys., 1981, 14, L643–L648 Search PubMed.
  3. A. J. Hurd, J. Phys. A: Math. Gen., 1985, 18, L1055–L1060 Search PubMed.
  4. L. Isa, K. Kumar, M. Müller, J. Grolig, M. Textor and E. Reimhult, ACS Nano, 2010, 4, 5665–5670 Search PubMed.
  5. L. Parolini, A. D. Law, A. Maestro, D. M. A. Buzza and P. Cicuta, J. Phys.: Condens. Matter, 2015, 27, 194119 Search PubMed.
  6. R. Aveyard, B. P. Binks, J. H. Clint, P. D. I. Fletcher, T. S. Horozov, B. Neumann, V. N. Paunov, J. Annesley, S. W. Botchway, D. Nees, A. W. Parker and A. D. Ward, Phys. Rev. Lett., 2002, 88, 246102 Search PubMed.
  7. K. Masschaele and J. Vermant, Soft Matter, 2011, 7, 10597–10600 Search PubMed.
  8. M. P. Boneva, N. C. Christov, K. D. Danov and P. A. Kralchevsky, Phys. Chem. Chem. Phys., 2007, 9, 6371–6384 Search PubMed.
  9. M. P. Boneva, K. D. Danov, N. C. Christov and P. A. Kralchevsky, Langmuir, 2009, 25, 9129–9139 Search PubMed.
  10. F. Ghezzi and J. C. Earnshaw, J. Phys.: Condens. Matter, 1997, 9, L517–L523 Search PubMed.
  11. J. Ruiz-Garcia, R. Gámez-Corrales and B. I. Ivlev, Physica A, 1997, 236, 97–104 Search PubMed.
  12. J. Ruiz-García, R. Gámez-Corrales and B. I. Ivlev, Phys. Rev. Lett., 1998, 58, 660–663 Search PubMed.
  13. D. Stamou, C. Duschl and D. Johannsmann, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 62, 5263–5272 Search PubMed.
  14. M. Quesada-Pérez, A. Moncho-Jordá, F. Martínez-López and R. Hidalgo-álvarez, J. Chem. Phys., 2001, 115, 10897–10902 Search PubMed.
  15. M. G. Nikolaides, A. R. Bausch, M. F. Hsu, A. D. Dinsmore, M. P. Brenner, C. Gay and D. A. Weitz, Nature, 2002, 420, 299–301 Search PubMed.
  16. D. Y. C. Chan, J. D. Henry and L. R. White, J. Colloid Interface Sci., 1981, 79, 410–418 Search PubMed.
  17. P. A. Kralchevsky and K. Nagyama, Adv. Colloid Interface Sci., 2000, 85, 145–192 Search PubMed.
  18. N. Sharifi-Mood, I. B. Liua and K. J. Stebe, Soft Matter, 2015, 11, 6768–6779 Search PubMed.
  19. D. Vella and L. Mahadevan, Am. J. Phys., 2005, 73, 814–825 Search PubMed.
  20. M. G. Nikolaides, A. R. Bausch, M. F. Hsu, A. D. Dinsmore, M. P. Brenner, C. Gay and D. A. Weitz, Nature, 2003, 424, 1014 Search PubMed.
  21. F. Bresme and M. Oettel, J. Phys.: Condens. Matter, 2007, 19, 413101 Search PubMed.
  22. D. Vella, P. D. Metcalfe and R. J. Whittaker, J. Fluid Mech., 2006, 549, 215–224 Search PubMed.
  23. M. Abkarian, S. Protière, J. M. Aristoff and H. A. Stone, Nat. Commun., 2013, 4, 1895 Search PubMed.
  24. D. Vella, Annu. Rev. Fluid Mech., 2015, 47, 115–135 Search PubMed.
  25. V. M. Pergamenshchik, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 011407 Search PubMed.
  26. V. M. Pergamenshchik, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 021403 Search PubMed.
  27. J. A. Domínguez, M. Oettel and S. Dietrich, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 011402 Search PubMed.
  28. J. Bleibel, S. Dietrich, A. Domínguez and M. Oettel, Phys. Rev. Lett., 2011, 107, 128302 Search PubMed.
  29. J. Bleibel, A. Domínguez, M. Oettel and S. Dietrich, Eur. Phys. J. E: Soft Matter Biol. Phys., 2011, 34, 125 Search PubMed.
  30. J. Bleibel, A. Domínguez, M. Oettel and S. Dietrich, Soft Matter, 2014, 10, 4091–4109 Search PubMed.
  31. J. Bleibel, A. Domínguez and M. Oettel, J. Phys.: Condens. Matter, 2016, 28, 244021 Search PubMed.
  32. M. M. Nicolson, Proc. Cambridge Philos. Soc., 1949, 45, 288–295 Search PubMed.
  33. D. Vella, PhD thesis, University of Cambridge, 2007.
  34. E. H. Mansfield, H. R. Sepangi and E. A. Eastwood, Philos. Trans. R. Soc., A, 1997, 355, 869–919 Search PubMed.
  35. P. L. H. M. Cooray, PhD thesis, University of Cambridge, 2013.
  36. S. G. Jones, N. Abbasi, B.-Y. Moona and S. S. H. Tsai, Soft Matter, 2016, 12, 2668–2675 Search PubMed.
  37. N. Vandewalle, N. Obara and G. Lumay, Eur. Phys. J. E: Soft Matter Biol. Phys., 2013, 36, 127 Search PubMed.
  38. S. Zeppieri, J. Rodríguez and A. L. López de Ramos, J. Chem. Eng. Data, 2001, 46, 1086–1088 Search PubMed.
  39. S. Coertjens, P. Moldenaers, J. Vermant and L. Isa, Langmuir, 2014, 30, 4289–4300 Search PubMed.
  40. H. Cooray, P. Cicuta and D. Vella, J. Phys.: Condens. Matter, 2012, 24, 284104 Search PubMed.
  41. L. L. Lo, J. Fluid Mech., 1983, 132, 65–78 Search PubMed.


Present address: Department of Nature-Inspired Nanoconvergence Systems, Korea Institute of Machinery and Materials, Daejeon 34103, Korea.

This journal is © The Royal Society of Chemistry 2017