Experimental methods
Preparation of particle suspensions for microscopy measurements
In our experiments we study the behaviour of three distinct types of microspheres with different surface chemistries. The three main classes of particle involve SiO2 particles (Bangs Laboratories), amine-derivatized silica particles (NH2–SiO2, referred to as ‘NH2’ or ‘\({\rm{N}}{{\rm{H}}}_{3}^{+}\)’ in the text) (microParticles) with an estimated NH2 group content of >30 µmol g−1, and carboxylated melamine formaldehyde (COOH–MF, referred to as ‘COOH’ or ‘COO−’ in the text) particles (microParticles) with a carboxyl group content of 400 µmol g−1. The particle size distributions provided by the manufacturer are shown in Supplementary Fig. 1.
For experiments on SiO2 and COOH particles in aqueous solution, particles were first rinsed (centrifuged and resuspended) in deionized water. Next, they were incubated in 5 mM NaOH (99.99%, Alfa Aesar) solution for 10 min. Following this, they were centrifuged and resuspended in aqueous electrolyte of the required ionic strength around six times until the measured electrical conductivity of the supernatant solution converged to that of the pure electrolyte. Note that in general NaOH treatment is not essential, and overnight exposure to deionized water with subsequent rinsing in deionized water is an equally effective treatment prior to clustering experiments (Supplementary Fig. 11). However, NaOH pretreatment was necessary to observe strong clustering in COOH particles.
NH2–SiO2 particles were first rinsed in deionized water and then resuspended multiple times in aqueous electrolyte solution until the measured supernatant conductivity converged to that of the pure electrolyte. The suspension was further sonicated for cases in which a large population of ‘sticking’ particles were observed. The presence of stuck particles in the experimental data gives rise to small ‘dimer’ peaks at interparticle separations 2R in the measured g(r) which cannot be entirely eliminated (indicated as ‘d.p.’ in Figs. 2b, 3c and 5c). In experiments on mixtures of SiO2 and COOH particles (Fig. 4), the two types of particles were initially mixed at a 1:1 ratio, then incubated in 10 mM NaOH solution for 10 min. Thereafter the procedures were the same as described above.
For experiments on colloidal dispersions in alcohols, NH2–SiO2 particles were first centrifuged and resuspended in deionized water, followed by resuspension in either ethanol (≥99.8%, Sigma-Aldrich) or IPA (≥99.5%, Sigma-Aldrich). The process of centrifugation and resuspension was repeated multiple times until the value of the supernatant conductivity converged to that measured for the pure alcohol. COOH particles were first centrifuged and resuspended in deionized water, followed by resuspension in ethanol, and final resuspension in either ethanol or IPA for measurements.
Preparation and characterization of electrolyte solutions
For experiments examining the dependence of interparticle interactions on the ionic strength of the electrolyte (Fig. 1), various concentrations of NaCl (99.998%, Alfa Aesar) solution were prepared in deionized water; the measured conductivity of these solutions corresponds to a background concentration of monovalent ions of c0 ≈ 5 μM. The ionic strength of the various electrolyte solutions in our experiments was determined from measurements of electrical conductivity, s, performed with a conductivity meter (inoLab Cond 7110). A calibration curve of standard solutions was used for this purpose (Supplementary Fig. 2).
To convert the measured electrical conductivity to a background salt concentration in alcohols, we used the same calibration relationship as for aqueous electrolytes as shown in Supplementary Fig. 2, but corrected the inferred concentrations for the viscosity of the alcohol as suggested in previous work51 (the viscosity values used for ethanol and IPA were 1.1 cP and 2.4 cP, respectively, see Supplementary Information, section 1.2).
In experiments exploring the pH dependence of interparticle interactions (Fig. 2), the pH of the electrolyte was adjusted to the desired value by adding either HCl (99.999%, Alfa Aesar) or a Tris buffer (≥99.9%, Carl Roth), to deionised water. Addition of acid or buffering agent to deionized water raised the conductivity of the solution to a value between 1 and 30 μS cm−1 (0.01–0.25 mM) depending on the target pH value. For experiments performed at variable pH, ionic strength in the electrolyte was maintained constant (to within ±0.02 mM) across the entire range of pH in a given experimental series via the addition of a variable amount of NaCl. The pH value of the aqueous solution was taken as the mean value of three consecutive measurements using a pH meter (Horiba PH-33). The pH of pure alcohol samples was inferred by extrapolation of pH values measured for water–alcohol mixtures (Supplementary Fig. 2b).
Layer-by-layer coating of silica particles with polypeptides and polyelectrolytes
In the experiments presented in Fig. 3, we used alternating coatings of positively and negatively charged polyelectrolytes on plain SiO2 particles. Coatings were applied in the following pairs of combinations of positively and negatively charged polyelectrolytes: poly-K (Mr, ≥300,000; Sigma) and poly-E (Mr, 50,000–100,000; Sigma); PDADMAC (Mr, 200,000–350,000; Aldrich) and PSS (Mr, ∼70,000; Aldrich), and finally, PEI (Mr, ∼750,000; Sigma) and PSS.
To coat the particles, plain SiO2 particles were first centrifuged and resuspended in deionized water, followed by incubation in 5 mM NaOH solution for 10 min, resuspension in deionized water, and repetition of the resuspension process until the supernatant conductivity converged to that of deionized water. The rinsed particles were then incubated in the polyelectrolyte solution at 0.1% w/v for 20 min with occasional vortexing to improve mixing. The coated particles were centrifuged and resuspended in deionized water to remove any excess polymer and the resuspension procedure repeated until the conductivity of the supernatant no longer changed. Subsequent layers of polymer coatings were applied by repeating the coating procedure described above with the corresponding oppositely charged polyelectrolyte. The sign of the surface charge of each coating layer was confirmed by zeta-potential measurements (Zetasizer Nano Z, Malvern Panalytical).
Cuvette preparation and sample loading
We used a glass cuvette with a polished flat-well of 1 mm depth (20/C/G/1, Starna Scientific), as shown in Supplementary Fig. 3b, for all video microscopy measurements. The cuvette was cleaned with piranha solution (3:1 mixture of concentrated sulfuric acid and 30 wt% hydrogen peroxide solution) and then rinsed thoroughly with deionized water. A glass cuvette naturally provides a negatively charged surface in water, ethanol and IPA for experiments with negatively charged particles. For experiments with positively charged particles, the entire cuvette was coated with 1% w/v PEI solution, rinsed and dried under nitrogen to provide a thin layer of positively charged polymer coating. To load the cuvette, the prepared particle solution was carefully pipetted into the well and sealed with the cover slide such that the device was free of air bubbles and held together by capillary force.
Microscopy
The optical microscope was constructed using a 470 nm light-emitting diode (LED) (M470L4, ThorLabs), a 10× objective (Olympus UPlanSApo) and a charge-coupled device camera (DCU223M, ThorLabs) for recording images (Supplementary Fig. 3). The sample holder was placed onto a carefully balanced pitch and roll platform (AMA027, ThorLabs). Following complete settling of particles in suspension to a plane near the bottom surface of the cuvette, which typically takes about 2 min, the focus was adjusted such that a clear intensity maximum was observed for all particles. All measurements were performed after complete settling. The intensity of the LED was adjusted such that the intensity maxima of illuminated particles did not exceed the saturation value of the camera, enabling accurate particle localization.
Video recording and data processing
Sequential images of the 2D suspension of colloidal particles were taken with ThorCam software at a constant frame rate of 5, 10 or 30 frames per second for 150–500 frames using an exposure time of ≈ 0.5 ms. The images were processed based on the radial symmetry method using the TrackNTrace particle-tracking framework, where the particle centre maximum is detected52,53. The localization precision for a static SiO2 particle during a 100 s measurement at an exposure time of ≈ 0.5 ms was found to be <20 nm, as shown in Supplementary Fig. 4. In the analysis of experimental images, coordinates of all particle centres were extracted from the recorded frames, and the radial distribution function curve g(r) calculated and averaged over all images. To clearly distinguish between experiments on particles with different signs of particle charge in three different solvents, the recorded images were digitized and false-coloured (Supplemenetary Information, section 1.6 and Supplementary Fig. 3). The average particle detection efficiency over all experiments was >98%.
Simulation methods
BD simulations of interparticle interactions
We performed BD simulations of a 2D distribution of spheres interacting via an appropriately chosen input potential using the BROWNIAN package in the LAMMPS software54. We inferred the required pair-interaction potentials, U(x), underpinning the experimental data by varying the input potential to the BD simulations. We thus generated simulated radial probability distribution functions, g(r)s, that matched the experimentally measurements. Validation and further discussion of the BD simulation set-up and approach are provided in Supplementary Information, section 2.1. Example input and necessary simulation files are provided in our Figshare repository, available at: https://doi.org/10.6084/m9.figshare.c.6132003.
We assumed a pairwise interaction potential U(x) of the form: \(U(x)=A{e}^{-{\kappa }_{1}x}+B{e}^{-{\kappa }_{2}x}+{U}_{{\rm{vdW}}}\). Here the first term represents the overall repulsive electrostatic free energy of interaction, \(\varDelta {F}_{{\rm{el}}}(x)=A\,\exp (-{\kappa }_{1}x)\), with A > 0 and the second term, \(\varDelta {F}_{{\rm{int}}}(x)=B\,\exp (-{\kappa }_{2}x)\), denotes the free-energy contribution arising from interfacial solvation2. Note that \({\kappa }_{2} < {\kappa }_{1}\approx \kappa\). The third term represents the vdW attraction between silica particles in solution, for which we have used the expression in ref. 55.
BD simulations were performed by taking into account the experimentally determined polydispersity in particle size as shown in Supplementary Fig. 1. This implies that at the simulation level a variable particle radius is taken into account to the lowest level of approximation (that is, the interaction potential remains fixed and independent of the size of the particles, which would not be true in practice). Using a value of the Hamaker constant AH = 2.4 zJ we found that UvdW made a small contribution (≈ −0.5kBT) to the total interaction energy at large separations, x ≥ 0.2 μm, that is, for the majority of experiments in this work56,57. However, for experiments at higher salt concentrations (c0 ≈ 1 mM; Fig. 1d), the vdW interaction can make a more substantial contribution (≈ −1kBT) to the interaction at separations x ≈ 0.1 μm, and for this reason it was included in our expression for U(x) when modelling these measurements.
The experimentally measured g(r) curve provides an estimate of the the location of the minimum in the pair potential xmin. In Supplementary Information equation (6), the screening length \({\kappa }_{1}^{-1}={\kappa}^{-1}\), which is known from the measured salt concentration. We then use a trial value of the interaction free energy at the minimum, \(U\left({x}_{\min }\right)=w < 0\), to obtain initial values for the parameters A and B as inputs for the pair-interaction potential, U(x), using Supplementary Information equation (4), where we have taken κ2/κ1 ≈ 0.95, as suggested in ref. 28. We note, however, that this ratio is not a strict requirement and that we may also treat it as a free parameter which yields an alternate set of parameters A, B, κ1 and κ2 that can provide equally good qualitative agreement with the experimental data (see, for example, Supplementary Table 5).
Particle configurations for the BD simulations were initialized via random particle placement in a 200 × 200 μm2 simulation box that reproduced the experimental particle density (≈ 0.008 particles μm−2). The polydispersity of the simulated colloids was drawn from the manufacturer’s size distribution for each particle type, as shown in Supplementary Fig. 1. Periodic boundaries were applied in the x,y dimensions while the z dimension was held finite. The z coordinate of the colloids were fixed at a constant height throughout the simulation, ensuring a 2D system and thus mimicking the experiment. BD simulations were performed assuming that the interactions between the particles may be regarded as pairwise additive (discussed further in the main text and in Supplementary Information, section 3.4)
Convergence of the potential energy per particle in our BD simulations was monitored over time (Supplementary Fig. 5). Particle positions used for the calculation of the final simulated g(r)s were collected once the value of the potential energy reached a stationary value, after ≈ 30 min of simulation time in a run involving a strongly attractive U(x) of a well depth of several kBT (Supplementary Fig. 5). Agreement between the simulated and the experimental g(r)s was assessed for a trial input pair-interaction potential U(x) and the value of the well depth w adjusted in subsequent BD simulations if required, to obtain a final simulated best match with the experimental data. This procedure permitted us to infer the functional form of an underlying pair-interaction potential U(x) capable of capturing the experimentally measured g(r).
Molecular dynamics simulations of alcohols at interfaces
The excess electrical potential due to the orientation of solvent molecules at an interface, φ0 or φint, is required as an input to the interfacial solvation model to calculate theoretical Utot(x) curves (Supplementary Information, section 3.2). To estimate φint(σ) as a function of surface electrical charge density, σ, we performed molecular dynamics (MD) simulations with the GROMACS MD code58. We examined the behaviour of a solvent phase in contact with a model surface composed of oxygen atoms in a parallel-plate capacitor set-up, as described extensively in previous work27,29. Example input files, force-field parameters and code for the analysis of the simulations performed in this study are available in our Figshare repository: https://doi.org/10.6084/m9.figshare.c.6132003.
Prior to running MD simulations in the capacitor set-up, we first ran preliminary simulations of a box of 7,500 IPA molecules, without the capacitor wall atoms, under constant pressure, maintained with the Parrinello–Rahman pressure-coupling method. The length of the box in z was allowed to fluctuate while keeping the x,y dimensions fixed to those of the capacitor walls of fixed area. This equilibrated slab of solvent was then sandwiched between capacitor plates comprised of positionally restrained oxygen atoms that only support Lennard–Jones interactions (Supplementary Fig. 7). In our simulations, IPA molecules were parametrized with the CHARMM36 force field59. As in previous work, the plates are ≈ 10×10 nm2 in area and are separated by ≈ 8 nm of solvent medium in the z direction27,29. This ensures that any oscillations in the solvent density or dipole moment profiles attain bulk-like properties at a location zmid in the middle of the capacitor. A subset of the atoms belonging to the first layer in each wall (in direct contact with the solvent) was randomly assigned a positive (left plate) or a negative charge (right plate) to generate an electric field of specific strength in the box while maintaining electroneutrality within the box. The capacitor system simultaneously yields estimates of φint(σ) for both positive and negative values of σ and provides a well-defined system for comparing solvation at a macroscopic surface with a continuum electrostatics model.
Next, a second round of equilibration was carried out for the entire capacitor system, including the capacitor walls, and consisted of a short NVT run with a velocity-rescaling thermostat, followed by 500 ps in an NPT ensemble where only the z dimension of the box was allowed to fluctuate, keeping the the x,y dimensions fixed. This ensured that solvent molecules were maintained at the correct density throughout the simulation box. Following this procedure, production MD runs of 20 ns duration were performed in an NVT ensemble with trajectory frames written every 20 ps. The particle mesh Ewald method was used to evaluate the long-range electrostatic interactions using a 1 Å grid spacing and a short-range cut-off of 12 Å. The Lennard–Jones interactions were smoothed over the range of 10–12 Å using the force-based switching function. We scaled the z dimension of the box by a factor 2 for Ewald summation only and applied the 3dc correction of Yeh and Berkowitz to remove artificial polarization induced by neighbouring image dipoles60. The orientation of solvent molecules as a function of distance from the walls was analysed to yield φint(σ) values that were used in the calculation of the interfacial free energy term, ΔFint, as previously described27,28,29.