Effects of charge and size on the coadsorption of counterionic colloids in Gibbs monolayers

This study uses a coarse-grained Monte Carlo algorithm to model and simulate the coadsorption of a binary mixture of counterionic colloids in Gibbs monolayers. These monolayers form at a idealized air-water interface, with one non-soluble species conﬁned at the interface and the second one partially soluble in the aqueous phase. The investigation focuses on the effect of colloidal size and charge on the thermodynamics and microstructure of the monolayer. We ﬁnd that the composition of the monolayer evolves non-trivially with surface coverage, depending on the balance of steric and electrostatic forces. When the electrostatic interactions are weak, the soluble species is expelled from the monolayer upon compression, yielding a phase behaviour particularly sensitive to the relative size of the soluble and non-soluble colloids. By contrast, strong electrostatic interactions favour the stabilization of the soluble particles in the monolayer and the formation of quasi-equimolar ﬂuids, with only a weak dependence on particle size. The combination of these phenomena results in the formation of a number of two-dimensional mesoscopic arrangements in the monolayer, ranging from diluted gas-phase behaviour to domains of aggregates and percolates, and to incipient crystalline structures.


Introduction
Computer simulation has consolidated as a fundamental tool for the rationalization and technological exploitation of the properties of molecular liquids, soft condensed matter, and materials in general.An important step in using computer simulation tools to investigate the properties of any system is to obtain a suitable model for the problem at hand.There is no single recipe for this task, and the level of complexity of the model to be used will depend on the properties and time and length scales to be studied.Very detailed microscopic models may provide quantitative information about some properties, but they usually make it difficult to explore the general physical characteristics that guide the behaviour of the system.This latter information is typically more accessible with simplified models, which may however compromise the acquisition of reliable quantitative and atomistic information.The choice and design of appropriate models for a specific study is an art in itself, which Professor Jaroslav Ilnytskyi has excelled throughout his career, in particular within the framework of complex molecular fluids [1,2].
We present a computer simulation study inspired by the work of prof.Ilnytsky, of the thermodynamics and structural properties of a multicomponent fluid of charged particles dispersed in an interfacial monolayer.More specifically, we investigate generic cases of Gibbs monolayers representing dispersed amphiphilic particles at an air-water interface [3,4].These systems began to be studied in modern times already in the 18th century, from the celebrated observations of films of olive oil on water carried out by Benjamin Franklin [5] and later revisited by Lord Rayleigh [6].Subsequently, Agnes Pockels [6,7] designed the famous trough to measure the surface tension of monolayers, which was eventually improved by Irving Langmuir and his assistant Katherine Blodgett, implementing the efficient transfer of the monolayer to a solid substrate [8].A monolayer is said to be of Langmuir-type if all the molecules are completely insoluble, while it is said to be a Gibbs monolayer if some of its constituents are partially soluble in the bulk of at least one of the phases that share the interface.In a general situation, the species belonging to the insoluble fraction will remain confined at the interface, while those belonging to the soluble fraction will be distributed between the fluid phases and the interface [9].Gibbs monolayers find a broad range of applications, often related to sensing devices, where interfacial changes are induced when species in solution selectively bind to non-soluble colloidal receptors hosted by the monolayer [10,11].A nice example of this type of application is the interfacial binding of biomacromolecules, such as DNA fragments, by cationic surfactants.For such systems, it has been shown that the compression of the monolayer affects the concentration of the soluble biomolecules in the monolayer, a thermodynamic effect that is not yet well understood [11,12].Gibbs monolayers have also a potential to be utilized in the synthesis of thin-film topologies, facilitated by counterionic template frameworks that are incorporated to the interface from solution.An example of this is the fabrication of Langmuir-Blodgett layered structures based on guanidinium surfactants employing soluble carboxylate and phosphate moieties [10].
Not surprisingly, Gibbs monolayers have attracted much attention and have motivated important developments from the field of molecular and colloidal simulation.In this context, different levels of approximation have been proposed, ranging from models with a highly detailed atomistic description of the molecules and particles involved, including water [13][14][15], to coarse-grained models where entire parts of the molecules are modelled as beads, using various effective force fields [16,17].In many studies, the monolayer is examined as a three-dimensional system, simultaneously modelling and simulating the interface and the two coexisting phases.In this way, it is possible to obtain valuable insights into the microscopic properties of the system, especially into the conformations that favour colloidal interactions [13,15,17].However, due to the high computational cost, these models are usually restricted to small systems and short timescales, which often precludes the access to thermodynamic and mesoscopic structural properties.
In an attempt to complement the information obtained with more detailed modelling studies, we have recently proposed a simple coarse-grained model for the study of Gibbs monolayers [18].Within our approach, only the interface itself is simulated as a two-dimensional system with interacting particles represented as flat disks.Amphiphilic and soluble particles can diffuse within the interface, but only the latter ones are capable of moving in and out of the interface.The interchange of soluble particles between the interface and the bulk phase is simulated with Grand Canonical Monte Carlo moves.Another aspect of our model is that the interaction between charged particles is modelled as a dipolar interaction, an assumption that resembles the different screening of colloidal charge induced by the two immiscible phases that conform the interface, e.g., the air-water interfaces [19,20].With these ingredients, we recently studied the phase behaviour and the structural properties of monolayers of mutually attractive soluble and non-soluble colloids with the same size and absolute dipolar charge (but opposite sign) [18].For that system, an interesting sequence of disorder fluid, cluster fluid, and crystalline phases was found, confirming the potential richness of the structural landscape of Gibbs monolayers.In this article we extend the investigation by introducing asymmetry in colloidal size in the two complementary cases where either the soluble or the non-soluble component of the mixture is larger in size.In this way, we intend to isolate the size effects on interfacial behaviour from those of other colloidal properties.

Model and methods
The colloidal model and the Monte Carlo simulation methodology are similar to those employed in our previous work [18].The Gibbs monolayer at the air-water interface is modelled using a simplified coarsegrained approach, which represents the system as a two-dimensional fluid, comprising a binary mixture of disk particles with different diameters, interacting through both steric and electrostatic interactions.One of the colloidal species (hereafter denoted as ) is assumed to be soluble, while its counterionic counterpart (denoted as ) is non-soluble (or amphiphilic) and it is therefore confined at the interface.Though the bulk of the liquid phase is not explicitly considered, it is taken into account on the assumption that it is in thermodynamic equilibrium with the interface.Thus, it is assumed that the chemical potential of the soluble  particles, denoted as   , is the same in the bulk and at the interface.This chemical potential is dependent on the concentration of the particles in solution,   , through the expression where  ref  is the chemical potential in a reference state and  is the activity coefficient ( = 1 for an ideal solution).Here,  B is the Boltzmann constant, and  is the temperature.The chemical potential modulates the incorporation of  particles from the solution to the monolayer and this is emulated by means of Monte Carlo simulations in the     Π ensemble, where   denotes the fixed number of non-soluble particles in the monolayer and Π is the surface pressure.The Monte Carlo algorithm samples the equilibrium configurations with a combination of random displacements of  or  particles, with the insertion or removal of  particles, and changes in the surface area of the monolayers.The acceptance or rejection of these move attempts will depend on the interaction between the particles, the chemical potential   and the surface pressure for volume changes, following general acceptance rules based on the Metropolis algorithm [21].
The simulations were run on systems of   = 4000 non-soluble  colloids, while the number of soluble  colloids is allowed to change throughout the simulation without constraints.Each Monte Carlo cycle involved   attempts for random displacements of a particle ( or ), for changes in the surface area, and for the insertion or deletion of a particle  in or from the interface.Using the values optimized in our previous work, the displacement of a random particle was attempted with frequency   = 0.765, the changes of surface with frequency   = 0.0005, while the insertion/deletion moves are both attempted with frequency   = 0.230.A typical simulation at a given surface pressure and chemical potential starts from a random initial configuration of low surface density.To equilibrate the system, about 1 − 2 • 10 6 Monte Carlo cycles were applied.Additional 5 • 10 5 Monte Carlo cycles were then run to obtain the averages of the observables of interest.
The choice of the interaction potential between the particles is a key aspect of the model.The interaction between charged particles at an air-water interface is a far from trivial problem that has attracted the interest of many researchers [19,20,[22][23][24].A long-range  −3 dipole-dipole interaction has been proposed to describe the interaction between charged interfacial particles, as a consequence of different shielding of the Coulombic forces in the air and water sides of the interface.Here, we consider such dipolar interaction, on top of a short-range steric repulsion between the particles.Other types of interaction, such as capillary interactions, have been pointed out as relevant [20], but are not considered in our model.With all these ingredients, the interfacial colloids interact through a pair potential  () =    () +   (), expressed as the sum of short-range steric repulsions and long-range effective dipole interactions [18].
On the one hand, the short-range interaction is represented by a truncated and shifted Lennard-Jones potentials of the form Here,    = 0.5(  +   ) is the contact distance between the interacting particles, with diameters   and   , respectively. is the distance between particles.Note that since the potential is truncated at the minimum of the Lennard-Jones potential well, it is a purely repulsive contribution to the net interaction.
On the other hand, the dipolar contribution to the interaction potential reads as follows where  and  are the unit of length and energy, respectively.  is the cutoff distance of the dipole interaction, which as in our previous article is set as   = 30.Δ  and Δ  are the dimensionless dipolar charge of particles  and .The values assigned to these dipole moments in this study are 3, 6 or 9 in reduced units, as commented below.

Results
In our previous work [18], we focused on the reference situation of  and  colloidal particles with equal dipolar charge modulus and diameter.Here, we consider situations in which the colloidal particles  and  are asymmetric in size.Two complementary cases are addressed in which either the non-soluble  particles or, alternatively, the soluble  particles are larger in size.We specifically considered the case with   =  and   = 3, as well as the reverse situation where   = 3 and   = .The temperature and the chemical potential of the soluble particles were set to  * =  B / = 1 and  * = / = −8, respectively, while surface pressure was varied in the range Π * = Π 2 / = 0.001 − 3. Our previous study showed that within this range of parameters the Gibbs monolayers display a rich phase diagram.

Thermodynamic study of the monolayers
The equations of state of the systems under study are characterized in figure 1, in terms of the evolution of the surface pressure Π * = Π 2 /, the molar fraction of soluble colloids   and the internal energy per particle in reduced units,  * = /, as a function of the inverse coverage,  −1 = /  .Note that  −1 represents the interfacial area per particle;  is the total area of the interface, while   =   •  +   •  is the effective area physically occupied by the particles at the interface, with   = π 2  /4 the area of a particle of species .
The main finding of the present study is that for weak dipolar interactions (Δ = 3), the monolayer tends to be equimolar in the two colloidal components in the high dilution limit (  ≈   ≈ 0.5, as  −1 → ∞), whereas it becomes progressively enriched in the non-soluble component  upon compression.The soluble component  tends to be expelled from the monolayer, so that at high coverage the monolayer becomes essentially a monocomponent system of  colloids (  →1,   → 0, as  −1 → 0, see figure 1d).A consequence of this behaviour is that the system displays a marked asymmetry with respect to the size ratio of the  and  colloids.Figures 1a and 1g show that when the smaller colloid is the insoluble component (  = 1), the monolayer displays a steep increase in pressure and internal energy with increasing coverage, at  −1 < 10.A similar increase occurs at much higher coverages, onset at  −1 < 2, when the insoluble colloid is the larger one (  = 3).In general terms, at any given surface coverage away from the dilution limit, much higher pressures and internal energies are associated with the   = 1 fluid in comparison to the   = 3 fluid.
The exclusion of the  particles at high coverage just described for Δ = 3, can be rationalized in terms of an increase in entropy associated with the larger area available to the intefacial particles.At a higher dipolar charge, such entropic effects are balanced by the attractive interaction between  and  particles, which favours the presence of  colloids in the monolayer.Hence, the thermodynamic behavior of the monolayer changes qualitatively as the dipolar interactions become significant.This also brings the system to a much more symmetric dependence on colloid size ratio, meaning that it becomes progressively less relevant whether the larger colloids belong to the soluble component or to the insoluble one.For the present study, the isotherms with   = 3 (  = 1) and   = 1 (  = 3) are similar for Δ = 6 and become almost overlapping for Δ = 9 (see figure 1).Interestingly, the strong dipole interactions largely favour the compression of the monolayer, which proceeds with small changes in pressure and with an internal energy  * , that remains negative throughout the full isotherms and becomes more negative upon compression.Note that this is in strong contrast with the increasingly net positive  * values attained for Δ = 3. Negative values of  * are associated with dominant attractive interactions, indicating that an appreciable fraction of  colloids remains at the interface at all surface coverages.
For Δ = 6 at intermediate dilutions ( −1 > 10),  * is weakly dependent on the surface coverage [figure 1(h)].In this range, entropic effects dominate, and the  particles are still sterically ejected from the system upon compression, yielding a decrease in   qualitatively similar to that found for Δ = 3 [figure 1(e)].Interestingly however, the -particles re-enter the monolayer at sufficently high surface coverages, correlating with a steep decrease in the internal energy.In the limit of high coverage, the molar fraction reaches the values similar to those obtained in the high dilution regime,   ≈ 0.4.For the strongly interacting Δ = 9 monolayer, the molar fraction of  colloids is larger than 0.5 at high dilution, indicating that the incorporation of an excess of soluble colloids to the monolayer is favored under these conditions.The composition of the monolayer approaches then equimolarity   ≈ 0.5 at high packing fractions.

Microscopic structure of the monolayers
The thermodynamic properties described so far are intimately related to the microscopic and mesoscopic structures adopted by the particles within the monolayer.The topology and relative stability of different interfacial structural arrangements are strongly conditioned by the incorporation of soluble colloids  from the bulk phase.Figure 2 provides an overview of different structures observed in our work.The ejection of S colloids from the monolayer at Δ = 3 can be appreciated in the sequences of configurations from 3A to 3C (  = 1,   = 3) and from 3D to 3F (  = 3,   = 1).In the high density configurations 3C and 3F, the monolayer consists of a monocomponent fluid of  colloids.A trend towards the hexatic crystalline order characteristic of the two-dimensional discotic particle fluid is observed at high pressures and surface coverages (configuration 3F), although a consolidated crystalline phase is not found within the thermodynamic range included in our study for Δ = 3.When the dipolar charge is increased, the particle aggregation and cluster formation tend to occur even in the high dilution limit.This phenomenon is particularly pronounced for Δ = 9, (configurations 9A and 9D) where linear aggregates of several particles are observed, whereas for Δ = 6 (panels 6A and 6D) mainly  −  dimers are found.As the concentration increases, the formation of clusters is more evident and percolated arrangements are formed.In the more dense monolayers of the strongly interacting colloids, the percolates colapse to an incipient square crystalline arrangements.Despite the abundance of structural defects, we obtain the evidence for square cristalline symmetry in the high coverage limit (e.g.configuration 9F in figure 2).The analysis of the corresponding radial distribution functions confirms these findings.Figure 3 depicts an illustrative set of radial  − ,  −  and  −  correlation functions (   ,   and   respectively).Configuration 9F is representative of the formation of a crystalline arrangement; the position of the peaks in the correlation functions indicates that the structure corresponds to a square lattice, in which the smaller colloids are located in the intersticial sites generated by the larger colloids.The distribution functions for configurations 6C, 6F and 9C display a strong correlation at short and medium distances, resembling the domains of crystalline structure, which however does not survive at long range.These configurations deserve further investigation to ensure that they correspond to true equilibrium states and are not trapped in a metastable configuration.Note that our simulations were performed by compressing the monolayer from the high dilution limit.Metastability could in particular explain the sharp drop in   observed for the Δ = 6 and 9 cases at high coverage, and the differences in the equation of state observed in this range while comparing the monolayers qith   = 3 vs.  = 1 (figure 1).
Colloidal aggregation in the strongly interacting monolayers is intriguing.Clustering is actually a common phenomenon in the fluids of charged colloids at air-water interfaces, as known from experiments [25,26].The transition from a regime dominated by small aggregates to either mesoscopic

13601-7
aggregation or long-range percolated arrangements is subtle and challenging to predict.Here, we assess the extent of aggregation by analyzing the interconnection between neighbouring colloids as a function of surface coverage.To this end, we considered that two particles belong to the same cluster if the distance between them is smaller than    + 0.5, where    is the contact distance between the particles.Within this framework, two types of ensemble-averaged magnitudes were employed to characterize the clustering of particles.Firstly, the distribution,  () =  •  ()/  , represents the fraction of particles belonging to a cluster of size , where  () is the average number of clusters of size  and   is the average number of particles at the interface.Secondly, the clustering paremeter,   , represents the fraction of colloids embedded in clusters of any size, including dimers and larger clusters (  = >1  ()).Figure 4 depicts the dependence of   on  −1 for the fluids under study.Consistently,   → 1 at high coverages, and it decreases steadily as the system approaches the high dilution limit.Cluster dissociation as the area accessible to the colloids increases, mainly driven by entropy, is clearly apparent for Δ = 3, with   reaching the values close to 0 at high dilution.For stronger dipolar interactions, clustering is still noticeable in the expanded monolayers.In the case of Δ = 6, even for  −1 > 100, around 10% of the colloidal particles are incorporated to a cluster.For Δ = 9, the particle aggregation at high dilution still reaches remarkable values for   above 80%.The domains of clustering and percolation are characterized in greater detail in figure 5, which shows the cluster distributions  () for relevant configurations of Δ = 9 monolayer.We find that colloids are distributed in clusters with a broad distribution of sizes.It can be noted that even at high dilution (configurations 9A and 9D), the fraction of isolated colloid monomers remains below 10%, while dimers amount to almost 20%, with the majority of the remaining 70% corresponding to clusters with sizes up to  ≈ 20.In the percolation domain (configurations 9B and 9E) not more than 0.5% of the colloids stay as monomers and the cluster distribution broadens significantly, extending to sizes  ≈ 50 − 100.Notably, a coincident aggregation behaviour is found for the monolayers with (  = 1,   = 3) and with (  = 3,   = 1) when cluster distributions are compared at similar surface coverages (9A vs. 9D and 9A vs. 9D, in figure 5).These results confirm the idea that for intermediate and high values of Δ, the presence of aggregates is one of the main characteristics of Gibbs monolayers.Similar results were obtained for colloidal particles of equal size [18].One question we leave open for future work is the stability of the clusters, whether they have long lifetimes or, on the contrary, they are continuously forming and fracturing.Answering this question demands simulation techniques other than the Monte Carlo algorithm used here, such as molecular dynamics simulations.(Colour online) Fraction of colloids in clusters of a given size  for some states of the monolayer with Δ = 9 (see labels in figure 4).Solid lines correspond to states with   = 1, while dashed lines correspond to states with   = 3.

Conclusions and final remarks
This work presents a computer simulation study of the structural and thermodynamic behaviour of Gibbs monolayers in the case where the soluble and non-soluble colloids have different diameters.As such, it is an extension of a previous work devoted to the case of colloids with equal size [18].This study also corroborates the diversity of phases and structures that are possible in Gibbs monolayers depending on colloidal parameters such as size and dipole charge, as well as thermodynamic conditions.Most of these phases and structures consist of disordered fluids with the particles forming aggregates.
It is interesting to see how the two complementary situations that we have studied, namely cases in which the non-soluble colloids  are smaller (  = 1 and   = 3), or larger (  = 3 and   = 1), display the most marked differences with each other in the case of weak colloid-colloid interactions (represented here by a dipolar moment Δ = 3).This finding can be traced back to the fact that the soluble colloids  are expelled from the monolayer with increasingly surface coverage, eventually leading to an interfacial a monocomponent fluid of  particles.Stronger colloid interactions (Δ= 6, 9) promote the stabilization of  colloids at the interface, leading to comparable fractions of soluble and nonsoluble colloids at the monolayer.Under these conditions, most of these phases and structures consist of disordered fluids with the particles forming aggregates.In very dilute systems we find gas phases of free particles and dimers.In denser systems, progressively larger clusters are formed, eventually leading to percolated arrangements.In the limit of high coverages, we find evidence of square crystalline structures, but also of fluids with short-range order which, due to the presence of defects, the correlation between particles is lost at long distances.
In short, in this work we have confirmed the potential complexity of Gibbs monolayers as a function of the system parameters, and the capacity of the outlined simulation methodology to explore it.We believe that the results obtained may be useful for a better understanding of the design and tailoring of this type of two-dimensional fluids for specific technological applications.

Figure 1 .
Figure 1.(Colour online) Monte Carlo isotherms for systems with colloids of diameters   = 3,   = 1 (black lines and circles) and   = 1,   = 3 (red lines and circles) and dipole moments Δ = 3 (panels a, d and g), 6 (panels b, e and h) and 9 (panels c, f and i).The top (a to c), middle (d to f) and bottom (g to i) panels plot the surface pressure Π * = Π 2 /, the mole fraction of  colloids   and the energy per particle  * = /, respectively, as a function of the inverse of the surface coverage,  −1 .The horizontal dash-dotted lines indicate the equimolarity of the oppositely charged colloids (  = 0.5).

Figure 2 .
Figure 2. (Colour online) Representative snapshots of microscopic configurations of the Gibbs monolayers explored in this study.Red and blue disks represent the nonsoluble () and the soluble () colloids, respectively.Snapshots with labels A, B and C correspond to situations where   = 1 and   = 3, while labels D, E and F indicate the reciprocal cases, with   = 3 and   = 1.The number in the labels indicates the dimensionless dipolar charge in each case, Δ = 3, 6 or 9.The composition of the monolayer in each state can be inferred from the representation of   in the bottom panel, which shows an extract of the data included in panels (c), (e) and (f) of figure 1.

Figure 4 .
Figure 4. (Colour online) Fraction of colloids in clusters   as a function of surface coverages  −1 , for Δ = 3 (black lines and circles), 6 (red lines and squares) and 9 (blue lines and triangles).Open symbols indicate the cases with   = 1 and   = 3 while solid symbols are the reciprocal case where   = 3 and   = 1.

Figure 5 .
Figure 5. (Colour online) Fraction of colloids in clusters of a given size  for some states of the monolayer with Δ = 9 (see labels in figure4).Solid lines correspond to states with   = 1, while dashed lines correspond to states with   = 3.